R Maurin

J Chauchat

B Chareyre

P Frey

A minimal coupled fluid-discrete element model for bedload transport

A minimal Lagragian two-phase model to study turbulent bedload transport focusing on the granular phase is presented, and validated with experiments. The model intends to describe bedload transport of massive particles in fully rough flows at relatively low Shields numbers, for which no suspension occurs. A discrete element method for the granular phase is coupled with a one dimensional volume-averaged two-phase momentum equation for the fluid phase. The coupling between the discrete granular phase and the continuous fluid phase is discussed, and a consistent averaging formulation adapted to bedload transport is introduced. An original simple discrete random walk model is proposed to account for the fluid velocity fluctuations. The model is compared with experiments considering both classical sediment transport rate as a function of the Shields number, and depth profiles of solid velocity, volume fraction, and transport rate density, from existing bedload transport experiments in inclined flume. The results successfully reproduce the classical 3/2 power law, and more importantly describe well the depth profiles of the granular phase, showing that the model is able to reproduce the particle scale mechanisms. From a sensitivity analysis, it is shown that the fluctuation model allows to reproduce a realistic critical Shields number, and that the influence of the granular parameters on the macroscopic results are weak. Nevertheless, the analysis of the corresponding depth profiles reveals an evolution of the depth structure of the granular phase with varying restitution and friction coefficients, which denotes the non-trivial underlying physical mechanisms.

I. INTRODUCTION

Historically studied by hydraulic engineers in relation to the management of river waterways [START_REF] Duboys | Le rhône et les rivières à lit affouillable[END_REF] , bedload represents the main contribution of sediment transport to the evolution of riverbeds. As such, it has major implications for environmental flows and associated risks like floods for example. In contrast to suspension, bedload transport is characterized by sediment transport for which the vertical gravity force is on average stronger than the upward fluid force, i.e. sediments rolling, sliding or in saltation over the bed. The paper focuses on bedload transport in turbulent flow conditions, which is the most common case in nature.

In this phenomenon, one of the main challenges is to link the sediment transport rate to the fluid flow rate. By making the problem dimensionless, it is equivalent to linking (ρ p -ρ f )gd ; where Q s is the sediment transport rate per unit width, d is the particle diameter, ρ p and ρ f are respectively the particle and fluid density, and g is the acceleration of gravity. The usual semi-empirical formulas established in this framework such as the Meyer-Peter and Müller 2 one, can differ by two orders of magnitude from what is observed in the field [START_REF] Bathurst | Effect of coarse surface layer on bed-load transport[END_REF] . This difference is usually explained by the difficulty of measurements, the complexity of the physical processes and the great variability of the situations encountered in the field (e.g. grain size segregation, particle shape, channel geometry) [START_REF] Recking | A comparison between flume and field bed load transport data and consequences for surface-based bed load transport prediction[END_REF] .

Bedload transport can be viewed as a granular medium in interaction with a fluid flow.

Following this two-phase decomposition, there are two major possibilities for numerical modelling: a continuous description for the two phases (Euler/Euler) or a continuous description for the fluid phase and a discrete one for the granular phase (Euler/Lagrange).

The former considers the momentum conservation of the two phases viewed as two continua in interaction, and is based on the two-phase Reynolds Averaged Navier-Stokes (RANS) equations [START_REF] Jackson | The dynamics of fluidized particles[END_REF][START_REF] Drew | Theory of Multicomponent Fluids[END_REF] . The averaged equations require different closures, and the main differences between the models pertain to the Reynolds stress tensor and the constitutive law for the intergranular stress. The Reynolds stress tensor models the effect of turbulence on the mean fluid flow, and ranges from simple descriptions such as mixing length formulations, to more complex ones such as k -. In the case of intense bedload transport, also termed sheet flow, a substantial number of particle layers are in motion. The Euler/Euler description has therefore been mainly used for this regime, with closures for the granular stress tensor according to the main theories for granular media, i.e. Bagnold formulation [START_REF] Hanes | A granular-fluid model for steady intense bed-load transport[END_REF] , the µ(I) rheology [START_REF] Revil-Baudard | A two-phase model for sheet flow regime based on dense granular flow rheology[END_REF][START_REF] Aussillous | Investigation of the mobile granular layer in bedload transport by laminar shearing flows[END_REF] , or the kinetic theory [START_REF] Jenkins | Collisional sheet flows of sediment driven by a turbulent fluid[END_REF][START_REF] Hsu | On two-phase sediment transport: sheet flow of massive particles[END_REF] .

The continuous approximation breaks down for the granular phase when considering bedload transport closer to the threshold of motion which is the common situation in mountain streams and the focus of this paper. Moreover, the constitutive equation for granular media is still a matter of debate and thus limits the analysis of the results of such models.

Euler/Lagrange models overcome these two limitations by resolving the motion of each grain. For bedload transport, the high concentration of particles inside the bed requires modelling the contact between grains, this is today commonly handled with the Discrete Element Method (DEM). The different scales of fluid description range from large scale average description, to Direct Numerical Simulation (DNS) resolving the fluid locally around the particles down to the smallest turbulence length scale. Euler/Lagrange approaches have been intensively developed in recent years for problems with particles in fluids such as fluidized bed [START_REF] Tsuji | Discrete particle simulation of two-dimensional fluidized bed[END_REF][START_REF] Zhu | Discrete particle simulation of particulate systems: Theoretical developments[END_REF] , particle suspension [START_REF] Kidanemariam | Direct numerical simulation of horizontal open channel flow with finite-size, heavy particles at low solid volume fraction[END_REF] , or sheet flow [START_REF] Jiang | Multiparticle simulation methods applied to the micromechanics of bed load transport[END_REF][START_REF] Yeganeh | Applicability of euler-lagrange coupling multiphaseflow model to bed-load transport under high bottom shear[END_REF][START_REF] Drake | Discrete particle model for sheet flow sediment transport in the nearshore[END_REF][START_REF] Calantoni | Simulation of sediment motions using a discrete particle model in the inner surf and swash-zones[END_REF][START_REF] Yeganeh-Bakhtiary | A three-dimensional distinct element model for bed-load transport[END_REF][START_REF] Kidanemariam | Interface-resolved direct numerical simulation of the erosion of a sediment bed sheared by laminar channel flow[END_REF] . Focusing on bedload transport, up to now only a few contributions have taken advantage of the Eulerian/Lagrangian approach. The work of Duran et al. [START_REF] Duran | Numerical simulation of turbulent sediment transport, from bed load to saltation[END_REF] used a simple average description of the fluid with a two dimensional discrete element method for the particles, to numerically study the dependence of sediment transport on the solid-fluid density ratio. Bedload transport was considered in this paper as an extreme case of low density ratio, the closures of the model being more adapted to aeolian transport. Using a DNS/DEM model, Ji et al. [START_REF] Ji | Direct numerical simulation of sediment entrainment in turbulent channel flow[END_REF] focused on the influence of particle transport on the turbulence. With a similar model, Fukuoka et al. [START_REF] Fukuoka | Effects of sizes and shapes of gravel particles on sediment transports and bed variations in a numerical movable-bed channel[END_REF] studied particle shape influence and size-segregation effects.

Bedload transport has mainly been studied focusing on the fluid phase. It is however clear that the granular behavior is important in this phenomenon and should be studied further [START_REF] Frey | How river beds move[END_REF][START_REF] Frey | Bedload : a granular phenomenon[END_REF] . The idea is therefore to analyze bedload transport at the particle scale in order to understand the behavior of the bed as a granular medium. Considering the complexity of the experimental technique for particle-scale three dimensional (3D) bedload transport analysis (e.g. index matching [START_REF] Aussillous | Investigation of the mobile granular layer in bedload transport by laminar shearing flows[END_REF] , or Magnetic Resonance Imaging 26 ), there are interests in developing a numerical approach of the problem. Focusing on the granular phase, the paper presents a model for bedload transport using a DEM Lagrangian description of the granular phase coupled with a one dimensional volume-averaged two-phase momentum equation for the fluid phase. Although this type of model is common, to our knowledge, no previous contribution focused on bedload transport at relatively low Shields number. Moreover, the usual experimental validations are limited to the classical macroscopic results of dimensionless transport rate as a function of the Shields number. In this paper, we present a model adapted to subaqueous bedload transport (section II) and perform a new particle-scale experimental comparison with solid depth profiles in quasi-2D bedload transport cases 27 (section III). In addition, the classical experimental comparison of the sediment transport rate as a function of the Shields number is considered in a more general 3D framework (section IV). The influence of the different model contributions is considered in terms of sediment transport rate and solid depth profiles.

II. NUMERICAL MODEL FORMULATION

The proposed model is based on a DEM Lagrangian description for the solid phase and an Eulerian description for the fluid phase. In the present approach, the fluid flow is not solved at the particle scale and the momentum coupling is ensured in an averaged sense via semi-empirical correlations. After presenting briefly the Discrete Element Method (section II A) and the fluid phase description (section II B), the coupling between both phases is discussed by detailing in particular the averaging procedure and the coupling forces in the framework of bedload transport (section II C).

A. Solid phase

The DEM, originally introduced by Cundall & Strack 28 for granular media, is based on the explicit resolution of Newton's equation of motion for each individual particle considering nearest neighbor contact forces f p c . For each particle p at position x p the equation of motion reads:

m d 2 x p dt 2 = f p c + f p g + f p f , (1)
where f p g is the force due to gravity and f p f represents the forces applied by the fluid on particle p. This last term arises from the DEM-fluid coupling and will be detailed in subsection II C. The application of the gravity force is straightforward. The contact forces are determined from the relative displacement of the neighboring particles using a defined contact law. In bedload transport, there is a sharp transition between rapidly sheared particles at the interface with the fluid, and almost quasi-static motion inside the bed. The so-called spring-dashpot contact law used in this paper, allows description of these two types of behavior and is classical in granular flows modelling. The contact law is based on a spring of stiffness k n in parallel with a viscous damper of coefficient c n for the normal contact, coupled with a spring of stiffness k s associated with a slider of friction coefficient µ for the tangential contact. For normal contact, the linear elastic spring and viscous damping define a constant restitution coefficient e n characteristic of the energy loss at collision, which can be evaluated experimentally.

B. Fluid phase

The fluid phase model is based on spatially averaged two-phase Navier-Stokes equations [START_REF] Jackson | The dynamics of fluidized particles[END_REF] , and inspired from the one-dimensional Euler-Euler model proposed by Revil-Baudard & Chauchat [START_REF] Revil-Baudard | A two-phase model for sheet flow regime based on dense granular flow rheology[END_REF] to deal with turbulent unidirectional sheet-flows. The simplifications of the general fluid phase equations 5 due to the incompressible and unidirectional character of the present bulk flow lead to the resolution of the same fluid phase momentum equation: In equation ( 2), omitting the fluid-particle interaction term, closure laws for the viscous shear stress τ xz f and the Reynolds shear stress R f xz need to be prescribed. In the present model, the fluid is considered as Newtonian, so that:

ρ f ∂ u x f ∂t = d ( τ xz f ) dz + dR f xz dz + ρ f g sin α -n f x s , (2)
τ xz f = ρ f ν f d u x f dz , (3)
where ν f is the clear fluid kinematic viscosity.

The Reynolds shear stress, representing the vertical turbulent mixing of horizontal momentum, is modeled based on the eddy viscosity concept (ν t ) with a mixing length formulation:

R f xz = ρ f ν t d u x f dz with ν t = l 2 m d u x f dz , (4)
in which the mixing length l m formulation proposed by Li & Sawamoto 29 is used:

l m (z) = κ z 0 φ max -φ(ζ) φ max dζ, (5)
where κ = 0.41 represents the von Karman constant. This simple formulation allows recovery of the two expected asymptotic behaviors: the mixing length is linear with z when the solid phase volume fraction vanishes (i.e. clear fluid), as in the law of the wall [START_REF] Prandtl | Bericht über neuere Turbulenzforschung[END_REF] , and the mixing length is zero when the solid phase is at its maximum packing fraction, i.e. the turbulence is fully damped inside the dense sediment bed. As explained in Revil-Baudard & Chauchat 8 , this formulation is well adapted for boundary layer flow above mobile rough beds. Indeed, the integral of the solid volume fraction predicts a non zero mixing length at the transition between the granular dominated and turbulent dominated layers. Also, with this formulation no virtual origin for the mixing length has to be prescribed.

C. DEM-fluid coupling

The key point in continuous-discrete models consists in the coupling of the two phases, which involves an averaging procedure and a parametrization of the fluid forces applied on the particles.

Averaging procedure

For this purpose, the spatial averaging operator for the solid phase needs to be defined consistently with the spatial averaging operator for the fluid phase [START_REF] Jackson | The dynamics of fluidized particles[END_REF][START_REF] Anderson | Fluid mechanical description of fluidized beds. equations of motion[END_REF] . According to Jackson [START_REF] Jackson | The dynamics of fluidized particles[END_REF] , the solid phase volume fraction φ( x 1 ) at a given position x 1 is defined as:

φ( x 1 ) = p Vp G(| x 1 -x|)dV (6)
where the sum is over all the particles, V p is the volume of particle p, and G( x) is a weighting function that must be normalized on the whole physical domain. Providing that the weighting function G is defined, the solid phase spatial average of a scalar quantity γ at a given position x 1 is defined as:

γ s = 1 φ( x 1 ) p Vp γ( x)G(| x 1 -x|)dV, (7)
In the general case, Jackson 5 uses a radial weighting function G. However in the present case, to match the discretization of the fluid resolution it is more convenient to define a cuboid weighting function. To fulfill the normalization property, a three-dimensional step function is chosen for the weighting function:

G( x) =      1 l x l y l z for |x| ≤ l x /2, |y| ≤ l y /2, |z| ≤ l z /2 0 otherwise (8)
In order to properly define the spatial averaging, the average should be independent from the length scales chosen for the weighting function: l x , l y and l z 5,31 . This is only possible if a separation of scales exists between the macroscopic length scale of the problem L, the length scales associated with the weighting function l x , l y , l z and the particle diameter d, i.e.

L >> l x , l y , l z >> d.

Due to the sharp transition occurring at the sediment bed interface in the wall-normal direction, the wall-normal macroscopic length scale of the problem L is lower than the particle diameter d. Therefore the vertical length scale of the weighting function l z should be lower than the particle diameter in order to accurately resolve the vertical gradients of the averaged solid phase variables. We postulate that this limited vertical length scale l z , can be compensated statistically by larger complementary horizontal length scales, l x and l y . The convergence analysis of the numerical results on the length scale l x presented in appendix A lends credibility to this hypothesis.

Fluid forces

The force applied by the fluid on a single particle f p f introduced in equation ( 1) is defined as the integral of the total fluid stress, pressure and shear stress, acting on the particle surface [START_REF] Jackson | The dynamics of fluidized particles[END_REF] . In the present model, the fluid flow is not resolved at the particle scale so that the hydrodynamic forces cannot be computed explicitly, and need to be prescribed through semi-empirical formulas based on average fluid variables. The main hydrodynamic forces in bedload transport reduce to the buoyancy, the drag and the shear-induced lift. Ji et al. [START_REF] Ji | Direct numerical simulation of sediment entrainment in turbulent channel flow[END_REF] numerical results exhibit a non-negligible importance of the lift force with respect to the other two. However, Schmeeckle et al. [START_REF] Schmeeckle | Forces on stationary particles in near-bed turbulent flows[END_REF] showed experimentally that the usual formulation of the lift 33 , derived using the inviscid flow assumption, is not valid close to the threshold of motion. Based on this observation and the absence of alternative formulation, it has been decided not to include the lift force at this stage.

Therefore, the force f p f induced by the fluid on a particle p appearing in the DEM model (equation ( 1)), reduces to buoyancy f p b and drag f p D :

f p f = f p b + f p D . (9)
According to Jackson 5 the generalized buoyancy force is defined as:

f p b = V p -∇ P f + ∇. τ f , (10)
where P f is the average fluid pressure and τ f is the average viscous shear stress tensor taken at a larger scale than the particle scale. This definition generalizes the so-called Archimedes buoyancy force for hydrostatic problems to cases where the fluid volume is submitted to a macroscopic deformation at a scale much larger than the particle scale i.e.

the fluid deformation viewed by the particles can be considered as constant. Similarly to Revil-Baudard & Chauchat 8 , we found that the viscous stress tensor contribution is however negligible with respect to the pressure contribution in bedload transport. The force applied on each particle can then be approximated by the usual buoyancy expression, which is equivalent to apply the buoyant weight in the vertical direction.

The drag force exerted by the fluid flow on a single particle is classically expressed as:

f p D = 1 2 ρ f πd 2 4 C D u f x p -v p u f x p -v p , (11)
where C D is the drag coefficient, and u f x p -v p is the relative velocity between the particle and the average fluid velocity taken at the position of the particle center. The Dallavalle formulation [START_REF] Dallavalle | Micrometrics : The technology of fine particles[END_REF] together with a Richardson-Zaki correction [START_REF] Richardson | Sedimentation and fluidization: Part i[END_REF] is used in the present model for the drag coefficient:

C D = 0.4 + 24.4 Re p (1 -φ) -ζ , (12)
where Re p = || u f x p -v p ||d/ν f is the particulate Reynolds number for particle p. This simple formulation has been used in different two-phase flow models for sediment transport applications [START_REF] Revil-Baudard | A two-phase model for sheet flow regime based on dense granular flow rheology[END_REF][START_REF] Jenkins | Collisional sheet flows of sediment driven by a turbulent fluid[END_REF][START_REF] Hsu | On two-phase sediment transport: sheet flow of massive particles[END_REF] . The Richardson-Zaki correction (1 -φ) -ζ accounts for the hindrance effect induced by the local particle concentration, and allows to recover realistic fluid velocity in the particle bed. The exponent has been set to ζ = 3.1 in reference to Jenkins & Hanes [START_REF] Jenkins | Collisional sheet flows of sediment driven by a turbulent fluid[END_REF] . Equations ( 11) and ( 12) are used to compute the drag force on each individual particles in the DEM model (equation ( 1)), while the effect of buoyancy is taken into account through the vertical buoyant weight.

The average effect of the particles on the fluid momentum balance does not simply consist in the solid averaging of the momentum transfer associated with the hydrodynamic forces. It also includes higher-order correlations which appear in the averaging process, and are due to perturbations of the flow by the presence of the particles. For the case of Stokesian particles at low concentration, Jackson showed analytically 36 that these higher-order correlations lead to a modification of the viscosity in the average viscous fluid stress tensor formulation, which takes the form of Einstein's effective viscosity [START_REF] Einstein | Eine neue bestimmung der molekül dimensionen[END_REF] . In the model, the clear fluid viscosity in equation ( 3) has been replaced by Einstein's effective viscosity ν e to take this effect into account:

ν e = ν f 1 + 5 2 φ . ( 13
)
The phase interaction term in the fluid momentum balance (eq. 2) reduces then to the momentum transfer associated with the hydrodynamic forces. In the present 1D fluid resolution, it is expressed as the average number of particles n = φ/V p = 6φ/πd 3 multiplied by the solid-phase average streamwise associated force. For drag force, it gives:

n f D x s = 3 4 φ ρ f d C D u f x p -v p u x f -v p x s . (14)
The drag coefficient C D depends on the relative velocity through the particle Reynolds number, so that it should be included in the spatial averaging.

Velocity fluctuation model

The proposed average model for the fluid phase does not account for the fluid turbulent velocity fluctuations, which are known in particular to influence the particle threshold of

u f x u f z = - R f xz ρ f , (15)
where the • represents an averaging operator in time.

From experimental measurements in open-channel flows [START_REF] Nezu | Turbulent structure in open-channel flows[END_REF][START_REF] Nezu | Turbulence in Open Channel Flows[END_REF] , it has been observed that the magnitude of the fluctuations in the streamwise direction is roughly two times larger than in the vertical direction. With this constraint the following DRW model for the streamwise component (u f x ) p and the normal component (u f z ) p of the fluid velocity fluctuation associated with each particle p is proposed:

(u f z ) p = λ 1 (16a) (u f x ) p = -(u f z ) p + λ 2 , (16b)
where λ 1 and λ 2 are two Gaussian random numbers of zero mean and of standard deviation σ. This standard deviation is obtained from the local value of the Reynolds shear stress at the position of the particle center σ = R f xz ρ f ( x p ). The velocity fluctuations are updated every τ t , defined as the turbulent eddy turn over time, which can be estimated as τ t = w d /U f where w d is the water depth, and U f is the average fluid velocity. These velocity fluctuations are added to the average fluid velocity in the drag force expression (equation ( 11)).

D. Numerical resolution strategy

The resolution of the fluid equation still needs to be clarified. For numerical reasons, it is necessary to express equation ( 2) linearly as a function of the average fluid velocity. The numerical treatment of the drag force is then handled as follows:

n f D x = β u x f -v p x s , ( 17
)
where β is computed according to equation ( 14) as:

β = 3 4 φ ρ f d C D u f x p -v p u x f x p -v p x s u x f x p -v p x s (18)
This formulation allows strictly the same average momentum transfer in the discrete solid phase problem and the continuous fluid phase one. With this definition the fluid phase momentum equation to be solved can be rewritten as:

ρ f ∂ u x f ∂t = ρ f ∂ ∂z ν e + ν t ∂ u x f ∂z + ρ f g sin α -β u x f -v p x s (19)
This equation is discretized using implicit finite differences for the diffusion and the drag terms. The resulting tridiagonal system is solved using a double-sweep algorithm [START_REF] Chauchat | Modelling sedimentation-consolidation in the framework of a one-dimensional two-phase flow model[END_REF] .

The fluid phase resolution period τ f should be small enough compared with the particle relaxation time τ D = β -1 . This characteristic time corresponds to the time needed by a particle initially at rest to reach its steady state velocity in a constant fluid flow.

The DEM solid phase model is solved using the open-source code Yade [START_REF] Šmilauer | The Yade Project[END_REF] . The time integration is explicit with a second order centered scheme [START_REF] Bathe | Numerical methods in finite element analysis[END_REF] to ensure energy conservation.

The time step has been estimated with a method similar to Catalano 44 (pp. 84, see also [START_REF] Catalano | Pore-scale modeling of fluid-particles interaction and emerging poromechanical effects[END_REF] ), considering the rigidity of the system of springs [START_REF] Chareyre | Dynamic spar elements and discrete element methods in two dimensions for the modeling of soil-inclusion problems[END_REF] and dampers as decoupled.

III. EXPERIMENTAL COMPARISON

The model is to be compared against experimental data. The declared aim of the present model is to focus on the granular phase behavior. We therefore reproduce the quasi-2D experiments of Frey [START_REF] Frey | Particle velocity and concentration profiles in bedload experiments on a steep slope[END_REF] , in which particle tracking allowed to obtain average solid depth profiles of bedload transport. After a brief presentation of the experiment, the numerical set-up and the comparison with the experimental results will be presented.

A. Description of the experiment

The experiment of Frey

(F r = U f / √ gw d 1)
. The particle settling velocity (w s = 0.54 m/s) and the suspension number S * = w s /u * are high, meaning that the particles are weakly influenced by the turbulent structures. A camera is placed perpendicular to the sidewall, filming a window of 25x8 cm 2 at 131.2 frames per second. Due to the one particle diameter width of the channel, image processing [START_REF] Böhm | Two-dimensional motion of a set of particles in a free surface flow with image processing[END_REF] enables particle trajectories to be followed inside the measurement window, and the average free-surface elevation to be evaluated. In each experiment, once bedload transport is at equilibrium, data acquisition time lasted 60 s. Experimental data are averaged in the same way as in the model using the definition of section II C. For more details on the experimental setup, refer to the original experimental article of Frey [START_REF] Frey | Particle velocity and concentration profiles in bedload experiments on a steep slope[END_REF] . The order of magnitude of the main dimensionless numbers associated with the experiment are shown in table I. The Stokes number comparing the inertia of the particle and the viscosity of the fluid is given by St = ρ p v p d/(9η f ).

B. Application to the model

To compare the model with the experiments, the simulation needs to match the experimental set-up. To focus on the bulk equilibrium properties of bedload transport, periodic boundary conditions are considered in the streamwise direction for the present 2D case.

The periodic characteristic of the granular phase, does not enable us to impose a feeding rate.

To have a situation equivalent to the experiment, the density of beads per unit length (equivalent to the number of layers of particle N l ) and the free-surface position h are instead imposed. Indeed, there is a unique couple, slope-water depth, corresponding to the trans- In the experiment the width to depth ratio is low, and we expect in consequence the fluid flow to have a complex 3D structure. However, experimental flow measurements in this particular channel showed that it still has a typical logarithmic profile [START_REF] Frey | Experimental study of narrow free-surface turbulent flows on steep slopes[END_REF] . This, together with the stated aim of the model to focus on the granular phase, made us consider only a correction for the fluid dissipation at the smooth lateral walls. The correction was included as a source term in the fluid averaged momentum balance resolution (eq. 19), taking the form of a dissipation term evaluated from the classical Einstein method with a Graf and Altinakar friction factor [START_REF] Graf | Fluvial hydraulics: flow and transport processes in channels of simple geometry[END_REF] . The method description can be found in Frey et al 51 . For each run the channel bottom is newly generated randomly, and particles are deposited under gravity. Once the system with fluid resolution is at equilibrium, the simulations last 100 seconds and measurements are made each 0.1 second. The latter corresponds to the particle relaxation time to the fluid velocity τ D = β -1 (eq. 18), and is characteristic from the evolution of the system. For the post-processing of both experimental and numerical results, the averaging definition is taken consistently with the numerical resolution from equation (A1).

To study the stability of the coupling, we performed a sensitivity analysis on the fluid resolution period, shown in appendix B, and set it to τ f = 10 -2 s regarding the results obtained. In agreement with Revil-Baudard & Chauchat [START_REF] Revil-Baudard | A two-phase model for sheet flow regime based on dense granular flow rheology[END_REF] , it has been found that the fluid effective rheology does not influence the fluid behavior as it is dominated by the turbulent shear stress. We therefore used a clear fluid viscosity. The restitution coefficient was set to e n = 0.5 based on measurements previously made in the experimental channel considered [START_REF] Bigillon | Etude du mouvement bidimensionnel d'une particule dans un courant d'eau sur forte pente[END_REF] .

In the present situation, the limited particle pressure allows artificial softening of the spheres stiffness in order to reduce computational costs. k n was set to 5.10 3 N/m which leads to an acceptable average overlap of the order of 10 -4 d and allows to be in the rigid grain limit [START_REF] Roux | Quasistatic rheology and the origins of strain[END_REF] .

The

C. Results

In bedload transport, one of the main challenges lies in the prediction of the integrated transport rate as a function of the flow rate. The experiment of Frey [START_REF] Frey | Particle velocity and concentration profiles in bedload experiments on a steep slope[END_REF] was designed to give more insight into the granular behavior, and to focus on the depth profile of bedload transport at the particle scale. It has been noted previously [START_REF] Frey | Particle velocity and concentration profiles in bedload experiments on a steep slope[END_REF][START_REF] Lajeunesse | Bed load transport in turbulent flow at the grain scale: Experiments and modeling[END_REF] that the integrated transport rate per unit width Q s can be expressed as a function of the average transport rate density q s , the product of the average solid velocity v p and solid volume fraction φ:

Q s = q s s dz = v p s φdz. (20)
Considering bulk equilibrium properties of bedload transport, v p s , φ, and q s s depends only on the depth z. The experimental comparison will then focus on the depth profiles of the solid volume fraction, the average solid velocity, and the transport rate density, which will be called for simplicity transport rate in the following. To complete this decomposition, we evaluated also the Shields number for each case. This was done with two different methods: from macroscopic parameters following Frey

* = ρ f u 2 * /[(ρ p -ρ f )gd],
where u * is given by the maximum turbulent shear stress u * = max(R f xz (z)). θ * was evaluated only in the simulation. This formulation avoids use of the macroscopic determination of the Shields numbers, which is sensitive to the water depth evaluation and the type of wall correction used.

In the previous subsections, we did not introduce any experimental or numerical error.

It appears that the dispersion is dominated by the limited measurement window length of the experiment (40 d). The order of magnitude of this dispersion has been evaluated numerically. Figure 3 exhibits the depth profiles of the solid volume fraction, the solid streamwise velocity, and the solid transport rate for a single simulation with different postprocessing averaging properties. The fluid mechanics convention is used, where the depth is represented on the y-axis while the quantities of interest are represented on the x-axis. The simulation corresponds to the case Sim20 in table III, considering a periodic length cell of l x = 1000 d. The figure shows the variability of the results when the averaging cell length is taken equal to the experimental one at different position in the channel. This dispersion is much greater than the evaluated experimental dispersion and the numerical variability due to the size of the periodic cell simulated. These latter two will be consequently ignored in the comparison, and the variability observed will be taken as error bars.

The three different experiments detailed in table III are considered for experimental comparison. The slope is the same and equilibrium transport rate ranges from 6 to 20 beads/s.

The differences in the input parameters between the runs lie in the water surface position h and the number of layers of particle N l . The different experimental cases represent a good test to evaluate the sensitivity to the parameters and the ability of the model to reproduce different experimental conditions. The macroscopic results presented in table III show that the integrated transport rates ṅ are in good agreement with the experiment even if slightly overestimated. Considering the Shields number, the two different methods of evaluation lead to an over-estimation using the macroscopic formulation θ, and an under-estimation using the formulation based on the turbulent shear stress θ * . This underlines the complexity to measure the Shields number in the present case, especially when using the macroscopic formulation which is very sensitive to the small water depth. The trends observed with both formulations are good, and the values have the same order of magnitude than the experiment. In the following, we will use θ * in order to avoid the somehow arbitrary determination of the water depth. Using this definition, the value observed for case 6 is below the classical critical Shields number. It should however be kept in mind that the present quasi-2D mono-disperse bed is less resistant, and the value of the critical Shields number is accordingly lowered. To summarize, the general trends observed for the macroscopic parameters are good and these results show that we are able to reproduce the experimental sensitivity to the free surface position and the number of bead layers.

To go further, figure 4 shows the solid depth profiles of velocity, volume fraction and transport rate density, for the three experimental comparisons. The global trends from one case to the other are well reproduced, and the shape of the simulated curves are close to the experimental ones. Focusing on the transport rate density profile, for each case the value of the peak is slightly overestimated, while the rest of the curve is in very good considered with respect to the average value around which it is oscillating. For the solid volume fraction profile, the agreement between simulation and experiments is excellent for case 6 and 20, while we note a slight discrepancy for case 14 at the interface. The solid velocity profiles show a good estimation of the maximal velocity, and of the depth structure.

The underestimation of the sediment transport rate peak is shown to correspond to an overestimation of the solid volume fraction in case 14, and of solid velocity in case 6 and 20.

For completeness, the fluid velocity profiles are presented in figure 5. No experimental data are available for comparison so that the simulated solid velocity profiles have been added

for reference. In the fixed bed, the fluid velocity exhibits some oscillations around a finite constant value and the solid velocity is negligible. The oscillations are due to the layering observed in the solid volume fraction profile (figure 4) that makes the drag coefficient oscillate accordingly (eq. 12-14). In the upper part, the velocity difference between the solid and the fluid phases is of the order of an isolated particle settling velocity w s ∼ 0.54m/s. The results are consistent with the drag coefficient formulation adopted and cannot be interpreted showed that the model is able to reproduce the particle-scale trends observed experimentally, and the variation between the three different runs.

IV. DISCUSSION

The experimental comparison gives credits to the model presented, and shows that the depth structure of the phenomenon is well reproduced. Starting from this point, the effect of the grains parameters (restitution and friction coefficients) and the fluid velocity fluctuations model are analyzed over a wide range of Shields numbers, in terms of dimensionless sediment transport rate versus Shields number, completed by solid depth profiles. The analysis aims at characterizing the influence of these parameters, in order to both study the influence of the different terms on the phenomenon, and the robustness of the experimental comparison. To extend the generality, a 3D bi-periodic (streamwise, spanwise) situation is considered. The random fixed bottom is generated from a gravity deposition, fixing all the particles contained in a slice of height d at a given elevation in the granular bed. The size of the periodic cell has been chosen from a convergence analysis similar to the one undertaken for the 2D case (see appendix A) and a cell size of l x = l y = 30d has been chosen to ensure statistical convergence and numerical stability. For each run, the DEM results are averaged over 100 seconds.

A. Macroscopic considerations

The dimensionless sediment transport rate as a function of the Shields number is presented in figure 6. The model results are compared with experimental data from Meyer-Peter and Müller 2 (+), and Wilson 55 (x). Simulation parameters for the reference configuration, represented as black squares ( ), are the same as the one used for the experimental comparison (see table II). The 3/2 power law is recovered by the numerical simulations and the results show a good agreement with experimental data for Shields number θ > 0.1. Near the threshold of motion, the model results differ from experimental measurements. The linear inset shows that the transition around the critical Shields number, characteristic of the onset of motion, is sharper in the numerical simulation results than in the experimen- Focusing on the linear plot, it is observed that the critical Shields number is changed from around 0.04 to around 0.09 in the case without fluid velocity fluctuations. The former is in the range of observed values under turbulent flow conditions [START_REF] Buffington | A systematic analysis of eight decades of incipient motion studies, with special reference to gravel-bedded rivers[END_REF] , while the latter is close to the value observed under laminar flow conditions [START_REF] Ouriemi | Determination of the critical shields number for particle erosion in laminar flow[END_REF] . The influence of the fluctuations on the critical Shields number can be associated with turbulent coherent structures (e.g. [START_REF] Papanicolaou | Stochastic incipient motion criterion for spheres under various bed packing conditions[END_REF][START_REF] Dwivedi | Role of turbulence and particle exposure on entrainment of large spherical particles in flows with low relative submergence[END_REF] ).

However, the present simple fluid velocity fluctuation model does not account for the spacetime correlations induced by turbulent boundary layer coherent structures. This partly explains that the fluctuations model does not allow to describe well the evolution of the sediment transport rate with Shields number close to the threshold of motion (figure 6).

Nevertheless it permits to successfully reproduce the onset of sediment transport motion in the turbulent regime, resulting in a good comparison with experimental depth profiles (section III C).

In the rigid grains limit, the granular interactions are characterized by the restitution and the friction coefficients. The restitution coefficient is representative of the energy loss during collisions, and has been shown experimentally to decrease with the impact velocity following a power law exponent lower than 1/4 [START_REF] Gondret | Bouncing motion of spherical particles in fluids[END_REF] . For a limited range of impact velocity and in a first approximation it can be considered as constant. As the fluid flow model does not allow to resolve the fluid at the particle scale, the local lubrication effect is included in the effective restitution coefficient e n . Following Gondret et al. [START_REF] Gondret | Bouncing motion of spherical particles in fluids[END_REF] , the effective restitution coefficient depends on the local Stokes number comparing the grain inertia to fluid viscous forces: St = ρ p vd/(9η f ) where η f is the fluid dynamic viscosity, and v is the impact velocity.

In the region where the collisions are dominant, the Stokes number is of order 10 2 -10 3 , corresponding to effective restitution coefficient in the range e n ∈ [0.6e d n -0.9e d n ], respectively, where e d n is the restitution coefficient for dry grains [START_REF] Gondret | Bouncing motion of spherical particles in fluids[END_REF] . In our model a constant restitution coefficient is adopted, taking into account the lubrication effect globally. It is therefore a characteristic of the material and of the lubrication effect.

The influence of the restitution coefficient is shown in figure 6 and no dissipation at contact (e n = 1). Focusing on the realistic range at high Shields number (θ ∼ 0.45), the effect on the sediment transport rate is negligible. A slight trend is observed, the sediment transport rate and the Shields number being respectively increasing and decreasing function of the restitution coefficient. The extreme case without dissipation at contact (e n = 1) follows the same trend but exhibits a more important transport rate increase. Quite surprisingly, the case e n = 0.01 shows an increase in transport rate with respect to case e n = 0.25. For the lower Shields number value (θ ∼ 0.1), while the dependency in restitution coefficient is limited in the realistic range, there is no associated clear trend.

The non-monotonous dependencies observed show non-trivial coupling between the granular phase characteristics and the sediment transport rate. The global weak dependency on the restitution coefficient is consistent with results obtained by Drake & Calantoni 17 under oscillatory flow conditions, and show that there is no need to include a lubrication model in the present condition. However, the relatively low importance of the restitution coefficient at such a high Shields number value is surprising. It is usually thought that collisional interactions are the dominant mechanism of momentum diffusion for such inertial particles [START_REF] Jenkins | Collisional sheet flows of sediment driven by a turbulent fluid[END_REF][START_REF] Armanini | Rheological stratification in experimental free-surface flows of granular liquid mixtures[END_REF] .

The effect of the particle friction coefficient is also shown in figure 6, represented by circles: µ = 0.2 (•), µ = 0.6 (•), µ = 0.8 (•). Unlike for the restitution coefficient, the trend observed is monotonous for all values, and similar at low (θ ∼ 0.09) and high Shields number (θ ∼ 0.35). The sediment transport rate and the Shields number decrease with increasing friction coefficient. The effect appears to be non-linear as the observed influence for a variation from µ = 0.2 to 0.4 is much greater than the one observed for a variation from µ = 0.4 and 0.8. This type of dependency is characteristic of dry dense granular flows [START_REF] Da Cruz | Rheophysics of dense granular materials: Discrete simulation of plane shear flows[END_REF] .

As a partial conclusion, the present analysis shows that (i) the 3/2 power law for the sediment transport versus Shields number relationship is well captured by the proposed model; (ii) the fluid velocity fluctuations model is essential to capture a realistic value for the critical Shields number under turbulent flow conditions (iii) the influence of the granular interaction parameters is low, when taken in a realistic range. These results underline the robustness of the model and strengthen the experimental validation. Extreme value of particle friction and restitution coefficient affects the results, and show complex behaviors.

In order to understand better the mechanisms at work, the sensitivity to granular interaction parameters will be further discussed by analyzing the results in terms of depth profiles.

B. Depth profiles analysis

Figure 7a shows the solid depth profiles for a Shields number value θ ∼ 0.45, and for the different restitution coefficient values. Such a high Shields number value enhances the effect of restitution coefficient. From the profile, a clear trend appears, the sediment transport rate density profile is broader with increasing restitution coefficient. Excepting case e n = 0.01, this is associated with an overall increase of particle velocity throughout the depth, the velocity profile being shifted with almost the same shape. The solid volume fraction profile shows an increase of the mobile layer thickness when the restitution coefficient is increased: the solid volume fraction is lowered close to the quasi-static bed while it is increased in the upper part of the flow. This can be explained considering predictions of the kinetic theory of granular flows, where particle phase normal stress is an increasing function of the restitution coefficient 10 . The increase in particle normal stress is logically associated with a decompaction of the bed, which is submitted only to gravity. The case e n = 0.01 is peculiar, the global trend is observed in the lower part of the flow while a higher particle velocity is predicted in the upper part. This reflects the coupling with the fluid phase, and features the complex mechanisms at work.

Figure 7b shows the influence at high Shields number θ ∼ 0.45 of the friction coefficient over the range µ ∈ [0.2; 0.8] with 0.2 steps. Interestingly, the solid volume fraction profile is not affected by the variation in friction coefficient. On the contrary, the particle velocity and thus the sediment transport rate density profiles are increasing when the friction coefficient is decreased. The increase of the velocity throughout the depth is mainly affecting the lower part of the sediment transport rate density profile, where the solid volume fraction is maximum. It corresponds to the denser part of the granular flow, for which the frictional interactions are dominant.

This analysis shows that, while affecting weakly the macroscopical results, the friction and restitution coefficient impact the depth structure of the granular flow differently. In addition, the non-monotonous behavior observed suggests the presence of non-trivial coupling between the solid and the fluid phases.

V. CONCLUSIONS

The model presented is a step toward a description of the granular processes of steady bedload transport. By adapting the closures to this particular case, it has been shown that focus on steady equilibrium results, we consider time-averaged results which are converged in time. The convergence analysis is conducted with respect to a large reference cell size for which we consider that the results at each fluid resolution are spatially converged. Indeed, a convergence with respect to this case ensures that the error made in the spatial averaging along the simulation are compensating each other. The analysis focuses on the transport rate depth profile. For both 2D and 3D cases, l z is taken at its minimal value in the problem

l z = d/30.
quasi-2D convergence analysis We present here the results of the time-averaged spatial convergence analysis for the quasi-2D case Sim20 detailed in section III B. In this configuration l y is fixed at the channel width, and the problem is considered only as a function of the streamwise length l x which determines the size of the averaging cells. The convergence analysis is made with respect to the reference state chosen as l x = 10000d, corresponding to a periodic length cell of 60m for particles of 6mm and about 80000 particles in the simulation. We performed different simulations with a periodic streamwise length cell l x of respectively 50, 100, 200, 300, 500, 1000, 2500 and 5000 d.

To quantitatively analyze the transport rate profile differences, an indicator representative of the deviation with respect to a reference case is defined. It is given as the root mean square (RMS) of the difference between the considered transport rate profile and the reference one: This variable effectively measures how close the results considered is from the reference case. experimental comparison of the paper, was chosen as l x = 1000 d, as it gives the best trade off between computational time and deviation observed. The main results are summarized in figure 9, expressing the RMS deviation with respect to the reference case as a function of the product l x l y /d 2 . The latter reflects the statistical representativity as it directly determines the size of the averaging cell V = l x l y l z (l z = dz fixed). The figure shows that cases with l x = l y give better convergence than the ones with l x = l y . There does not seem to be a convergence with increasing cell area. To our opinion, this reflects the fact that the results are already converged. Interestingly, the relative periodic cell size for convergence is less important in the 3D configuration, as a cell size of l x = l y = 10d (to compare with l x = 1000d and l y = 1d) is almost already converged with respect to l x = l y = 100d. This can be explained by the better randomness of the 3D packing, and suggests that the statistical representativity was not the limiting parameter in the quasi-2D convergence analysis. For 3D cases, a higher Shields number increases the numerical instability of the coupling, so that it has been necessary to consider cell sizes up to l x = l y = 50d for the highest Shields numbers simulations presented in the paper.

∆Q rms i Q ref = 1 N N z=0 ( Q i z -Q ref z ) 2 1 N N z=0 Q ref z , (

  the dimensionless sediment transport rate Q * s = Qs d √ (ρ p /ρ f -1)gd , to the Shields number which compares the fluid bed shear stress τ f b to the buoyant weight θ = τ f b

  where is the fluid phase volume fraction, ρ f is the fluid density, u x f is the averaged fluid velocity, τ xz f is the averaged fluid viscous shear stress, R f xz is the Reynolds shear stress, α is the channel inclination angle, and n f x s is the averaged fluid-particle general interaction term. A schematic picture representing the main fluid model variables is shown in figure 1. The operator . s denotes a spatial averaging over the solid phase while the operator . f denotes a spatial averaging over the fluid phase. The major difference with the continuous two-phase model proposed by Revil-Baudard & Chauchat 8 and with Euler/Euler models in general, stands in the average fluid-particle interaction n f x s and solid volume fraction φ = 1 -, obtained in the present model from a spatial averaging of the DEM solution, whereas otherwise obtained by solving a continuous momentum balance equation. All the details concerning the averaging process and the fluid-particle interaction term will be given in subsection II C.

FIG. 1 .

 1 FIG. 1. Sketch of the problem representing the axes and the variables used in the model: the water free-surface position h, the water depth w d , the slope S 0 = tan α, the gravity vector g, as well as the profiles of average streamwise fluid and solid velocities (resp. u x f and v p x s ), and solid volume fraction φ. Streamwise periodic boundary conditions (BC) are used for the solid phase DEM description as indicated on the scheme.

  motion. In order to account for these turbulent processes in the average fluid model, a Discrete Random Walk (DRW) model for the fluid velocity fluctuations inspired from Zannetti 38 is therefore introduced. It consists in associating a random velocity fluctuation with each particle for a given duration, as a function of the local turbulent intensity and turbulent time scale. The fluctuations are not correlated in space, nor in time, and the model is built so that the Reynolds shear stress definition is consistent between the average fluid model and the DRW model:

  [START_REF] Frey | Particle velocity and concentration profiles in bedload experiments on a steep slope[END_REF] consisted in a quasi-2D ideal case of mountain stream bedload transport on a steep slope. The setup is depicted in figure2, it is composed of a 2 m long inclined channel of slope S 0 = 0.1, and width 6.5d/6. Water (ρ f = 1000 kg/m 3 ) flows inside the open-channel and entrains the spherical glass particles (ρ p = 2500 kg/m 3 ) of diameter d = 6 mm. Particles are introduced at the inlet and create an erodible bed thanks to the obstacle placed at the outlet. The number of particle layers is controlled by the height of this obstacle. The channel bottom is made of metal half-cylinders of diameter d, fixed at a random elevation between -2.75 mm and 2.75 mm with steps of 0.5 mm to break

FIG. 2 .

 2 FIG. 2. Experimental setup scheme, modified from Böhm et al 47 and Hergault et al 48 . The inclined channel width of 6.5d/6 implies a quasi-2D bead flow, permitting particle tracking of each spherical bead in the observation window filmed by the camera placed perpendicularly.

FIG. 3 .

 3 FIG. 3. Depth profiles of average solid velocity (m/s), volume fraction and sediment transport rate density (m/s) for the case Sim20 with a periodic length cell of 1000 d. The different black lines (-) correspond to different post-processing averaging performed in the experimental condition, i.e. over boxes of streamwise length 40 d and time-averaged over 60 s. The full red dots (•) correspond to the averaging as performed for the simulation in general, with a period of averaging of 100 s and a streamwise length of the size of the periodic cell. The figure shows the order of magnitude of the variability of the experimental results due to the limited spatio-temporal window of averaging.

  [START_REF] Frey | Particle velocity and concentration profiles in bedload experiments on a steep slope[END_REF] , θ = ρ f Rh b S 0 /[(ρ p -ρ f )d], with Rh b the corrected water depth; and from the fluid bottom shear stress defined by the friction velocity u * : θ

FIG. 4 .FIG. 5 .

 45 FIG. 4. Experimental comparison for the different cases presented in table III: (a) Case 20, (b) Case 14, (c) Case 6. The figure shows for each case the depth profiles of the average streamwise solid velocity (m/s), solid volume fraction, and sediment transport rate density (m/s). The full symbol (•, , ) represents the experimental results from Frey 27 while the empty linked one represents the simulation (+, x,•). The black line represents the imposed free surface position. The error bars show the variability of the experimental results as evaluated from figure 3. The results show a good general agreement for the three profiles with well-reproduced trends in each case.

  further in the absence of experimental data. Considering the comparison for the three different cases, with respect to the simplicity of the fluid description and the goal of describing the average solid behavior, the agreement with the experiments is good. The values of the integrated transport rate are close to the experimental ones and the sensitivity to the experimental parameters such as the free-surface position or the number of bead layers has been well reproduced. The comparison of the averaged depth profile of the solid velocity, the solid volume fraction and the transport rate

25 FIG. 6 .

 256 FIG. 6. Dimensionless sediment transport rate Q * s as a function of Shields number θ for different configuration. Classical runs ( ) with e n = 0.5 and µ = 0.4 are shown together with the exact same runs without turbulent fluctuation model ( ) for different Shields number. Triangle symbols represent the change in restitution coefficient with e n = 0.01 ( ), e n = 0.25 ( ), e n = 0.75 ( ), e n = 1 ( ). Variations of the particle friction coefficient are represented by bullet points: µ = 0.2 (•),µ = 0.6 (•),µ = 0.8 (•). The experimental data of Meyer-Peter and Müller 2 (+) and Wilson 55 (x) synthesized in Yalin 56 , show the experimental trend in power law 3/2, with the dispersion of the data. The gray line corresponds to Q * s = 11.8(θ * -θ * c ) 3/2 as found asymptotically by Wilson 57 .The inset in linear scale shows the behavior near the threshold of motion.

  for two different Shields numbers, by keeping the free-surface position and number of particle layers constant. The restitution coefficient has been varied in the range e n ∈ [0.01, 1]. It corresponds to a realistic range e n ∈ [0.25, 0.75], complemented by two extreme cases: no rebounds (e n = 0.01)

FIG. 7 .

 7 FIG. 7. Effect of the restitution (a) and friction coefficient (b) on the average solid depth profiles for a Shields number θ ∼ 0.45. The solid velocity v p x and sediment transport rate density q s are given in m/s, while the solid volume fraction φ is dimensionless. To give a global picture of the trend, the color of the lines are representative of the friction and restitution coefficient values. The position of the free surface in both case is situated at 20d. The curve are plotted only for values of solid volume fraction larger than 10 -3 .

  the model is able to reproduce the classical macroscopic validation in term of sediment transport rate and Shields number. In addition, an original detailed validation with existing bedload transport experiments has been performed, comparing simultaneous measurements of average solid velocity and volume fraction. The good agreement with experiments together with the rather low sensitivity of the results to the granular parameters show the relevancy and the robustness of the proposed model, which reproduces not only the evolution of the sediment transport rate as a function of the Shields number, but also the depth structure of the granular phase. The influence of the different model contributions have been studied. In particular, the discrete random walk fluid velocity fluctuations model has been shown to be sufficient to reproduce the reduction of the critical Shields number due to turbulent fluctuations. A weak impact of the restitution and friction coefficients variations has been observed on the macroscopic sediment transport rate versus Shields number curve. The analysis of the depth profiles variations shows however that the granular parameters influence the depth structure of the granular flows and induce non-trivial coupling with the fluid phase. The rigorous development of the model and the experimental validations demonstrate the potential of this modeling approach to deal with granular processes in bedload transport. Future work will take advantage of the description of the depth structure to analyze the effective granular rheology under bedload transport conditions.

  A2) where the RMS ∆Q rms i is normalized by the average transport rate of the reference configuration Q ref , N is the number of averaging cell in the depth, Q i z and Q ref z are the values of transport rate in the cell z for the considered case and the reference case respectively.

Figure 8 1 2xFIG. 8 . 5 x(

 8185 Figure 8 shows the normalized deviation with respect to the reference configuration defined by equation (A2) as a function of the streamwise periodic cell length of the simulation considered. The time-averaged results show a convergence as a function of cell length l x of the order of l -1 2x . The size of the periodic cell used for the simulations presented in the

3D

  convergence analysis A similar analysis has been undertaken for the three dimensional bi-periodic configuration at a Shields number around 0.1. The reference case has been chosen as l x = l y = 100d, i.e. the relative cell size V = l x l y l z being the same as the quasi-2D case. Different simulations have been performed with a cell of l x = l y = 10 d, 20 d, 30 d,50 d, and three cases with l x = l y : (l x , l y ) = (500 d, 5 d), (300 d, 5 d), and (10 d, 100 d) . It has not been possible to consider smaller cell sizes, as the coupled model becomes unstable.

1 FIG. 9 .

 19 FIG. 9. Convergence of the sediment transport rate profile as a function of the periodic cell size considered for the 3D case. The vertical axis represents the deviation with respect to the reference configuration (l x = l y = 100d) as defined in equation (A2). Cases with l x = l y are represented with •, while the triangles denotes cases with l x = l y : (l x , l y ) = (500 d, 5 d) ( ), (300 d, 5 d) ( ), and (10 d, 100 d) ( ). The logarithmic scale inset shows that the results are already converged.

TABLE I .

 I Characteristic values of the main dimensionless numbers. The particle feeding rate is controlled, and the flow rate is adjusted in order to reach transport equilibrium, i.e. feeding rate equal to the sediment transport rate at the outlet without having aggradation and degradation of the bed. The free-surface fluid flow is turbulent (Re = U f w d /ν f ∼ 10 4 ), hydraulically rough (Re p ∼ 10 3 ), and supercritical

	θ	Re	Re p	F r	ρ p /ρ f	St	S *
	0.05-0.1	10 4	10 3	1	2.5	10 2 -10 3	2 -10
	clusterization.

TABLE II .

 II Model input parameters for the contact law and the fluid resolution. k n and k s are respectively the normal and tangential contact stiffness, e n and µ denote the restitution and friction coefficient, κ is the Von Karman constant, ζ the Richardson-Zaki exponent, φ max the bed solid volume fraction, and τ f the fluid resolution period.

	k n (N/m)	k s (N/m)	e n	µ	κ	ζ	φ max	τ f (s)
	5 10 3	2.5 10 3	0.5	0.4	0.41	3.1	0.51 (2D)/0.61 (3D)	10 -2

TABLE III .

 III Experimental and numerical run characteristics. The free surface position h and the number of bead layers N l are both measured in the experiment and imposed in the simulation. ṅ is the measured transport rate, θ and θ * the Shields numbers respectively based on macroscopic flow parameters and turbulent shear stress profile. The latter has only been determined in the

	simulations.
	Run	h (cm)	N l	ṅ (b/s)	θ	θ *
	Exp6	5.3	7.08	6.67	0.076	-
	Sim6	5.3	7.08	10.15	0.083	0.031
	Exp14	5.7	7.37	13.68	0.100	-
	Sim14	5.7	7.37	18.13	0.120	0.048
	Exp20	5.9	7.30	19.74	0.106	-
	Sim20	5.9	7.30	26.38	0.130	0.061

TABLE IV .

 IV Results of the sensitivity analysis to the fluid resolution period τ f , for the case Sim20.The measured Shields number θ * , transport rate ṅ, and RMS deviation with respect to the case Ref. ∆Q rms i Q ref , are given for each case.

	Case	τ f (s)	θ *	ṅ (b/s)	∆Q rms Q ref	i
	Ref.	10 -2	0.067	16.83	0.0
	τ f = 10 -3	10 -3	0.067	16.74	0.06
	τ f = 10 -1	10 -1	0.067	16.91	0.04
	τ f = 1	1	0.066	17.26	0.13

ACKNOWLEDGMENTS

We are grateful to Michael Church for general remarks and English corrections. This research was supported by Irstea (formerly Cemagref) and the French Institut National des Sciences de l'Univers program EC2CO-BIOHEFECT, and EC2CO-LEFE MODSED .

Appendix A: Averaging expression and convergence analysis

The averaging of the solid phase takes a central part in the coupling between the Lagrangian solid phase and the Eulerian fluid phase. The important wall-normal gradient requires the length scale of the weighting function in this direction to be lower than particle diameter (l z d/30 for the lowest Shields number) in order to define a rigorous averaging.

We postulated that the complementary length scales l x and l y can statistically compensate the limited l z .

In the present model, the average fluid description is 1D so that it depends only on the wall-normal component, z. The solid averaging can therefore be performed on the full width and length cell. With the cuboidal formulation of the weighting function defined previously (eq. 8), the solid averaging of a scalar particle quantity γ at a wall-normal position z can be rewritten as:

Where l z is the defined wall-normal weighting function length scale, and Ṽ p is the fraction of the particle volume contained in the slice of height l z at elevation z. We recover here the averaging formulation of Hill et al [START_REF] Hill | Structure and kinematics in dense free-surface granular flow[END_REF] , which is convenient to compute since the volume of a slice of spheres can be evaluated analytically. The averaging box height l z is imposed by the vertical grid size of the fluid problem and no overlapping between the different slices is allowed.

For each fluid resolution, so at each given time step, the statistical representativity of the averages is a requirement for a consistent definition of the averaging process (section II C 1). The spatial convergence of the averages with increasing complementary length scales l x and l y includes the effect of the bottom boundary conditions and particles arrangement, in addition to the statistical representativity. The results are required to be independent from the three effects, and consequently the spatial convergence of the results with respect to l x and l y is analyzed in the present appendix.

There are two different convergence scale in the problem. The first one is associated to the spatial convergence at each given time step, and the second one is associated to the temporal convergence of a simulation with a given cell size. In the present analysis, as the paper

Appendix B: fluid resolution period

The DEM time step needs to be particularly low and the evolution of the granular medium over this time is limited. Consequently the fluid resolution period τ D does not need necessarily to be equal to the solid time step. The stability of the coupling however depends on this period of resolution and it is important to study this parameter in order to have a meaningful model. The fluid resolution period should be defined to be smaller than the characteristic time of evolution of the granular medium. The fluid resolution is 1D and the equation is influenced only by the streamwise particle velocity and the wall-normal particle position (through respectively f p D,x and φ). For a single particle the evolution of these properties depends on the collisions and the entrainment by the fluid. As seen previously (section IV), collisions do not significantly influence the phenomenon so that we consider only the characteristic time of entrainment. As explained in section II C, the characteristic time of relaxation to accelerate a particle to the fluid velocity is given by τ D = β -1 , with β expressed from equation (18). In the present case it is of order τ D ∼ 10 -1 s. However, the characteristic time associated with each independent particle is not in general representative of the evolution of a complex many-body problem.

We therefore performed a sensitivity analysis on the period of fluid resolution. The results are shown in table IV in term of RMS deviation (as defined by equation (A2)) with respect to the reference configuration Sim20 for which τ f = 10 -2 s. It includes τ f = 10 -3 s, 10 -1 s, and 1s. The results exhibit no dependence on the fluid period resolution in the range 10 -3 s to 1s. The values of the RMS deviation with respect to the reference configuration τ f = 10 -2 s is for all cases below 0.1, i.e. below the reproducibility deviation value. The fluid resolution period therefore does not have an influence on the averaged equilibrium results within the range considered. This means that the solid average quantities transmitted to the fluid do not vary importantly during the simulation after reaching equilibrium. These results confirm that the simulations considered are at transport equilibrium. The fluid resolution period could however be important for unsteady conditions.