Raphael Maurin

Julien Chauchat
email: julien.chauchat@grenoble-inp.fr

Philippe Frey

Dense granular flow rheology in turbulent bedload transport

Keywords: Sediment transport, granular media, rheology

The local granular rheology is investigated numerically in turbulent bedload transport. Considering spherical particles, steady uniform configurations are simulated using a coupled fluid-discrete-element model. The stress tensor is computed as a function of the depth for a series of simulations varying the Shields number, the specific density and the particle diameter. The results are analyzed in the framework of the µ(I) rheology and exhibit a collapse of both the shear to normal stress ratio and the solid volume fraction over a wide range of inertial numbers. Contrary to expectations, the effect of the interstitial fluid on the granular rheology is shown to be negligible, supporting recent work suggesting the absence of a clear transition between the free-fall and turbulent regime. In addition, data collapse is observed up to unexpectedly high inertial numbers I ∼ 2, challenging the existing conceptions and parametrization of the µ(I) rheology. Focusing upon bedload transport modelling, the results are pragmatically analyzed in the µ(I) framework in order to propose a granular rheology for bedload transport. The proposed rheology is tested using a 1D volume-averaged two-phase continuous model, and is shown to accurately reproduce the dense granular flow profiles and the sediment transport rate over a wide range of Shields numbers. The present contribution represents a step in the upscaling process from particle-scale simulations toward large-scale applications involving complex flow geometry.

Introduction

Among the different regimes of sediment transport, bedload transport is of major importance as it represents an important contribution to river morphology evolution. Accordingly, prediction of the sediment transport rate under turbulent bedload conditions is fundamental for preventing environmental risks associated with floods and scouring. From a physical point of view, bedload transport corresponds to the dynamic response of a granular bed submitted to a fluid shear stress. In contrast to suspended load, bedload is defined as the part of the sediment load occuring close to the granular bed in which particles are in permanent or intermittent contact with the bed [START_REF] Fredsøe | Mechanics of Coastal Sediment Transport[END_REF].

The sediment transport rate is usually quantified in dimensionless form using the so-called Einstein number [START_REF] Einstein | Formulas for the transport of bed sediment[END_REF]

, Q * s = Q s /[d (ρ p /ρ f -1)gd],
where Q s represents the volumetric sediment transport rate per unit width, ρ p and ρ f are the particle and fluid densities, g is the acceleration of gravity and d is the particle diameter. In the hydraulic community the classical way to model bedload transport consists of assuming a power law relation between the Einstein number and the Shields number θ. The latter corresponds to the ratio of the traction force exerted by the fluid at the bed τ f d 2 and the buoyant weight of a single particle (ρ p -ρ f )gd 3 , i.e. θ = τ f /[(ρ p -ρ f )gd]. Given the lack of accuracy of the classical formulae (e.g. [START_REF] Meyer-Peter | Formulas for bed-load transport[END_REF]), renewed attention has recently been given to the granular phase behaviour in bedload transport [START_REF] Frey | How river beds move[END_REF], 2011) from both an experimental [START_REF] Mouilleron | Inside the moving layer of a sheared granular bed[END_REF][START_REF] Hergault | Image processing for the study of bedload transport of two-size spherical particles in a supercritical flow[END_REF][START_REF] Lajeunesse | Bed load transport in turbulent flow at the grain scale: Experiments and modeling[END_REF][START_REF] Frey | Particle velocity and concentration profiles in bedload experiments on a steep slope[END_REF][START_REF] Aussillous | Investigation of the mobile granular layer in bedload transport by laminar shearing flows[END_REF][START_REF] Houssais | Onset of sediment transport is a continuous transition driven by fluid shear and granular creep[END_REF] and a numerical point of view [START_REF] Duran | Numerical simulation of turbulent sediment transport, from bed load to saltation[END_REF][START_REF] Ji | Direct numerical simulation of sediment entrainment in turbulent channel flow[END_REF][START_REF] Kidanemariam | Interface-resolved direct numerical simulation of the erosion of a sediment bed sheared by laminar channel flow[END_REF][START_REF] Maurin | A minimal coupled fluid-discrete element model for bedload transport[END_REF].

The granular rheology characterises the response of the granular medium, in terms of deformation rate, to a given external stress and vice-versa. It is usually given as the relation between the granular stress and strain rate tensors. In bedload transport, the granular rheology governs both the response to, and the interaction with the fluid shear stress. Its knowledge is required for two-phase flow models in which both the fluid and solid phases are considered as continuous (Eulerian-Eulerian) (e.g. [START_REF] Aussillous | Investigation of the mobile granular layer in bedload transport by laminar shearing flows[END_REF]). As such it represents an important step in the upscaling of sediment transport processes from Eulerian-Lagrangian models to Eulerian-Eulerian models, which enable numerical simulations of larger-scale problems and/or complex flow forcings.

Under bedload transport conditions, the granular medium experiences all of the granular flow regimes, from quasi-static in the sediment bed to dilute rapid granular flows in the upper sediment transport layer. The dilute rapid and dense granular flow regimes have mainly been described using the kinetic theory of granular flows and the µ(I) rheology respectively. The former is based on the analogy of dilute granular flows with molecular gases, assuming binary collisions. It has been proven to accurately describe rapid dilute granular flows in different configurations [START_REF] Campbell | Rapid granular flows[END_REF][START_REF] Goldhirsch | Rapid granular flows[END_REF], and has recently been extended to the dense granular flow regime [START_REF] Jenkins | Dense shearing flows of inelastic disks[END_REF][START_REF] Jenkins | Dense inclined flows of inelastic spheres[END_REF][START_REF] Berzi | Constitutive relations for steady, dense granular flows[END_REF]. Alternatively, the µ(I) local rheology accurately describes the dense granular flow regimes and is based on the dimensional analysis of the simple shear configuration (Midi 2004;[START_REF] Forterre | Flows of dense granular media[END_REF][START_REF] Jop | Rheological properties of dense granular flows[END_REF]. In the latter case, the unique dimensionless number controlling the system is the so-called inertial number [START_REF] Da Cruz | Rheophysics of dense granular materials: Discrete simulation of plane shear flows[END_REF]:

I = γd ρ p /P p ,
(1.1) with P p the confining granular pressure and γ the granular shear rate. This dimensionless number can be interpreted as the ratio between a macroscopic time scale of deformation t macro = 1/ γ and a microscopic time scale of rearrangement t micro = d/ P p /ρ p . From the dimensional analysis, the dimensionless solid volume fraction φ and shear to normal granular stress ratio µ = τ p /P p are unique functions of the inertial number I. From experiments and numerical simulations, they can be expressed as [START_REF] Da Cruz | Rheophysics of dense granular materials: Discrete simulation of plane shear flows[END_REF][START_REF] Jop | A constitutive law for dense granular flows[END_REF])

τ p P p = µ(I) = µ s + µ 2 -µ s 1 + I 0 /I , (1.2) φ(I) = φ max -aI, (1.3)
where φ max is the maximum packing fraction, µ s is the static effective granular friction coefficient, and µ 2 , I 0 and a are phenomenological constants. By combining the different expressions (eq. 1.1-1.3), the granular stress tensor can be expressed as a function of the shear rate, defining the granular rheology.

Based on the simple shear configuration, the µ(I) rheology has been extended to account for the presence of interstitial fluids [START_REF] Courrech Du Pont | Granular avalanches in fluids[END_REF][START_REF] Cassar | Submarine granular flows down inclined planes[END_REF]. Considering rearrangement time scales according to the dominant mechanism, three regimes can be defined by introducing two dimensionless numbers [START_REF] Courrech Du Pont | Granular avalanches in fluids[END_REF][START_REF] Cassar | Submarine granular flows down inclined planes[END_REF][START_REF] Andreotti | Granular media: between fluid and solid[END_REF], with C D the drag coefficient,

St = d √ ρ p P p η f , (1.4) r = ρ p ρ f C D , (1.5)
the free-fall regime (St >> 1, r >> 1) corresponding to negligible influence of the interstitial fluid (i.e. dry granular media), the viscous regime (St << 1, r << 1) corresponding to a rearrangement time scale dominated by viscous drag and the turbulent regime (St >> 1, r << 1) corresponding to a rearrangement time scale dominated by fluid inertial effects. In the different regimes, the results therefore scale with the dry inertial number I dry = γd ρ p /P p , the turbulent one I turb = γd ρ f C D /P p and the viscous one I visc = η f γ/P p respectively. This approach has been applied with some success to various complex immersed configurations from avalanches (Courrech du [START_REF] Courrech Du Pont | Granular avalanches in fluids[END_REF][START_REF] Cassar | Submarine granular flows down inclined planes[END_REF][START_REF] Doppler | Relaxation dynamics of water-immersed granular avalanches[END_REF], to granular collapses [START_REF] Rondon | Granular collapse in a fluid: Role of the initial volume fraction[END_REF][START_REF] Izard | Simulation of an avalanche in a fluid with a soft-sphere/immersed boundary method including a lubrication force[END_REF], annular shear cells [START_REF] Boyer | Unifying suspension and granular rheology[END_REF][START_REF] Trulsson | Transition from the viscous to inertial regime in dense suspensions[END_REF]) and sediment transport [START_REF] Ouriemi | Sediment dynamics. part 1. bed-load transport by laminar shearing flows[END_REF][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][START_REF] Chiodi | A two-phase flow model of sediment transport: transition from bedload to suspended load[END_REF].

To the best of our knowledge, there are only a few contributions in the literature on the dense granular flow rheology in the context of bedload transport. [START_REF] Ouriemi | Sediment dynamics. part 1. bed-load transport by laminar shearing flows[END_REF] and [START_REF] Aussillous | Investigation of the mobile granular layer in bedload transport by laminar shearing flows[END_REF] have studied laminar bedload in closed conducts using refractive index-matching experiments and theoretical two-phase continuous models. Assuming a constant solid volume fraction and a limited transitional layer from dense to dilute granular flow, the authors have shown that the µ(I) rheology gives excellent agreement with experimental data. It is worth noting that in the laminar case, the lag between the fluid and solid velocities is negligible [START_REF] Mouilleron | Inside the moving layer of a sheared granular bed[END_REF][START_REF] Aussillous | Investigation of the mobile granular layer in bedload transport by laminar shearing flows[END_REF], the two phases being tightly coupled. In intense turbulent bedload transport, termed sheet flow, the transitional layer from dense to dilute granular flows and the fluid-particle velocity lag are more important [START_REF] Sumer | Velocity and concentration profiles in sheet-flow layer of movable bed[END_REF][START_REF] Cowen | An insitu borescopic quantitative imaging profiler for the measurement of high concentration sediment velocity[END_REF][START_REF] Revil-Baudard | Investigation of sheet-flow processes based on novel acoustic high-resolution velocity and concentration measurements[END_REF]. This regime has been studied experimentally [START_REF] Capart | Transport layer structure in intense bed-load[END_REF] and numerically using the two-phase flow continuous modelling framework with the kinetic theory of granular flows [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]) and 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] Chauchat | An improved two-phase flow model for steady uniform sheet-flow based on dense granular flow rheology[END_REF]. While the agreement of the numerical simulations with the available experimental data gives some credit to both granular rheologies, recent experimental investigations by [START_REF] Revil-Baudard | Investigation of sheet-flow processes based on novel acoustic high-resolution velocity and concentration measurements[END_REF] showed that the solid volume fraction profiles observed were substantially different from the predicted ones. The observed effective friction coefficient also differs importantly from the classical values used in the µ(I) rheology, possibly due to the strong intermittency induced by the turbulent coherent structures [START_REF] Revil-Baudard | Investigation of sheet-flow processes based on novel acoustic high-resolution velocity and concentration measurements[END_REF]. At the transition from sheet flows to debris flow, the granular rheology has been used to model steep slope configurations [START_REF] Armanini | Rheological stratification in experimental free-surface flows of granular-liquid mixtures[END_REF][START_REF] Larcher | Set of measurement data from flume experiments on steady uniform debris flows[END_REF][START_REF] Berzi | A theoretical analysis of free-surface flows of saturated granular-liquid mixtures[END_REF][START_REF] Armanini | Submerged granular channel flows driven by gravity[END_REF]. Combining experimental particle tracking at the wall [START_REF] Armanini | Rheological stratification in experimental free-surface flows of granular-liquid mixtures[END_REF][START_REF] Larcher | Set of measurement data from flume experiments on steady uniform debris flows[END_REF]) and numerical analysis [START_REF] Berzi | A theoretical analysis of free-surface flows of saturated granular-liquid mixtures[END_REF][START_REF] Armanini | Submerged granular channel flows driven by gravity[END_REF], both the µ(I) rheology [START_REF] Berzi | A theoretical analysis of free-surface flows of saturated granular-liquid mixtures[END_REF]) and hybrid µ(I)-kinetic theory models [START_REF] Armanini | Submerged granular channel flows driven by gravity[END_REF] have been shown to accurately describe the granular behaviour in this configuration.

The literature review underlines the relevance of both the kinetic theory and the dense granular flow rheology for bedload transport modelling. However, the numerical studies are restricted to continuous two-phase analysis and the experimental evaluations of the granular stress tensor rely on solid volume fraction measurements that are highly uncertain. In addition, while direct numerical simulations (DNS) have been performed at the particle scale in the laminar regime [START_REF] Kidanemariam | Interface-resolved direct numerical simulation of the erosion of a sediment bed sheared by laminar channel flow[END_REF], this type of numerical model is not affordable for turbulent bedload transport due to the typical values of bulk Reynolds numbers explored. Therefore, to go beyond the existing twophase continuous numerical works, the present paper analyses the local granular rheology in turbulent bedload transport using a coupled fluid-discrete-element model (DEM) with a volume-averaged fluid description [START_REF] Maurin | A minimal coupled fluid-discrete element model for bedload transport[END_REF]. This approach permits local computation of the granular stress tensor and the granular shear rate, bringing new insights into the local granular behaviour in turbulent bedload transport. Based on these Eulerian-Lagrangian simulations a dense granular flow rheology for bedload transport is proposed and further tested in a two-phase continuous (Eulerian-Eulerian) model. Besides, bedload transport configurations enable local analysis of the granular rheology in complex immersed granular flows, and the DEM results are shown to be of interest for granular media in general.

The paper is structured as follows. First, the coupled fluid-DEM and two-phase continuous models are briefly presented (section 2). Then, fluid-DEM simulations are analyzed in the framework of the dense granular flow µ(I) rheology (section 3). Finally, extending the simulations to realistic conditions for bedload transport, a parametrization of the µ(I) rheology is proposed from the DEM results, and is tested with the two-phase continuous model (section 4).

Model formulation

The Eulerian-Lagrangian and Eulerian-Eulerian models are based on the same volumeaveraged two-phase flow equations for the fluid phase. The solid phase is modelled on the one hand as a continuum (Eulerian-Eulerian) and on the other hand using a 3D DEM in which each particle motion is computed explicitly (Eulerian-Lagrangian). Both two-phase flow models have already been described and validated with experiments [START_REF] Maurin | A minimal coupled fluid-discrete element model for bedload transport[END_REF][START_REF] Chauchat | An improved two-phase flow model for steady uniform sheet-flow based on dense granular flow rheology[END_REF] and only a brief description will be given. The reader is referred to the abovementioned references for a complete description of the model formulations and validations.

Considering steady uniform conditions, the problem is unidirectional so that the average fluid and solid velocities depend only on the wall-normal direction and reduce to their streamwise components u f = u x f (z) e x and v p s = v p x s (z) e x . Therefore, the streamwise and wall-normal volume-averaged fluid phase momentum balances are given by [START_REF] Anderson | Fluid mechanical description of fluidized beds. equations of motion[END_REF][START_REF] Jackson | The dynamics of fluidized particles[END_REF][START_REF] Revil-Baudard | A two-phase model for sheet flow regime based on dense granular flow rheology[END_REF]

0 = ∂S f xz ∂z + ∂R f xz ∂z + ρ f (1 -φ)g sin α -n f x p , (2.1) 0 = - ∂P f ∂z + ρ f (1 -φ)g cos α -n f z p , (2.2)
where σ f ij = -P f δ ij + S f ij is the volume-averaged effective viscous stress tensor, R f xz is the Reynolds shear stress tensor, φ is the solid phase volume fraction, f k p is the volume-averaged fluid-particle interaction force, n = φ/(πd 3 /6) is the number density of particles and α is the channel inclination angle. Omitting the model-dependent fluid-particle interaction force, the solution of the fluid momentum balance requires closure laws for the viscous shear stress tensor and the Reynolds stress tensor.

Considering a Newtonian fluid, the viscous shear stress is classically expressed as

S f xz = ρ f (1 -φ)ν f d u x f dz , (2.3)
with ν f the clear fluid kinematic viscosity and u x f the volume-averaged streamwise fluid phase velocity. The Reynolds shear stress is based on the eddy viscosity concept (ν t ) using a mixing length formulation:

R f xz = ρ f ν t d u x f dz with ν t = (1 -φ) l 2 m d u x f dz , (2.4)
where the mixing length is taken similarly to [START_REF] Li | Multi-phase model on sediment transport in sheet-flow regime under oscillatory flow[END_REF] as

l m (z) = κ z 0 φ max -φ(ζ) φ max dζ, (2.5)
with κ = 0.41 the von Karman constant. The formulation adopted allows one to recover the law of the wall [START_REF] Prandtl | Bericht über neuere Turbulenzforschung[END_REF] in clear fluid, while the turbulence is completely damped inside the granular bed at maximum packing fraction (φ max ).

Eulerian-Lagrangian model

The Eulerian-Lagrangian model is based on the explicit solution of the dynamic equation for each individual particle using the 3D DEM open-source code YADE ( Šmilauer et al. 2015). The DEM solution is spatially averaged and explicitly coupled with the fluid phase momentum balance equations (2.1-2.2) through the drag term n < f D

x > p and the solid volume fraction φ.

For each particle p at position x p , the Newton equations are solved considering nearestneighbours interactions:

m d 2 x p dt 2 = k∈N f pk c + f ext = k∈N f pk c + f p g + f p b + f p D (2.6)
where the sum of the contact forces f pk c is made over the ensemble of nearest neighbours N , f p g is the gravity force, f p b is the buoyancy force and f p D is the 3D drag force applied by the fluid on particle p. For each contact, the contact force is computed explicitly using the classical spring-dashpot contact law [START_REF] Schwager | Coefficient of restitution and linear spring-dashpot model revisited[END_REF]:

F n = -k n δ n -c n δ n (2.7a) F t = -min(k s δ t , µ p F n ), (2.7b)
where F n and F t are the normal and tangential contact forces between particle p and k, δ n and δ t are the normal and tangential overlaps, k n and k s are the normal and tangential contact stiffnesses, c n is the normal viscous damping and µ p is the tangential friction coefficient. This contact law is well suited for granular flow analysis and allows one to define a unique restitution coefficient e n [START_REF] Schwager | Coefficient of restitution and linear spring-dashpot model revisited[END_REF]. Similarly, the rotation of the particles is solved from the Newton equations of motion.

The interaction with the fluid phase is restricted to buoyancy and drag forces:

f p b = - πd 3 6 ∇P f , (2.8) f p D = 1 2 ρ f πd 2 4 C D u f x p -v p u f x p -v p , (2.9)
where the average fluid velocity and the fluid pressure are taken at the center of particle p. The drag coefficient C D depends on the particle Reynolds number

Re p = u f x p -v p d/ν f and takes into account hindrance effects (DallaValle 1948; Richardson & Zaki 1954): C D = (0.4 + 24.4/Re p ) (1 -φ) -3.1 .
Considering steady uniform turbulent flows, the buoyancy force reduces to its wallnormal component and equation (2.2) leads to a hydrostatic fluid pressure distribution, the wall-normal average drag force being negligible. The solution of the streamwise fluid phase momentum balance (eq. 2.1) requires the evaluation of both the spatially averaged solid phase volume fraction and the momentum transmitted from the fluid to the particles through drag forces:

n < f D x > s = φ πd 3 /6 < f D x > s = 3 4 φ ρ f d C D u f -v p u x f -v p x s
.

(2.10)

In the DEM processing, the solid phase averaging < • > s is used instead of the particle phase averaging < • > p . While the two are identical provided that there is scale separation [START_REF] Anderson | Fluid mechanical description of fluidized beds. equations of motion[END_REF][START_REF] Jackson | The dynamics of fluidized particles[END_REF], the importance of the wall-normal gradients in bedload transport requires a weighting function length scale lower than the particle diameter in order to define an independent averaging [START_REF] Maurin | A minimal coupled fluid-discrete element model for bedload transport[END_REF][START_REF] Maurin | Investigation of granular behavior in bedload transport using an Eulerian-Lagrangian model[END_REF]. Therefore, a small wall-normal weighting function length scale has been adopted (typically d/30), and this choice has been validated through an experimental comparison [START_REF] Maurin | A minimal coupled fluid-discrete element model for bedload transport[END_REF].

The 3D DEM and the fluid model are solved as transient problems applying a fixed bottom boundary condition for both the fluid ( u f (z = 0) = 0) and the particle phase (fixed random particles) and imposing the position of the water free-surface (d u x f /dz(z = h) = 0). In order to achieve a stable integration, the DEM time step is bounded by the propagation time of the fastest wave over a particle diameter [START_REF] Maurin | Investigation of granular behavior in bedload transport using an Eulerian-Lagrangian model[END_REF][START_REF] Maurin | A minimal coupled fluid-discrete element model for bedload transport[END_REF]. The fluid resolution time step corresponds to a typical characteristic evolution time scale of the granular medium and is taken much larger than the DEM one [START_REF] Maurin | A minimal coupled fluid-discrete element model for bedload transport[END_REF]: ∆t f = 10 -2 s with respect to ∆t p ∼ O(10 -4 -10 -5 ) s. Therefore, the coupling between the fluid and granular phases ensures momentum conservation of the system on average. The model has been compared with experiments and has shown its ability to describe accurately the granular depth structure in turbulent bedload transport [START_REF] Maurin | A minimal coupled fluid-discrete element model for bedload transport[END_REF].

Eulerian-Eulerian model

The Eulerian-Eulerian model is based on the numerical solution of the fluid phase momentum balance equations (eqs 2.1 and 2.2) coupled with the spatially averaged granular phase momentum balance equations in the streamwise and wall-normal directions. For steady uniform flow conditions the granular phase momentum equations read [START_REF] Chauchat | An improved two-phase flow model for steady uniform sheet-flow based on dense granular flow rheology[END_REF]

0 = ∂τ p xz ∂z + ρ p φg sin α + n f x p , (2.11) 0 = - ∂P p ∂z + ρ p φg cos α + n f z p , (2.12)
where τ p xz is the spatially-averaged granular phase shear stress. Like in the Eulerian-Lagrangian model, the buoyancy is applied along the wall-normal direction and leads to hydrostatic pressure distribution for the granular phase. In the continuous formalism the spatially averaged streamwise drag force is given by

n f D s = 3 4 φρ f d C D u f -v p s + u d u f -v p s + u d , (2.13)
where the extra term u d is the drift velocity which represents the dispersion effect due to the averaged fluid-particle velocity fluctuations. This term is responsible for sediment transport suspension and allows one to reproduce the Rouse profile provided that the wall normal component is taken as: (u d ) z = -ν t /φ dφ/dz [START_REF] Chauchat | An improved two-phase flow model for steady uniform sheet-flow based on dense granular flow rheology[END_REF]. The streamwise component is taken as (u d ) x = 0 in this paper.

The granular shear stress tensor is modeled using the µ(I) rheology presented in the introduction. Combining equations (1.1) and (1.2), the shear and normal components of the granular stress tensor can be expressed as

P p = 1 I 2 ρ p d 2 d v p x s dz 2 , (2.14) τ p = µ(I) I 2 ρ p d 2 d v p x s dz 2 .
(2.15)

Inverting equation 1.3, the inertial number is obtained as a function of the solid volume fraction and the two equations can be directly expressed as a function of the solid volume fraction and the average solid velocity gradient.

The numerical algorithm is based on a transient solution of the governing equations (2.1; 2.2; 2.11 and 2.12), together with the continuity equation for the solid phase. Fixed boundary conditions are imposed at the bottom of the sample at z = 0 ( u f = 0) where the granular phase is considered to be at rest ( v p s = 0, dφ/dz = 0), while the positions of the water and granular free-surfaces are imposed in z = h (d u x f /dz = 0, u z f = 0, v p z s = 0, dφ/dz = 0). The governing equations are discretised using a finite volume technique for the mass conservation equation and a finite difference technique for the momentum balance. A staggered grid is used with the velocities located at the cell face and the scalar quantities (e.g. volume fractions or viscosities) located at the cell center. x (z) ex, solid volume fraction φ, and solid velocity v s = vx s (z) ex depth profiles. The inclined 3D bi-periodic granular description is coupled with a unidirectional fluid momentum balance using imposed fixed random bottom boundary condition and water free-surface position. The particle color is representative of the velocity intensity.

θ * Re Rep ρ p /ρ f -1 F r S * St r [0.04, 0.6] 10 3 -10 6 10 3 -10 4 [0.375, 1.5] 1 [2 , 5] 10 2 -10 4 0.1 -2
Table 1. Main characteristic dimensionless numbers of the configurations considered.

For the pressure-velocity coupling a projection method is used. The numerical schemes are Euler implicit for the time derivative, upwind for the advection terms and central difference for the diffusion terms. For additional informations on the numerical model, the reader is referred to [START_REF] Chauchat | Modelling sedimentation-consolidation in the framework of a one-dimensional two-phase flow model[END_REF] and [START_REF] Chauchat | An improved two-phase flow model for steady uniform sheet-flow based on dense granular flow rheology[END_REF].

Results and methodology

The numerical setup of the Eulerian-Lagrangian model is presented in figure 1. It consists of a channel flow tilted with an inclination angle α = 0.05 rad, partially filled with monodisperse spherical particles of diameter d. The three-dimensional granular sample is biperiodic (streamwise and spanwise) and the periodic cell size is taken as l x = l y = 30d to ensure statistical convergence of the spatial averaging operator [START_REF] Maurin | A minimal coupled fluid-discrete element model for bedload transport[END_REF]. A 10-15 diameter thick granular layer is deposited under gravity over a rough random bottom, and the system evolves under gravity with a fixed water freesurface elevation. The water flow is turbulent (Re = U f h/ν f ∼ 10 4 ), hydraulically rough (Re p = U f d/ν f ∼ 10 3 ) and supercritical (F r = U f / √ gh 1), with U f the average fluid velocity within the water depth h. The results are independent of the granular bottom boundary conditions in the range of parameters investigated [START_REF] Maurin | Investigation of granular behavior in bedload transport using an Eulerian-Lagrangian model[END_REF] and the simulations are performed in the rigid grain limit [START_REF] Roux | Quasistatic rheology and the origins of strain[END_REF] with a tangential stiffness set to half the normal one. The friction coefficient is set to a realistic value for glass beads, µ p = 0.4, and the restitution coefficient is taken as e n = 0.5 to account for the lubrication effect consistently with the experimental validation [START_REF] Maurin | A minimal coupled fluid-discrete element model for bedload transport[END_REF].

In order to define properly an equivalent continuous medium for the granular phase, the analysis is restricted, in a first approach, to cases with a non-negligible number of particle

τ f = max(R f xz ) = ρ f u 2 * .
Each specific density is associated with a color, which intensity reflects the Shields number. The symbol associated with the run is characteristic of the particle diameter.

layers in motion. The range of Shields numbers investigated has been chosen between θ * ∼ 0.2 and θ * ∼ 0.6, to prevent suspended load and stay in the bedload regime by keeping high suspension numbers (S * = w s /u * ∈ [2 , 5]). For mountain stream bedload, this represents intense rare events with high transport capacity. It can also be seen as a limit case of sheet flow without suspension. In order to investigate the scaling laws, the Shields number, the specific density and the particle diameter are varied. For each parameter, three values are considered starting from a realistic case for intense bedload transport (θ * ∼ 0.4, ρ p /ρ f -1 = 1.5, d = 6mm), corresponding to the particle density and diameter used for the experimental validation [START_REF] Maurin | A minimal coupled fluid-discrete element model for bedload transport[END_REF]. The specific density is lowered starting from ρ p /ρ f -1 = 1.5, considering that experimental data are often made with plastic in sheet flows (e.g. [START_REF] Capart | Transport layer structure in intense bed-load[END_REF][START_REF] Revil-Baudard | Investigation of sheet-flow processes based on novel acoustic high-resolution velocity and concentration measurements[END_REF]). The main characteristic dimensionless numbers of the problem are shown in table 1 and the exact parameters investigated are given in table 2. For each run, once the system is at steady state, the data measured every 0.1s are averaged over time for 300s for post-processing.

Methodology

In order to study the granular rheology, the granular stress tensor is computed from the 3D DEM results by applying the spatial averaging operator to the granular stress tensor. Therefore, σ p ij s is obtained for each slice of volume V by computing [START_REF] Goldhirsch | Stress, stress asymmetry and couple stress: from discrete particles to continuous fields[END_REF][START_REF] Andreotti | Granular media: between fluid and solid[END_REF])

σ p ij s = -P p δ ij + τ p ij = - 1 V β∈V m β v β i v β j - 1 V (m,n)∈V f m,n c,i b m,n j , (3.1)
where the sums are respectively over the particles and the contacts contained in the volume V , v β k = v β k -v k is the k th component of the spatial velocity fluctuation associated with particle β of mass m β , f m,n c is the contact force applied by particle m on particle n and b m,n is the branch vector from particle m to particle n.

Considering the unidirectional character of the problem, the weighting function of the spatial averaging operator extends over the whole width and length of the periodic cell. Therefore, one obtains for each run time-averaged depth profiles of the shear and normal granular stress components: τ p (z) = τ p xz (z) and P p (z) = τ p zz (z). Negligible stress asymmetry and compensated normal stress difference have been observed [START_REF] Maurin | Investigation of granular behavior in bedload transport using an Eulerian-Lagrangian model[END_REF]. Figure 2 presents typical depth profiles for a given representative simulation (run r2d6s2 in table 2). The last panel shows the different components of the streamwise mixture momentum balance, as derived in appendix A. Similarly to Revil-Baudard & Chauchat (2013) in the sheet flow regime, the viscous stress tensor contribution is found to be negligible in turbulent bedload transport. In addition, the momentum balance is closed at each elevation, as the streamwise projection of the gravity contribution is balanced by the sum of the other terms, showing the consistency of the model formulation and stress tensor evaluation. The depth profile of the inertial number (figure 2) shows that a very important range of inertial numbers (approximately five orders of magnitude) is sampled as a function of the depth in a single simulation. This is due to the nature of bedload transport in which the granular flow evolves from quasi-static in the bed to very dynamic at the granular free-surface. As the shear to normal stress ratio and the solid volume fraction vary accordingly throughout the depth (see figure 2), this allows us to obtain the granular rheology for a wide range of inertial numbers from a unique simulation. Performing simulations with variation of ρ p , d and h the water depth, the results take the form of a set of µ(I)/φ(I) curves sampling the parameter space. In the present paper, each simulation will be represented by a colored symbol, where the symbol is associated with the particle diameter d, its color is associated with the value of the specific density ρ p /ρ f -1 and its shading represents the Shields number (the darker the colour, the higher the Shields number) (see table 2).

Results

The influence of the interstitial fluid on the granular rheology varies between and within the different configurations investigated and one might expect transitions between the different regimes of the µ(I) rheology (free-fall, turbulent, viscous). While generally of minor importance, the choice for the computation of the Stokes and r dimensionless numbers (eq. 1.4 and 1.5) becomes crucial when approaching a transition between two regimes, as in the present analysis. In order to define accurately the different regimes and to impose that the transition from the viscous to the turbulent regime depends only on the particle Reynolds number (Re p = St/r), the limit turbulent drag coefficient C ∞ D = 0.44 should be considered and both turbulent and viscous drag forces should be taken without the contribution from hindrance effects. However, the classical picture of the µ(I) rheology associates the inertial number with a time-scale ratio representative of the local rearrangement process [START_REF] Andreotti | Granular media: between fluid and solid[END_REF], which is affected by hindrance effects and the particle Reynolds number in the low-Re p limit of the turbulent regime (Re p ∈ [1, 10 4 ]). Focusing on the present case (St >> 1 and r ∼ 1 so that Re p >> 1), the r and turbulent inertial numbers have been computed from a constant drag coefficient taking into account the local hindrance effects

C D (z) = C ∞ D (1 -φ(z)) -3.1
, putting aside the considerations on the transition from the viscous to the turbulent regime.

Figure 3a shows the parameter variation in the St/r plane. According to the estimation of Courrech du [START_REF] Courrech Du Pont | Granular avalanches in fluids[END_REF] (r c ∼ 6), all of the data belong to the turbulent regime and should show a collapse as a function of the turbulent inertial number. As can be seen from figure 3b, this would also be the case if the r number was evaluated without taking into account the local hindrance effect. Figure 4 and 5 show the shear to normal stress ratio and the solid volume fraction as a function of the turbulent inertial number and the dry inertial number respectively. Both figures exhibit a data collapse over the range I ∈ [10 -2 , 0.1]. At higher inertial numbers however, the data split apart as a function of the turbulent inertial number, while the collapse persists up to I dry ∼ 2 for the dry inertial number case. Therefore, the configurations sampled are better described by the dry inertial number and belong at first order to the free-fall regime, in contradiction with the predicted transition from the free-fall to the turbulent regimes of Courrech du [START_REF] Courrech Du Pont | Granular avalanches in fluids[END_REF]. One still notes a slight dependence on the specific density in figure 5, suggesting a second-order influence of the fluid inertial rearrangement mechanism. This absence of clear transition from the free-fall to the turbulent regimes is consistent with recent global analysis of immersed granular collapse using a DNS-DEM model [START_REF] Izard | Simulation of an avalanche in a fluid with a soft-sphere/immersed boundary method including a lubrication force[END_REF]. Similarly to the transition from the viscous to the free-fall regime [START_REF] Trulsson | Transition from the viscous to inertial regime in dense suspensions[END_REF], one might expect the transition region to be described by a combination of the turbulent and dry inertial numbers. In the present paper, the second-order effects are neglected for simplicity and the dry inertial number will be adopted in the following.

Coming back to figure 5 and considering the semi-logarithmic scale insets, the solid volume fraction curves are seen to collapse down to the lowest inertial number sampled (I dry ∼ 10 -5 ). For the shear to normal stress ratio, the different curves split apart below I dry ∼ 10 -2 , following different branches depending on both the specific density and the Shields number of the run considered. No variation with the particle diameter is observed. In all cases, the stress ratio values in this region are below the expected static effective friction coefficient of the granular medium (µ s ∼ 0.38 for monodisperse glass beads [START_REF] Andreotti | Granular media: between fluid and solid[END_REF]). This low-inertial-number region corresponds to the lower parts of the different simulations where a creeping regime is observed with the solid velocity following an exponential decrease [START_REF] Maurin | A minimal coupled fluid-discrete element model for bedload transport[END_REF][START_REF] Houssais | Onset of sediment transport is a continuous transition driven by fluid shear and granular creep[END_REF] similar to the one observed in dry granular flows on a heap [START_REF] Komatsu | Creep motion in a granular pile exhibiting steady surface flow[END_REF][START_REF] Richard | Rheology of confined granular flows: Scale invariance, glass transition, and friction weakening[END_REF]. These two features are characteristic of non-local effects (e.g. [START_REF] Kamrin | Nonlocal constitutive relation for steady granular flow[END_REF]; [START_REF] Bouzid | Nonlocal rheology of granular flows across yield conditions[END_REF]), where the granular flow is influenced by the far field, i.e. by the top granular flow in the present configuration. Despite the interest for granular media and out-of-equilibrium configurations, the quasi-static part of bedload transport does not contribute significantly to the sediment transport rate and will not be investigated further in this paper.

Toward the upper limit, the data are seen to collapse up to I dry ∼ 2 for both the solid volume fraction and the shear to normal stress ratio (figure 5). At higher inertial numbers, the different curves -corresponding to the different simulations -progressively exhibit a transition to a different behaviour characterised by a decrease of the shear to normal stress ratio, and a slope break on the solid volume fraction versus inertial number curves.

The observed decrease in the shear to normal stress ratio as a function of the inertial number is characteristic of the transition from dense to dilute granular behaviour [START_REF] Forterre | Flows of dense granular media[END_REF]. In the present results, the position of the transition in terms of inertial number is particularly high (I ∼ 2-3) and depends on the configurations sampled (figure 5). At such a high inertial number, the restitution coefficient is seen to have a non-negligible effect on the rheological curves (see figure 6). 2015)), the µ(I) approach is considered to break down above I dry ∼ 0.5 for glass beads, when collisional mechanisms comes into play [START_REF] Da Cruz | Macroscopic friction of dry granular materials[END_REF][START_REF] Lois | Emergence of multi-contact interactions in contact dynamics simulations of granular shear flows[END_REF][START_REF] Forterre | Flows of dense granular media[END_REF]. Therefore, it is unable to explain both the observed variation of the transition position at constant restitution coefficient and the collapse observed in regions where collisional mechanisms contribute non-negligibly to the rheology. Recently, similar persistence of dense granular flow behaviour at high inertial numbers has been observed in inclined plane configurations. [START_REF] Holyoake | High-speed granular chute flows[END_REF] observed experimentally dense granular flows up to inertial numbers I dry ∼ 2 in steady non-uniform dry granular flows of sand in a steep channel with lateral walls and a bumpy base. Moreover, analyzing granular roll wave instabilities using the DEM, [START_REF] Börzsönyi | Patterns in flowing sand: Understanding the physics of granular flow[END_REF] computed the local granular rheology and obtained a collapse of their rheological data up to I dry ∼ 1. In these two configurations as well as in the present study, the complexity of the granular flow considered can lead to differences with respect to the simple shear picture. In particular, the local analysis of spatially non-homogeneous granular flows made in the present paper and in [START_REF] Börzsönyi | Patterns in flowing sand: Understanding the physics of granular flow[END_REF], might lead to secondary gradient effects that are not taken into account in the µ(I) rheology. Therefore, the local character of the analysis and/or the specificity of the configuration seem to affect the persistence of a dense granular flow at high inertial numbers. As a consequence of the high-inertial-numbers collapse, the classical µ(I)/φ(I) expressions (eq. 1.2 and 1.3) and parametrisation do not fit the present results well at high inertial numbers (see the fit of [START_REF] Jop | A constitutive law for dense granular flows[END_REF] in figure 5). Indeed, the collisional contribution becomes important in this region and one might need to extend the µ(I)/φ(I) relationships as done by [START_REF] Holyoake | High-speed granular chute flows[END_REF] for the inclined plane configuration.

The results have evidenced the ability of the µ(I) framework to describe the granular flow rheology in bedload transport in the dense granular flow region. They have shown that the interstitial fluid does not influence the dense granular rheology importantly and that the fluid flow acts mainly as an external forcing in turbulent bedload transport. In addition, the analysis has underlined the variety and the complexity of the granular behaviours observed in bedload transport, which challenge the existing formulation of the µ(I) rheology.

Granular rheology in bedload transport

The rest of this paper focuses on the analysis of turbulent bedload transport, considering realistic conditions typical of gravel-bed rivers. Therefore, the specific density is taken as ρ p /ρ f -1 = 1.5 and additional DEM simulations are performed at Shields numbers down to the onset of motion. The analysis of the DEM results allows us to propose a µ(I)based granular rheology for turbulent bedload transport. The relevance of the proposed granular rheology is further tested in a second time, using the Eulerian-Eulerian model presented in section 2.

Parametrisation of the µ(I) rheology

Figure 7 shows the granular shear to normal stress ratio and the solid volume fraction as a function of the inertial number, for simulations with Shields numbers from θ * = 0.04 to θ * = 0.6 and for particle diameters d ∈ {3; 6; 12} mm. In these new configurations closer to realistic bedload transport, most of the data are seen to collapse for both the solid volume fraction and the granular stress ratio up to high inertial numbers. Above I ∼ 3 and under I ∼ 10 -2 , the curves split apart similarly to the general case. No dependence upon the particle diameter is observed over the whole range of inertial numbers, pointing out the negligible effect of the Re p and St numbers in the problem. The lowest-Shields-number cases (lightest points) exhibit some specific behaviour with oscillations, lower shear to normal stress ratio and solid volume fraction at given inertial numbers. These cases correspond to the sharpest transitions between a quasi-static granular bed and a few particles in motion on top of it. Therefore, in these configurations the continuous assumption is not strictly valid and the limits of the rheological description are reached.

Overriding these limitations close to the critical Shields number and assuming that the high-inertial-numbers granular behaviour is not dominant in bedload transport, the µ(I) constitutive laws are pragmatically best-fitted on the Eulerian-Lagrangian simulation results over the range of inertial numbers I dry ∈ [10 -2 , 3]. The classical solid volume fraction expression (eq. 1.3) predicts negative values at high inertial numbers and is unable to describe the present results. Consequently, the expression is modified and is chosen following [START_REF] Aussillous | Investigation of the mobile granular layer in bedload transport by laminar shearing flows[END_REF]Revil-Baudard &Chauchat (2013) as

φ(I) = φ max 1 + bI . (4.1)
The best fits of equations 1. (2011) and[START_REF] Trulsson | Transition from the viscous to inertial regime in dense suspensions[END_REF] in the viscous regime. The values of µ 2 and I 0 are logically higher than in the original µ(I) rheology, as no real saturation of the stress ratio is observed in our data. The best fits describe very well the results obtained within the collapsing range of I dry , as can be seen in figure 7. Outside this interval, they do not reproduce the non-local effects at low inertial numbers (I dry < 10 -2 ) and the drop in both the shear to normal stress ratio and the solid volume fraction for I dry > 3. Assuming the dense granular behaviour to be dominant, the two relationships (eq. 1.2-4.1) together with the values determined from the best fits are proposed as a µ(I)-based granular rheology for turbulent bedload transport.

Eulerian-Eulerian simulations

In order to evaluate the relevance of the proposed granular rheology it has been implemented in the Eulerian-Eulerian model presented in section 2. Simulations have been performed using exactly the same parameters for both the Eulerian-Lagrangian and the Eulerian-Eulerian models (ρ p /ρ f -1 = 1.5 ; d = 6 mm, α = 0.05, h). The particle diameter has been varied (d ∈ 3, 6, 12mm) but does not influence the results so that only the case d = 6mm is presented for clarity. experimental data is obtained using both models. The proposed granular rheology enables the Eulerian-Eulerian model to accurately reproduce the critical Shields number (θ * c ∼ 0.04) and the global trend over a wide range of Shields numbers (θ * ∈ [θ * c , 0.7]). Quite surprinsingly, the Eulerian-Eulerian model better predicts the smooth transition around the critical Shields number, in a region where the continuous assumption breaks down. Conversely, at high Shields numbers the increased importance of the dilute region, which is not well described by the proposed rheology (section 4.1), induces a slight discrepancy between the Eulerian-Eulerian and Eulerian-Lagragian models.

To go further in the analysis, the comparison is extended to the granular depth profiles by decomposing the sediment transport rate into the integral of the sediment transport rate density q s (z) = v x s (z) φ(z) [START_REF] Bagnold | The flow of cohesionless grains in fluids[END_REF]:

Q s = ∞ 0 q s (z)dz = ∞ 0 v x s (z) φ(z)dz. (4.2)
Figure 9 shows the comparison of the solid volume fraction, the solid velocity and the transport rate density profiles obtained with the Eulerian-Eulerian and Eulerian-Lagrangian models for three different Shields numbers θ * = {0.05; 0.1; 0.2}. The Eulerian-Eulerian model reproduces accurately the granular depth profiles obtained with the Eulerian-Lagrangian model for solid volume fractions approximately higher than 0.3 (φ 0.3). This value roughly corresponds to the domain for which the proposed granular rheology describes the Eulerian-Lagrangian results very well (see figure 7). Only the lower quasi-static velocity profile -negligible in terms of sediment transport rate -is not reproduced, as non-local effects have not been taken into account in the granular rheology. In the dilute region (φ 0.3), the particle velocity is overestimated by the Eulerian-Eulerian model for the three different Shields numbers, showing an increasing discrepancy toward the lower solid volume fraction region. This is associated with an under-estimation of the solid volume fraction close to incipient motion and an overestimation at high Shield numbers. These behaviours point out both the ill-posedness of the µ(I) rheology as the solid volume fraction goes to zero and the inability of the proposed rheology to describe the dilute region accurately. Despite these limitations, the good lower boundary condition given by the accurate description of the dense granular flow constrains the dilute behaviour and enables a good description of the sediment transport rate over a wide range of Shields numbers with the Eulerian-Eulerian model.

The present results show that a µ(I) rheology is able to reproduce the main features of turbulent bedload transport in terms of sediment transport rate and dense granular depth profiles, in steady uniform configurations. Considering the much lower computational cost of Eulerian-Eulerian simulations compared with Eulerian-Lagrangian ones, it will allow an important step to be made in the upscaling process towards applications. In this perspective, the fluid phase description might have to be upgraded, using a two-equation turbulence model as k -or k -ω [START_REF] Hsu | On two-phase sediment transport: sheet flow of massive particles[END_REF][START_REF] Amoudry | Extension of turbulence closure to two-phase sediment transport modelling: Application to oscillatory sheet flows[END_REF][START_REF] Lee | Multi-dimensional rheologybased two-phase model for sediment transport and applications to sheet flow and pipeline scour[END_REF], or eddy resolving simulations.

Summary and conclusion

The granular rheology in idealised turbulent bedload transport has been studied using a coupled fluid-discrete-element model. Computing the granular stress tensor locally as a function of the depth and performing simulations for various Shields numbers, specific densities and particle diameters, the data have been shown to collapse over an important range of inertial numbers, demonstrating the relevance of the µ(I) rheology for turbulent bedload transport description. No clear transition from the free-fall to the turbulent regime of the µ(I) rheology has been observed, supporting recent numerical observations and suggesting that the interstitial fluid does not importantly influence the granular rheology in turbulent bedload transport. In addition, consistently with recent literature, the data have been observed to collapse in terms of both the solid volume fraction and the shear to normal stress ratio, up to inertial numbers as high as I ∼ 2. The persistence of a dense granular flow behaviour at such a high inertial number together with the observed progressive transition from dense to dilute behaviour challenges the classical conceptions and parametrizations of the µ(I) rheology, and opens interesting perspectives for a better understanding of granular media rheology at high inertial numbers. Beyond these observations, a µ(I)-based granular rheology has been fitted to the present Eulerian-Lagrangian bedload transport simulations. The proposed rheology has been sucessfully tested using a Eulerian-Eulerian model. It has been shown to accurately reproduce the dense granular depth profiles and the classical experimental results in terms of sediment transport rate from the onset of motion up to a Shields number of θ * ∼ 0.7. Further work is needed to better understand the high-inertial-numbers behaviour and extend the granular rheology to the description of the dilute region. At this stage it is not clear whether the µ(I) rheology could be extended to handle this regime transition. Besides, the proposed rheology can already be used in three dimensional Eulerian-Eulerian numerical models and represents a step in the upscaling process from particlescale simulations toward full-scale problem modelling.

Figure 1 .

 1 Figure 1. Scheme of the numerical setup and its equivalent average unidirectional picture with typical fluid velocity u f = u fx (z) ex, solid volume fraction φ, and solid velocity v s = vx s (z) ex depth profiles. The inclined 3D bi-periodic granular description is coupled with a unidirectional fluid momentum balance using imposed fixed random bottom boundary condition and water free-surface position. The particle color is representative of the velocity intensity.

Figure 2 .

 2 Figure 2. Representative example (case r2d6s2 in table2) of solid volume fraction, shear to normal granular stress ratio, inertial number, and total stress repartition depth profiles. The last panel shows the different components of the mixture momentum balance as a function of the depth, as derived in appendix A: the granular shear stress (-), the Reynolds shear stress (-), the viscous shear stress (--), the sum of the three (-), and the slope contribution (last term of equation A 4, --). The water free-surface elevation is represented in each different panel for reference (-).

Figure 3 .

 3 Figure 3. Position of all the computed granular rheology points in the (St,r) plane, considering the drag coefficient taking into account the local hindrance effects (a) and a constant drag coefficient respectively. The dashed lines represent the transition between the different regimes, with a shaded area between the expected transition from turbulent to free-fall regime rc = 1 and the one estimated by Courrech du Pont et al. (2003) rc = 6.

Figure 4 .

 4 Figure4. Shear to normal stress ratio µ = τ p /P p and solid volume fraction φ as a function of the turbulent inertial number I turb = γd ρ f CD/P p , for all the cases presented in table 2 with variation of Shields number, specific density and particle diameter.

Figure 5 .

 5 Figure 5. Shear to normal stress ratio µ = τ p /P p and solid volume fraction φ as a function of the dry inertial number I dry = γd ρ p /P p for all the cases presented in table 2 with variation of Shields number, specific density and particle diameter. The parameters of the simulation sampled and the corresponding symbols are shown in table 2. The dashed lines (-) represents the classical expression of µ(I)/φ(I) (eq. 1.2 and 1.3) with the parametrization of Da Cruz et al. (2005) and Jop et al. (2006): µs = 0.38, µ2 = 0.64, I0 = 0.279, a = 0.31.

  Focusing on the classical dry granular flow literature (e.g. Da Cruz et al. (2005); Jop et al. (2006); Forterre & Pouliquen (2008); Jop (

Figure 6 .

 6 Figure6. Effect of the restitution coefficient on the shear to normal stress ratio µ = τ p /P p and solid volume fraction φ as a function of the dry inertial number I dry = γd ρ p /P p for a representative case (case r2d6s2 in table 2).

Figure 7 .

 7 Figure 7. Shear to normal stress ratio and solid volume fraction as a function of the dry inertial number, for realistic bedload transport simulations with specific density ρ p /ρ f -1 = 1.5, Shields number from incipient motion to 0.7, and particle diameter of d = 3mm, d = 6mm, and d = 12mm. The black lines represent the best fit obtained with equations 1.2 and 4.1 (µ1 = 0.38, µ2 = 0.64, I0 = 0.279, φ max = 0.61, and b = 0.31). The symbols colors are function of the simulation Shields number: the darkest the color, the highest the Shields number.

Figure 8 .

 8 Figure 8. Dimensionless sediment transport rate as a function of the Shields number for two-phase continuous simulations (Eulerian-Eulerian) and coupled fluid-discrete element simulations (Eulerian-Lagrangian), with realistic parameters for turbulent bedload transport: ρ p /ρ f -1 = 1.5 and d = 6mm. The black crosses and the red line represent the data and the fit of Meyer-Peter & Müller (1948) respectively.

Figure 8 Figure 9 .

 89 Figure 8 presents the dimensionless sediment transport rate versus the Shields number for the Eulerian-Eulerian model, the Eulerian-Lagrangian model and the experimental data and empirical law of Meyer-Peter & Müller (1948). A good agreement with the

Table 2 .

 2 Parameters of the simulations studied and symbol correspondence. The specific density, particle diameter, and Shields number have been varied. The latter is evaluated from the maximum of the turbulent shear stress, using

	Run	ρ p /ρ f -1 d(mm) θ * symbol
	r0d3s1	0.375	3	0.214	•
	r0d3s2	0.375	3	0.444	•
	r0d3s3	0.375	3	0.611	•
	r0d6s1	0.375	6	0.215	+
	r0d6s2	0.375	6	0.44	+
	r0d6s3	0.375	6	0.598	+
	r0d12s1	0.375	12	0.216	x
	r0d12s2	0.375	12	0.434	x
	r0d12s3	0.375	12	0.593	x
	r1d3s1	0.75	3	0.188	•
	r1d3s2	0.75	3	0.378	•
	r1d3s3	0.75	3	0.593	•
	r1d6s1	0.75	6	0.193	+
	r1d6s2	0.75	6	0.381	+
	r1d6s3	0.75	6	0.598	+
	r1d12s1	0.75	12	0.191	x
	r1d12s2	0.75	12	0.379	x
	r1d12s3	0.75	12	0.596	x
	r2d3s1	1.5	3	0.205	•
	r2d3s2	1.5	3	0.443	•
	r2d3s3	1.5	3	0.694	•
	r2d6s1	1.5	6	0.205	+
	r2d6s2	1.5	6	0.455	+
	r2d6s3	1.5	6	0.692	+
	r2d12s1	1.5	12	0.21	x
	r2d12s2	1.5	12	0.451	x
	r2d12s3	1.5	12	0.696	x

Aknowledgement

We are greatful to the three anonymous reviewers for relevant and constructive questions, and to Ashley Dudill for English corrections. RM would like to thank Thierry Faug for fruitful discussions. This research was supported by Irstea (formerly Cemagref), the labex OSUG@2020 and the French Institut National des Sciences de l'Univers programs EC2CO-BIOHEFECT and EC2CO-LEFE MODSED.

Appendix A. Momentum balance derivation

In order to validate the stress tensor profile computed, the momentum balance is analyzed in the framework of the continuous two-phase equations [START_REF] Jackson | The dynamics of fluidized particles[END_REF]. In the present appendix, the equations are derived and integrated to show the exact terms computed in the body of the paper. Using the steady and uniform character of the problem considered, the momentum balances for both the fluid and the granular phases read along the streamwise direction [START_REF] Jackson | The dynamics of fluidized particles[END_REF] Combining these two equations together, the mixture momentum balance can be written as

In order to study the stress repartition, the equation is integrated between a given position z in h > z > 0 and the water free-surface elevation h, where the viscous, turbulent and particle shear stresses vanish:

From one simulation, it is possible to evaluate all of the terms of the equation at each given elevation. Provided that the stress tensor formulation is appropriate and the system is at equilibrium, the equality should be satisfied.