2.1 Discussion on the classical critical Shields number derivation

The classical derivation of the critical Shields number reproduced in the introduction relies on the following expression of the buoyancy force applied to a particle:

(2.1) $$\begin{eqnarray}\boldsymbol{f}_{b}=-\unicode[STIX]{x1D70C}^{f}\frac{\unicode[STIX]{x03C0}}{6}d^{3}\boldsymbol{g}.\end{eqnarray}$$

While this expression is classically used, it does not apply a priori to all the different configurations explored in fluid–particle flow, as stressed by Christensen (Reference Christensen1995).

The buoyancy force is defined as the force a fluid element would undergo if it was occupying the position of the particle (Maxey & Riley Reference Maxey and Riley1983). It can be derived explicitly in the Stokes flow case and leads to the following formulation (Maxey & Riley Reference Maxey and Riley1983):

(2.2) $$\begin{eqnarray}\boldsymbol{f}_{b}=-\frac{\unicode[STIX]{x03C0}}{6}d^{3}\unicode[STIX]{x1D735}\boldsymbol{\cdot }\unicode[STIX]{x1D70E}_{u}^{f},\end{eqnarray}$$

where $\unicode[STIX]{x1D70E}_{u}^{f}$ is the undisturbed fluid velocity stress tensor, i.e. the fluid stress tensor based on the fluid velocity field undisturbed by the presence of the particles. This formulation is used in Jackson (Reference Jackson2000) at the level of the continuous medium, including the Reynolds stresses contribution ( $R^{f}$ ) inside the fluid stress tensor, i.e. replacing $\unicode[STIX]{x1D70E}_{u}^{f}$ by $\unicode[STIX]{x1D70E}_{u}^{f}+R^{f}$ inside (2.2). However, the formulation of the buoyancy force should not a priori include the contribution from the Reynolds stresses considering that they are associated with advection and related to the time averaging of the fluid velocity fluctuations. Indeed, considering a single particle in a turbulent Newtonian fluid flow, it experiences a time-averaged buoyancy force depending on the average fluid velocity field but independent of the fluctuations with time. This is the case because the expression of the buoyancy is linear with the second spatial derivative of the fluid velocity for a Newtonian fluid,

(2.3) $$\begin{eqnarray}\boldsymbol{f}_{b}=-\frac{\unicode[STIX]{x03C0}}{6}d^{3}\unicode[STIX]{x1D702}^{f}\unicode[STIX]{x1D735}\boldsymbol{\cdot }[\unicode[STIX]{x1D735}\boldsymbol{u}^{\boldsymbol{f}}+(\unicode[STIX]{x1D735}\boldsymbol{u}^{\boldsymbol{f}})^{\text{T}}],\end{eqnarray}$$

with $\unicode[STIX]{x1D702}^{f}$ the dynamic fluid viscosity, so that the time averaging of the buoyancy force reduces to the contribution from the average fluid velocity field. From these considerations, the expression of the buoyancy force should follow (2.2) in all the different fluid flow regimes (D. Lhuillier, 2011, Personal communication).

As a consequence, the above-mentioned derivation is not valid in uniform gravity-driven turbulent bedload transport as the undisturbed fluid equation, $0=-\unicode[STIX]{x1D735}R^{f}-\unicode[STIX]{x1D735}\unicode[STIX]{x1D70E}_{u}^{f}+\unicode[STIX]{x1D70C}^{f}\boldsymbol{g}$ , together with the regular expression of the buoyancy force (2.2), does not allow us to write the buoyancy force as (2.1).

Expressing the buoyancy directly from (2.2), it follows that $(f_{b})_{z}=\unicode[STIX]{x03C0}d^{3}/6\,\unicode[STIX]{x1D70C}^{f}g\cos \unicode[STIX]{x1D6FC}$ and $(f_{b})_{x}=-\unicode[STIX]{x03C0}d^{3}/6\,\unicode[STIX]{x2202}_{z}\unicode[STIX]{x1D70E}_{xz}^{f}$ , by using the hydrostatic characteristics of the fluid along $z$ . Considering that the viscous shear stress contribution is negligible in turbulent bedload transport $\unicode[STIX]{x2202}_{z}\unicode[STIX]{x1D70E}_{xz}^{f}\simeq 0$ (Revil-Baudard & Chauchat Reference Revil-Baudard and Chauchat2013; Maurin, Chauchat & Frey Reference Maurin, Chauchat and Frey2016), the force balance performed on a particle lying on the granular bed at the onset of motion (1.3) reads:

(2.4) $$\begin{eqnarray}\frac{\unicode[STIX]{x03C0}}{8}\unicode[STIX]{x1D70C}^{f}d^{2}C_{D}u_{\ast }^{2}-\unicode[STIX]{x1D707}_{s}(\unicode[STIX]{x1D70C}^{p}-\unicode[STIX]{x1D70C}^{f})g\cos \unicode[STIX]{x1D6FC}\frac{\unicode[STIX]{x03C0}}{6}d^{3}+\unicode[STIX]{x1D70C}^{p}g\sin \unicode[STIX]{x1D6FC}\frac{\unicode[STIX]{x03C0}}{6}d^{3}=0,\end{eqnarray}$$

where only the particle weight projection along the streamwise axis is non-negligible. This leads to the following modification of the critical Shields number due to the slope effect:

(2.5) $$\begin{eqnarray}\unicode[STIX]{x1D703}_{c}(\unicode[STIX]{x1D6FC})=\unicode[STIX]{x1D703}_{c}^{0}\cos \unicode[STIX]{x1D6FC}\left[1-\frac{\tan \unicode[STIX]{x1D6FC}}{\unicode[STIX]{x1D707}_{s}}\frac{1}{1-\unicode[STIX]{x1D70C}^{f}/\unicode[STIX]{x1D70C}^{p}}\right],\end{eqnarray}$$

where the density ratio $\unicode[STIX]{x1D70C}^{p}/\unicode[STIX]{x1D70C}^{f}$ modifies the classical slope effect characterised by (1.4). This difference is of importance for underwater bedload transport, where the density ratio is typically of order one ( $\unicode[STIX]{x1D70C}^{p}\gtrsim \unicode[STIX]{x1D70C}^{f}$ ).

Adopting the regular expression for the buoyancy force also affects the derivation of the slope influence on the critical Shields number in the pressure-driven configuration. In this case, the streamwise component of the buoyancy force is completed by the term due to the presence of a longitudinal pressure gradient:

(2.6) $$\begin{eqnarray}\frac{\unicode[STIX]{x03C0}}{8}\unicode[STIX]{x1D70C}^{f}d^{2}C_{D}u_{\ast }^{2}-\unicode[STIX]{x1D707}_{s}(\unicode[STIX]{x1D70C}^{p}-\unicode[STIX]{x1D70C}^{f})g\cos \unicode[STIX]{x1D6FC}\frac{\unicode[STIX]{x03C0}}{6}d^{3}+\unicode[STIX]{x1D70C}^{p}g\sin \unicode[STIX]{x1D6FC}\frac{\unicode[STIX]{x03C0}}{6}d^{3}-\frac{\unicode[STIX]{x2202}P^{f}}{\unicode[STIX]{x2202}x}=0,\end{eqnarray}$$

which does not enable us to recover the classical modification of the critical Shields number derived in the introduction (1.4). The latter equation remains only valid for laminar flows and in the case where the fluid mass column stays at a constant level perpendicular to the gravity, i.e. for underwater avalanches or coastal sediment transport.

2.2 Vertical granular flow structure and entrainment mechanisms

Instead of focusing on the critical Shields number modification in a discrete particle framework, we adopt a more general approach based on the continuous two-phase flow framework to analyse the slope influence from the onset of motion up to intense turbulent bedload transport. In order to better understand the vertical structure of the granular flow and the local mechanisms at play, let us express the shear-to-normal granular stress ratio as a function of the depth $\unicode[STIX]{x1D70F}_{xz}^{p}(z)/P^{p}(z)$ . The yield criterion $\unicode[STIX]{x1D70F}_{xz}^{p}/P^{p}>\unicode[STIX]{x1D707}_{s}$ is characteristic of the onset of granular flow, so that the positive contributions to the shear-to-normal stress ratio at a given height $z$ reflects the local entrainment mechanism of the granular medium.

The two-phase volume-averaged equations read for steady uniform flows (Anderson & Jackson Reference Anderson and Jackson1967; Jackson Reference Jackson2000; Chauchat Reference Chauchat2017):

(2.7) $$\begin{eqnarray}\displaystyle & \displaystyle 0=\frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D70F}_{xz}^{f}}{\unicode[STIX]{x2202}z}+\frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D619}_{xz}^{f}}{\unicode[STIX]{x2202}z}+\unicode[STIX]{x1D70C}^{f}(1-\unicode[STIX]{x1D719})g\sin \unicode[STIX]{x1D6FC}-n\langle \,{f_{f}^{p}}_{x}\rangle ^{p}, & \displaystyle\end{eqnarray}$$

(2.8) $$\begin{eqnarray}\displaystyle & \displaystyle 0=\frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D70F}_{xz}^{p}}{\unicode[STIX]{x2202}z}+\unicode[STIX]{x1D70C}^{p}\unicode[STIX]{x1D719}g\sin \unicode[STIX]{x1D6FC}+n\langle \,{f_{f}^{p}}_{x}\rangle ^{p}, & \displaystyle\end{eqnarray}$$

(2.9) $$\begin{eqnarray}\displaystyle & \displaystyle 0=-\frac{\unicode[STIX]{x2202}P^{f}}{\unicode[STIX]{x2202}z}+\unicode[STIX]{x1D70C}^{f}(1-\unicode[STIX]{x1D719})g\cos \unicode[STIX]{x1D6FC}-n\langle \,{f_{f}^{p}}_{z}\rangle ^{p}, & \displaystyle\end{eqnarray}$$

(2.10) $$\begin{eqnarray}\displaystyle & \displaystyle 0=-\frac{\unicode[STIX]{x2202}P^{p}}{\unicode[STIX]{x2202}z}+\unicode[STIX]{x1D70C}^{p}\unicode[STIX]{x1D719}g\cos \unicode[STIX]{x1D6FC}+n\langle \,{f_{f}^{p}}_{z}\rangle ^{p}, & \displaystyle\end{eqnarray}$$

where $\unicode[STIX]{x1D70E}_{ij}^{f}=-P^{f}\unicode[STIX]{x1D6FF}_{ij}+\unicode[STIX]{x1D70F}_{ij}^{f}$ is the effective fluid stress tensor with $\unicode[STIX]{x1D6FF}_{ij}$ Kronecker’s delta, $\unicode[STIX]{x1D70E}_{ij}^{p}=-P^{p}\unicode[STIX]{x1D6FF}_{ij}+\unicode[STIX]{x1D70F}_{ij}^{p}$ is the granular stress tensor, $\unicode[STIX]{x1D619}_{ij}^{f}$ is the Reynolds stress tensor, $\unicode[STIX]{x1D719}$ and $\unicode[STIX]{x1D716}$ are the solid and fluid volume fraction respectively and $n\langle \,{f_{f}^{p}}_{k}\rangle ^{p}$ is the $k$ component of the fluid–solid momentum transfer, with $n$ the particle density and $\langle \,{f_{f}^{p}}_{k}\rangle ^{p}$ the volume-averaged fluid–particle interaction force. These equations represents the volume-averaged momentum balance over the streamwise and wall-normal directions of the fluid and solid phases, and the brackets with superscript $p$ denotes the classical volume averaging over the particle phase defined in Anderson & Jackson (Reference Anderson and Jackson1967).

Integrating these equations between a position $z<h_{p}$ in the moving granular layer and $h_{p}$ , the top of it, allows us to simplify the equations and to express the shear-to-normal stress ratio as a function of the vertical position $z$ . Focusing on the granular equations (2.8) and (2.10), the granular stress tensor is zero at $h_{p}$ and the integrated streamwise component of the granular momentum balance leads to:

(2.11) $$\begin{eqnarray}\unicode[STIX]{x1D70F}_{xz}^{p}(z)=\unicode[STIX]{x1D70C}^{p}g\sin \unicode[STIX]{x1D6FC}\int _{z}^{h_{p}}\unicode[STIX]{x1D719}(\unicode[STIX]{x1D701})\,\text{d}\unicode[STIX]{x1D701}+\int _{z}^{h_{p}}n\langle \,{f_{f}^{p}}_{x}\rangle ^{p}(\unicode[STIX]{x1D701})\,\text{d}\unicode[STIX]{x1D701}.\end{eqnarray}$$

Assuming that the fluid wall-normal pressure is hydrostatic, the fluid–particle interaction term $n\langle \,{f_{f}^{p}}_{z}\rangle ^{p}$ can be expressed from (2.9). Replacing it in (2.10), we obtain:

(2.12) $$\begin{eqnarray}P^{p}(z)=g\cos \unicode[STIX]{x1D6FC}(\unicode[STIX]{x1D70C}^{p}-\unicode[STIX]{x1D70C}^{f})\int _{z}^{h_{p}}\unicode[STIX]{x1D719}(\unicode[STIX]{x1D701})\,\text{d}\unicode[STIX]{x1D701}.\end{eqnarray}$$

Lastly, the fluid streamwise momentum balance (2.7) can be integrated between $z$ and $h_{p}$ in order to express the last term on the right-hand side of (2.11):

(2.13) $$\begin{eqnarray}0=\unicode[STIX]{x1D70F}_{xz}^{f}(h_{p})-\unicode[STIX]{x1D70F}_{xz}^{f}(z)+R_{xz}^{f}(h_{p})-R_{xz}^{f}(z)+\unicode[STIX]{x1D70C}^{f}g\sin \unicode[STIX]{x1D6FC}\int _{z}^{h_{p}}[1-\unicode[STIX]{x1D719}(\unicode[STIX]{x1D701})]\,\text{d}\unicode[STIX]{x1D701}-\int _{z}^{h_{p}}n\langle \,{f_{f}^{p}}_{x}\rangle ^{p}(\unicode[STIX]{x1D701})\,\text{d}\unicode[STIX]{x1D701}.\end{eqnarray}$$

It has been shown previously in turbulent bedload transport and sheet-flow simulations that the effective viscous shear stress tensor is negligible throughout the depth with respect to the Reynolds stresses and to the slope contribution (Revil-Baudard & Chauchat Reference Revil-Baudard and Chauchat2013; Maurin et al. Reference Maurin, Chauchat and Frey2016). Furthermore, the Reynolds stresses are fully damped in the granular bed, and we assume that their contribution is negligible in the moving granular layer $R_{xz}^{f}(z)\sim 0$ . Also, the fluid bed shear stress used to define the Shields number is taken here as the maximum of the turbulent shear stress (Revil-Baudard & Chauchat Reference Revil-Baudard and Chauchat2013; Maurin et al. Reference Maurin, Chauchat and Frey2016), assuming that the latter is located at the top of the granular layer, $\unicode[STIX]{x1D70F}_{b}=R_{xz}^{f}(h_{p})$ . Therefore, equation (2.13) reads:

(2.14) $$\begin{eqnarray}\int _{z}^{h_{p}}n\langle \,{f_{f}^{p}}_{x}\rangle ^{p}(\unicode[STIX]{x1D701})\,\text{d}\unicode[STIX]{x1D701}=\unicode[STIX]{x1D70F}_{b}+\unicode[STIX]{x1D70C}^{f}g\sin \unicode[STIX]{x1D6FC}\int _{z}^{h_{p}}[1-\unicode[STIX]{x1D719}(\unicode[STIX]{x1D701})]\,\text{d}\unicode[STIX]{x1D701}.\end{eqnarray}$$

Combining equations (2.11), (2.13) and (2.14) and defining the average solid volume fraction between $z$ and $h_{p}$ , $\bar{\unicode[STIX]{x1D719}}_{z}$ as:

(2.15) $$\begin{eqnarray}\bar{\unicode[STIX]{x1D719}}_{z}(h_{p}-z)=\int _{z}^{h_{p}}\unicode[STIX]{x1D719}(\unicode[STIX]{x1D701})\,\text{d}\unicode[STIX]{x1D701},\end{eqnarray}$$

the shear-to-normal stress ratio $\unicode[STIX]{x1D70F}_{xz}^{p}(z)/\unicode[STIX]{x1D70F}_{zz}^{p}(z)$ can be expressed as:

(2.16) $$\begin{eqnarray}\displaystyle \frac{\unicode[STIX]{x1D70F}_{xz}^{p}(z)}{\unicode[STIX]{x1D70F}_{zz}^{p}(z)}=\frac{\unicode[STIX]{x1D70C}^{p}}{\unicode[STIX]{x1D70C}^{p}-\unicode[STIX]{x1D70C}^{f}}\tan \unicode[STIX]{x1D6FC}+\frac{\unicode[STIX]{x1D70F}_{b}}{g\cos \unicode[STIX]{x1D6FC}(\unicode[STIX]{x1D70C}^{p}-\unicode[STIX]{x1D70C}^{f})\bar{\unicode[STIX]{x1D719}}_{z}(h_{p}-z)}+\frac{\unicode[STIX]{x1D70C}^{f}}{\unicode[STIX]{x1D70C}^{p}-\unicode[STIX]{x1D70C}^{f}}\frac{1-\bar{\unicode[STIX]{x1D719}}_{z}}{\bar{\unicode[STIX]{x1D719}_{z}}}\tan \unicode[STIX]{x1D6FC}. & & \displaystyle\end{eqnarray}$$

Three contributions to the granular phase shear-to-normal stress ratio can be identified. The first term on the right-hand side of (2.16) is constant within the depth, and represents the slope effect on the granular phase. It is analogous to the classical configuration of a dry granular medium on an inclined plane, where the shear-to-normal stress ratio is constant within the depth and equal to the tangent of the inclination angle $\unicode[STIX]{x1D6FC}$ (Andreotti et al. Reference Andreotti, Forterre and Pouliquen2013). In bedload transport, the presence of a buoyancy force along the vertical axis leads to a modification of this term by the density ratio. The two last terms on the right-hand side correspond to the fluid contribution on the granular phase, which has been split into two contributions coming from the fluid flow above the bed and the fluid flow inside the granular layer respectively. The Shields number can be identified in the former and the term reduces to:

(2.17) $$\begin{eqnarray}\frac{\unicode[STIX]{x1D70F}_{b}}{g\cos \unicode[STIX]{x1D6FC}(\unicode[STIX]{x1D70C}^{p}-\unicode[STIX]{x1D70C}^{f})\bar{\unicode[STIX]{x1D719}}_{z}(h_{p}-z)}=\frac{\unicode[STIX]{x1D703}^{\ast }}{\bar{\unicode[STIX]{x1D719}}_{z}}\frac{d}{h_{p}-z}.\end{eqnarray}$$

This represents the contribution from the fluid bed shear stress to the entrainment of particles. Its surface nature is characterised by a decrease of this contribution inside the granular bed, i.e. as $z$ decreases. It is accounted for in bedload transport studies, by assuming that the Shields number is representative of the total imposed fluid shear stress. However, the contribution from the fluid flow inside the granular layer, corresponding to the third term on the right-hand side of (2.16), is usually neglected while it is seen to be of importance when increasing the slope angle. This contribution depends on the slope and characterises the effect of the slope on the fluid flow inside the granular layer. The latter affects the granular phase through the fluid–particle interaction, and in consequence depends also on the specific density. In addition, the different dependence on the vertical position $z$ of the two fluid terms implies that a modification of the slope at given Shields number and specific density would induce a change in the vertical repartition of the shear-to-normal stress ratio, and therefore on the vertical structure of the granular flow. This could explain why the velocity and solid volume fraction profiles are not self-similar in gravity-driven bedload transport (Armanini et al. Reference Armanini, Capart, Fraccarollo and Larcher2005; Larcher et al. Reference Larcher, Fraccarollo, Armanini and Capart2007; Capart & Fraccarollo Reference Capart and Fraccarollo2011; Frey Reference Frey2014; Maurin Reference Maurin2015; Maurin et al. Reference Maurin, Chauchat, Chareyre and Frey2015).

In order to gain more insight into the physical meaning of the fluid flow contribution inside the granular bed, let us consider the simple case of a gravity-driven fluid flow through a quasi-static porous granular bed. The Shields number contribution is negligible in this configuration and the granular bed is at rest so that the average solid volume fraction is equal to the maximum packing fraction $\bar{\unicode[STIX]{x1D719}}_{z}=\unicode[STIX]{x1D719}^{max}$ . This leads to a shear-to-normal stress ratio independent of the vertical position  $z$ :

(2.18) $$\begin{eqnarray}\displaystyle \frac{\unicode[STIX]{x1D70F}_{xz}^{p}(z)}{P^{p}(z)}=\frac{\unicode[STIX]{x1D70C}^{p}}{\unicode[STIX]{x1D70C}^{p}-\unicode[STIX]{x1D70C}^{f}}\tan \unicode[STIX]{x1D6FC}+\frac{\unicode[STIX]{x1D70C}^{f}}{\unicode[STIX]{x1D70C}^{p}-\unicode[STIX]{x1D70C}^{f}}\frac{1-\unicode[STIX]{x1D719}^{max}}{\unicode[STIX]{x1D719}^{max}}\tan \unicode[STIX]{x1D6FC}=\frac{\tan \unicode[STIX]{x1D6FC}}{\unicode[STIX]{x1D70C}^{p}/\unicode[STIX]{x1D70C}^{f}-1}\left[\frac{\unicode[STIX]{x1D70C}^{p}}{\unicode[STIX]{x1D70C}^{f}}+\frac{1-\unicode[STIX]{x1D719}^{max}}{\unicode[STIX]{x1D719}^{max}}\right].\qquad & & \displaystyle\end{eqnarray}$$

As a consequence, there exists a critical angle above which the entire granular layer will be entrained independently from its thickness. This avalanche angle, $\unicode[STIX]{x1D6FC}_{0}$ , corresponds to the configuration for which the shear-to-normal stress ratio exceeds the static friction coefficient of the granular medium  $\unicode[STIX]{x1D707}_{s}$ :

(2.19) $$\begin{eqnarray}\tan \unicode[STIX]{x1D6FC}_{0}=\frac{\unicode[STIX]{x1D707}_{s}}{1+[(\unicode[STIX]{x1D70C}^{p}/\unicode[STIX]{x1D70C}^{f}-1)\unicode[STIX]{x1D719}^{max}]^{-1}}.\end{eqnarray}$$

This corresponds to the critical angle for the onset of debris flow as derived by Takahashi (Reference Takahashi1978), and is well known in the debris flow community (Takahashi Reference Takahashi2007). The derivation presented here makes the link between debris flows and gravity-driven turbulent bedload transport, and evidences the importance of the slope influence through the fluid flow inside the granular layer when considering steep slope configurations. In addition, it shows that the expression of the critical angle of debris flow $\unicode[STIX]{x1D6FC}_{0}$ is characteristic of the slope influence on the granular phase through both its direct impact and the fluid-induced contribution in gravity-driven bedload transport, two mechanisms that depend on the specific density $\unicode[STIX]{x1D70C}^{p}/\unicode[STIX]{x1D70C}^{f}-1$ .

2.3 Sediment transport rate scaling

Taking advantage of the developed analysis, we look for a scaling law of the sediment transport rate accounting for the slope effect. Following Bagnold (Reference Bagnold1956), the sediment transport rate per unit width $Q_{s}$ can be expressed as a function of the solid volume fraction $\unicode[STIX]{x1D719}$ and the average particle velocity $\langle v_{x}\rangle ^{p}$ profiles:

(2.20) $$\begin{eqnarray}Q_{s}=\int _{0}^{\infty }\unicode[STIX]{x1D719}(z)\langle v_{x}\rangle ^{p}(z)\,\text{d}z=\int _{h_{c}}^{h_{p}}\unicode[STIX]{x1D719}(z)\langle v_{x}\rangle ^{p}(z)\,\text{d}z=\overline{\unicode[STIX]{x1D719}\langle v_{x}\rangle ^{p}}\unicode[STIX]{x1D6FF}_{s},\end{eqnarray}$$

where $h_{c}$ is the lower bound of the mobile granular layer, $\unicode[STIX]{x1D6FF}_{s}=h_{p}-h_{c}$ is the mobile granular layer thickness and $\bar{\bullet }$ denotes the average over $z$ along the mobile granular layer thickness. Considering $z=h_{c}$ in (2.16), the mobile layer thickness can be expressed as:

(2.21) $$\begin{eqnarray}\frac{\unicode[STIX]{x1D6FF}_{s}}{d}=\frac{\unicode[STIX]{x1D703}^{\ast }}{\bar{\unicode[STIX]{x1D719}}[\unicode[STIX]{x1D707}_{s}-(1+[\bar{\unicode[STIX]{x1D719}}(\unicode[STIX]{x1D70C}^{p}/\unicode[STIX]{x1D70C}^{f}-1)]^{-1})\tan \unicode[STIX]{x1D6FC}]}\simeq \frac{\unicode[STIX]{x1D703}^{\ast }}{\unicode[STIX]{x1D719}^{max}\unicode[STIX]{x1D707}_{s}[1-\tan \unicode[STIX]{x1D6FC}/\!\tan \unicode[STIX]{x1D6FC}_{0}]},\end{eqnarray}$$

where (2.19) has been used and it has been assumed that $\bar{\unicode[STIX]{x1D719}}\simeq \unicode[STIX]{x1D719}^{max}$ as the average solid volume fraction $\bar{\unicode[STIX]{x1D719}}$ is only weakly varying in the problem.

The term $\overline{\unicode[STIX]{x1D719}\langle v_{x}\rangle ^{p}}$ has the dimension of a velocity, i.e. the square root of the product of an acceleration and a length scale. In the problem, the only acceleration scale is the buoyant reduced gravity $\tilde{g}=(\unicode[STIX]{x1D70C}^{p}/\unicode[STIX]{x1D70C}^{f}-1)g$ (Duran et al. Reference Duran, Andreotti and Claudin2012), and we have seen that the relevant length scale is the mobile granular layer thickness. Thus, the volume flux should scale as:

(2.22) $$\begin{eqnarray}\overline{\unicode[STIX]{x1D719}\langle v_{x}\rangle ^{p}}\sim \sqrt{(\unicode[STIX]{x1D70C}^{p}/\unicode[STIX]{x1D70C}^{f}-1)g\unicode[STIX]{x1D6FF}_{s}}.\end{eqnarray}$$

Compiling (2.20)–(2.22), the scaling law of the dimensionless sediment transport rate $Q_{s}^{\ast }$ can be written as:

(2.23) $$\begin{eqnarray}Q_{s}^{\ast }=\frac{Q_{s}}{\sqrt{(\unicode[STIX]{x1D70C}^{p}/\unicode[STIX]{x1D70C}^{f}-1)gd^{3}}}\sim \frac{\sqrt{(\unicode[STIX]{x1D70C}^{p}/\unicode[STIX]{x1D70C}^{f}-1)g\unicode[STIX]{x1D6FF}_{s}^{3}}}{\sqrt{(\unicode[STIX]{x1D70C}^{p}/\unicode[STIX]{x1D70C}^{f}-1)gd^{3}}}\sim \left(\frac{\unicode[STIX]{x1D703}^{\ast }}{\unicode[STIX]{x1D719}^{max}\unicode[STIX]{x1D707}_{s}[1-\tan \unicode[STIX]{x1D6FC}/\!\tan \unicode[STIX]{x1D6FC}_{0}]}\right)^{3/2}.\end{eqnarray}$$

This approach not only allows us to recover the scaling law of the dimensionless sediment transport rate with the Shields number and the power $3/2$ in the high Shields number limit, but also allows us to express the influence of the slope through an effective granular friction coefficient, $\unicode[STIX]{x1D707}_{eff}=\unicode[STIX]{x1D707}_{s}[1-\tan \unicode[STIX]{x1D6FC}/\!\tan \unicode[STIX]{x1D6FC}_{0}]$ . The proposed scaling law is consistent with the original definition of the Shields number as a ratio between the driving and resistive forces applied to the granular phase. Indeed, it represents the ratio between the tangential fluid bed shear stress $\unicode[STIX]{x1D70F}_{b}$ and the tangential shear stress associated with the friction between granular layers $\unicode[STIX]{x1D70F}=\unicode[STIX]{x1D707}_{eff}P^{p}=\unicode[STIX]{x1D719}\unicode[STIX]{x1D707}_{eff}(\unicode[STIX]{x1D70C}^{p}-\unicode[STIX]{x1D70C}^{f})gd$ . The dependence of the effective friction coefficient on the distance to the critical debris flow angle, $1-\tan \unicode[STIX]{x1D6FC}/\!\tan \unicode[STIX]{x1D6FC}_{0}$ , is consistent with the analysis of § 2.2, and allows us to take into account both the direct and indirect slope effects on the granular phase in gravity-driven bedload transport. It suggests that the key parameter is not the slope in itself but the ratio between the slope and the critical angle of debris flow. This is of particular importance as the latter is not constant and depends on the maximum solid volume fraction and on the specific density.

3.1 Numerical model

The model is based on a three-dimensional (3-D) discrete element method (DEM) coupled with a 1-D volume-averaged fluid momentum balance. The fine resolution of the granular phase and the two-way coupling ensures an accurate description of the vertical depth profiles and the momentum conservation of the system on average. The model has been described in detail and validated experimentally in Maurin et al. (Reference Maurin, Chauchat, Chareyre and Frey2015), so that only the main characteristics are recalled here. The typical configuration considered here is shown in figure 2.

The bi-periodic 3-D DEM is based on the explicit resolution of the particle phase, solving Newton’s equations of motion for each particle $p$ at position $x_{p}$ :

(3.1) $$\begin{eqnarray}m\frac{\text{d}^{2}\boldsymbol{x}^{p}}{\text{d}t^{2}}=\boldsymbol{f}_{ext}^{p}+\mathop{\sum }_{k\in {\mathcal{N}}}\boldsymbol{f}_{c}^{pk}=\boldsymbol{f}_{g}^{p}+\boldsymbol{f}_{f}^{p}+\mathop{\sum }_{k\in {\mathcal{N}}}\boldsymbol{f}_{c}^{pk},\end{eqnarray}$$

where the sum of the contact forces $\boldsymbol{f}_{c}^{pk}$ is made over the ensemble of nearest neighbours ${\mathcal{N}}$ , $\boldsymbol{f}_{g}^{p}$ is the gravity force, $\boldsymbol{f}_{f}^{p}$ is the force applied by the fluid on particle $p$ . Similarly, the rotation of the particles are solved from Newton’s equations of motion. The contact force between particles is determined from the particles overlap using the classical spring-dashpot contact law (Schwager & Pöschel Reference Schwager and Pöschel2007) which defines a unique normal restitution coefficient $e_{n}$ , and considers tangential friction characterised by a friction coefficient $\unicode[STIX]{x1D707}_{p}$ . The latter two are taken as respectively $e_{n}=0.5$ and $\unicode[STIX]{x1D707}_{p}=0.4$ , as determined from experimental comparisons (Maurin et al. Reference Maurin, Chauchat, Chareyre and Frey2015). The interaction with the fluid phase ( $\boldsymbol{f}_{f}^{p}$ ) is mainly restricted to buoyancy, drag and lift forces in turbulent bedload transport (Nino & Garcia Reference Nino and Garcia1994; Niño & García Reference Niño and García1998). While it is clear that the lift force plays a non-negligible role in turbulent bedload transport (Ji et al. Reference Ji, Munjiza, Avital, Ma and Williams2013), its expression has only been derived in the limits of Stokes flow (Saffman Reference Saffman1965) or an inviscid fluid (Schmeeckle et al. Reference Schmeeckle, Nelson and Shreve2007), and does not apply to the fluid flow regimes associated with turbulent bedload transport (Schmeeckle et al. Reference Schmeeckle, Nelson and Shreve2007). As a consequence, it has been decided to avoid including a controversial expression of the lift force, which appeared to be unnecessary to reproduce accurately turbulent bedload transport experiments (Maurin et al. Reference Maurin, Chauchat, Chareyre and Frey2015). Therefore, the fluid–particle interactions are restricted in the present work to the three-dimensional buoyancy ( $\boldsymbol{f}_{b}^{p}$ ) and drag forces ( $\boldsymbol{f}_{D}^{p}$ ) (Jackson Reference Jackson2000; Maurin et al. Reference Maurin, Chauchat, Chareyre and Frey2015):

(3.2) $$\begin{eqnarray}\displaystyle & \displaystyle \boldsymbol{f}_{b}^{p}=-\frac{\unicode[STIX]{x03C0}d^{3}}{6}\unicode[STIX]{x1D735}P^{f}, & \displaystyle\end{eqnarray}$$

(3.3) $$\begin{eqnarray}\displaystyle & \displaystyle \boldsymbol{f}_{D}^{p}=\frac{1}{2}\unicode[STIX]{x1D70C}^{f}\frac{\unicode[STIX]{x03C0}d^{2}}{4}C_{D}\Vert \langle \boldsymbol{u}\rangle _{\boldsymbol{x}^{\boldsymbol{p}}}^{f}-\boldsymbol{v}^{\boldsymbol{p}}\Vert (\langle \boldsymbol{u}\rangle _{\boldsymbol{x}^{\boldsymbol{p}}}^{f}-\boldsymbol{v}^{\boldsymbol{p}}), & \displaystyle\end{eqnarray}$$

where the average fluid velocity and the fluid pressure are taken at the centre of particle $p$ , and $\boldsymbol{v}^{\boldsymbol{p}}$ represents the velocity of particle $p$ . The drag coefficient $C_{D}$ depends on the particle Reynolds number $Re_{p}=\Vert \langle \boldsymbol{u}\rangle _{\boldsymbol{x}^{\boldsymbol{p}}}^{f}-\boldsymbol{v}^{\boldsymbol{p}}\Vert d/\unicode[STIX]{x1D708}^{f}$ and takes into account hindrance effects (DallaValle Reference DallaValle1948; Richardson & Zaki Reference Richardson and Zaki1954): $C_{D}=(0.4+24.4/Re_{p})(1-\unicode[STIX]{x1D719})^{-3.1}$ . Knowing the position and velocity of the particles, the fluid pressure and velocity fields, Newton’s equations of motion are solved for the ensemble of particles using the open-source DEM code YADE (Šmilauer et al. Reference Šmilauer2015).

The coupled fluid description solves the volume-averaged momentum balance of the fluid phase, corresponding to (2.7) and (2.9):

(3.4) $$\begin{eqnarray}\displaystyle & \displaystyle 0=\frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D70F}_{xz}^{f}}{\unicode[STIX]{x2202}z}+\frac{\unicode[STIX]{x2202}R_{xz}^{f}}{\unicode[STIX]{x2202}z}+\unicode[STIX]{x1D70C}^{f}(1-\unicode[STIX]{x1D719})g\sin \unicode[STIX]{x1D6FC}-n\langle \,{f_{f}^{p}}_{x}\rangle ^{p}, & \displaystyle\end{eqnarray}$$

(3.5) $$\begin{eqnarray}\displaystyle & \displaystyle 0=-\frac{\unicode[STIX]{x2202}P^{f}}{\unicode[STIX]{x2202}z}+\unicode[STIX]{x1D70C}^{f}(1-\unicode[STIX]{x1D719})g\cos \unicode[STIX]{x1D6FC}-n\langle \,{f_{f}^{p}}_{z}\rangle ^{p}. & \displaystyle\end{eqnarray}$$

Due to the steady uniform character of the problem, the fluid velocity field reduces to its streamwise component and depends only on the vertical position: $\langle \boldsymbol{u}\rangle ^{f}=\langle u_{x}\rangle ^{f}(z)\boldsymbol{e}_{\boldsymbol{x}}$ , as sketched in figure 2. Therefore, the buoyancy force (3.2) is restricted to its wall-normal component and (3.5) leads to a hydrostatic fluid pressure distribution, the wall-normal average drag force being negligible. The solution of the streamwise fluid phase momentum balance (3.4) requires closures for the viscous and Reynolds stress tensors, as well as the determination of the solid volume fraction $\unicode[STIX]{x1D719}$ and the momentum transfer associated with the granular phase interactions. Considering a Newtonian fluid, the effective viscous shear stress is expressed as:

(3.6) $$\begin{eqnarray}\unicode[STIX]{x1D70F}_{xz}^{f}=\unicode[STIX]{x1D70C}^{f}(1-\unicode[STIX]{x1D719})\unicode[STIX]{x1D708}^{f}\frac{\text{d}\langle u_{x}\rangle ^{f}}{\text{d}z},\end{eqnarray}$$

with $\unicode[STIX]{x1D708}^{f}$ the clear fluid kinematic viscosity and $\langle u_{x}\rangle ^{f}$ the volume-averaged streamwise fluid phase velocity. The Reynolds shear stress is based on the eddy viscosity concept ( $\unicode[STIX]{x1D708}^{t}$ ) using a mixing length formulation:

(3.7) $$\begin{eqnarray}R_{xz}^{f}=\unicode[STIX]{x1D70C}^{f}\unicode[STIX]{x1D708}^{t}\frac{\text{d}\langle u_{x}\rangle ^{f}}{\text{d}z}\quad \text{with }\unicode[STIX]{x1D708}^{t}=(1-\unicode[STIX]{x1D719})l_{m}^{2}\left|\frac{\text{d}\langle u_{x}\rangle ^{f}}{\text{d}z}\right|,\end{eqnarray}$$

where the mixing length is taken similarly to Li & Sawamoto (Reference Li and Sawamoto1995) as:

(3.8) $$\begin{eqnarray}l_{m}(z)=\unicode[STIX]{x1D705}\int _{0}^{z}\frac{\unicode[STIX]{x1D719}^{max}-\unicode[STIX]{x1D719}(\unicode[STIX]{x1D701})}{\unicode[STIX]{x1D719}^{max}}\,\text{d}\unicode[STIX]{x1D701},\end{eqnarray}$$

with $\unicode[STIX]{x1D705}=0.41$ , the von Kármán constant. The formulation adopted allows us to recover the law of the wall (Prandtl Reference Prandtl1926) in a clear fluid, while completely damping the turbulence inside the granular bed at maximum packing fraction ( $\unicode[STIX]{x1D719}^{max}$ ).

The solid volume fraction ( $\unicode[STIX]{x1D719}$ ) and the fluid–particle interaction term ( $n\langle \,f_{f}^{p}\rangle ^{p}$ ) are determined from spatial averaging of the discrete solid phase. Considering cubic boxes of finite wall-normal length scale $l_{z}$ , these two terms are averaged over the whole length and width of the granular bi-periodic cell. In order to solve the important wall-normal gradients present in turbulent bedload transport, a small wall-normal weighting function length scale has been adopted (typically $l_{z}=d/30$ ), and this choice has been confirmed by the experimental validation (Maurin et al. Reference Maurin, Chauchat, Chareyre and Frey2015).

The 3-D DEM and the fluid model are solved as transient problems applying a fixed bottom boundary condition for both the fluid ( $\langle \boldsymbol{u}\rangle ^{f}(z=0)=0$ ) and the particle phase (fixed random particles) and imposing the position of the water free surface ( $d\langle u_{x}\rangle ^{f}/\text{d}z(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 (Maurin Reference Maurin2015; Maurin et al. Reference Maurin, Chauchat, Chareyre and Frey2015). 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 (Maurin et al. Reference Maurin, Chauchat, Chareyre and Frey2015): $\unicode[STIX]{x0394}t_{f}=10^{-2}s$ with respect to $\unicode[STIX]{x0394}t_{p}\sim O(10^{-4}{-}10^{-5})~\text{s}$ .

Table 1.

Parameters of the simulations studied and symbol correspondence. Dimensionless number definitions are given in appendix B.

The model has been compared with experiments and has shown its capability to describe accurately both the sediment transport rate and the granular depth structure in turbulent bedload transport (Maurin et al. Reference Maurin, Chauchat, Chareyre and Frey2015).