2. Governing equations and semi-analytical solutions

The saltating particles are assumed to be identical spheres of diameter d and mass density ρ^s. A shearing flow of a fluid of mass density ρ^f and molecular viscosity $\mu$^f drives the flow in the presence of gravity, with g the gravitational acceleration. The mean horizontal velocities of the particles and the fluid are u and U, respectively. We assume that the flow is steady and uniform, so that the velocities are only functions of the vertical distance from the bed y. We imagine that the bed is made bumpy by gluing identical particles of diameter d_w, in close contact to each other, over a flat plate. Hence, the ratio of d_w to d is a natural measure of the bed roughness. We characterize the particles through the fall particle Reynolds number $R = {\rho ^f}\sqrt {g(r - 1)/r} {d^{3/2}}/{\mu ^f}$, where r = ρ^s/ρ^f is the density ratio. A sketch of the flow configuration is shown in figure 1. In what follows, all quantities are made dimensionless using the diameter and mass density of the particles and the reduced gravitational acceleration, $g(r - 1)/r$. Then, lengths, velocities and stresses are expressed in units of d, ${[g(r - 1)d/r]^{1/2}}$ and ${\rho^s}g(r - 1)d/r$, respectively.

Figure 1. Sketch of a particle driven into saltation over a rigid, bumpy bed by a shearing flow.

The balance equations for the y and x-momenta of the particles, under steady and fully developed conditions, are, respectively,

(2.1)\begin{equation}\frac{{\textrm{d}{\sigma _y}}}{{\textrm{d}y}} = - c;\end{equation}

and

(2.2)\begin{equation}\frac{{\textrm{d}s}}{{\textrm{d}y}} = - c{C_D}(U - u),\end{equation}

where ${\sigma _y}$ and s are the particle normal stress in the y-direction and the particle shear stress, respectively; c is the particle volume concentration; and C_D is the drag coefficient that, for a single particle, has a component independent of the relative velocity between the single particle and the fluid and a component proportional to the absolute value of that velocity difference. Here, we assume perhaps the simplest form of C_D (Dallavalle Reference Dallavalle1943), which we make independent of y by evaluating it at the bed in the case of saltation in viscous flows,

(2.3a)\begin{equation}{C_D} = \frac{{18}}{{St}} + \frac{{0.3}}{r}|{U_0} - {u_0}|,\end{equation}

where St = rR is the fall Stokes number and the subscript 0 indicates quantities evaluated at the bed. We employ the terms Stokes drag to refer to situations in which ${C_D} \simeq 18/St$, form drag for situations in which ${C_D} \simeq 0.3|{U_0} - {u_0}|/r$ and nonlinear drag for the generic case in which both Stokes and form drag are present. A more appropriate expression for the average C_D should also involve the strength of the particle velocity fluctuations (Jenkins & Hanes Reference Jenkins and Hanes1998). However, including it would have only a small quantitative effect on the results. Also permitting the drag coefficient to vary along y does not significantly alter the solution to the continuum model. We anticipate that, while (2.3a) works well for saltation in viscous shearing flows over rigid, bumpy beds, where the fluid and particle velocity profiles are approximately linear, a more appropriate expression for saltation in turbulent flows should incorporate the difference in the concentration-weighted average velocities between the two phases, as proposed by Pähtz et al. (Reference Pähtz, Liu, Xia, Hu, He and Tholen2021)

(2.3b)\begin{equation}{C_D} = \frac{{18}}{{St}} + \frac{{0.3}}{r}|\bar{U} - \bar{u}|,\end{equation}

where $\bar{U} - \bar{u} = \int_0^h {c(U - u)\,\textrm{d}y} /\int_0^h {c\,\textrm{d}y}$, with h the depth of the saltation layer.

For flows over rigid beds, we take the no-slip boundary condition for the fluid

(2.4)\begin{equation}{U_0} = 0.\end{equation}

If saltating particles with a large value of the coefficient of sliding friction impact a rigid, bumpy bed at a small mean angle θ with respect to the horizontal (figure 1), then, using the analysis of Lämmel et al. (Reference Lämmel, Dzikowski, Kroy, Oger and Valance2017) in Appendix C,

(2.5)\begin{equation}{u_0} = {\alpha _u}\frac{1}{{{C_D}\theta }},\end{equation}

where α_u is a strictly positive coefficient of order unity that is weakly dependent on the rebound properties of the particles and the flow regime of the fluid. The reason for the latter dependence is that the probability distributions of impact angles are qualitatively different for saltation in viscous and turbulent shearing flows (Appendix C).

In a steady, uniform flow, in the absence of a horizontal pressure gradient, the sum of the fluid and particle shear stress, S and s, respectively, is constant and equal to the far-field fluid shear stress, which, in dimensionless terms, is the Shields parameter

(2.6)\begin{equation}s + S = Sh.\end{equation}

To close the problem, we require constitutive relations for the particle and fluid stresses.

In the absence of particle collisions above the bed, the only mechanism responsible for the particle stresses is the transfer of momentum associated with the particles crossing a reference surface. Hence, the particle normal stress in the y-direction is simply given by the average vertical flux of y-momentum

(2.7)\begin{equation}{\sigma _y} = c{T_y},\end{equation}

where ${T_y}$ is the mean square of the vertical velocity fluctuations of the particles . There is a corresponding mean square of the horizontal velocity fluctuations T_x, where T_x and T_y are the diagonal components of the second moment of the particle velocity tensor.

The distribution of the additional hydrodynamic field, T_y, along y is governed by the balance of kinetic energy associated with the particle vertical motion, that is, the yy-component of the particle second-moment tensor, which, in a steady and fully developed flow, and in the absence of mid-trajectory collisions and fluid velocity fluctuations, reduces to (Saha & Alam Reference Saha and Alam2016, Reference Saha and Alam2017)

(2.8)\begin{equation}{-}\frac{{\textrm{d}{Q_{yyy}}}}{{\textrm{d}y}} - 2{C_D}{\sigma _y} = 0.\end{equation}

Here, ${Q_{yyy}}$ is the y-component of the flux of the kinetic energy associated with the vertical particle velocity fluctuations and $2{C_D}{\sigma _y}$ is its dissipation due to the fluid drag.

Assuming that the fluid motion is only horizontal and that the drag coefficient is independent of y, we can obtain approximate analytical expressions for the trajectories of the saltating particles even in the turbulent case (Appendix A). In doing so, we distinguish between ascending and descending particles. Then, the horizontal, $\xi _x^ -$, and vertical, $\xi _y^ -$, velocities of any descending particle at a given distance y from the bed can be obtained analytically from the horizontal, $\xi _x^ +$, and vertical, $\xi _y^ +$, velocities of the same particle at the same distance y from the bed during its ascending motion. In the classical framework of statistical mechanics, we introduce a velocity distribution function, ${f^ + }$, for the ascending particles, so that $({\rm \pi}/6)\int_{all\;{\xi ^ + }} {{f^ + }\,{\textrm{d}^2}{\boldsymbol{\xi} ^ + }}$ gives the concentration of the ascending particles. With this, ${Q_{yyy}} = ({\rm \pi}/6)\int_{all\;{\xi ^ + }} {(\xi _y^{+2} - \xi _y^{-2}){\xi _y^{+}} {f^ + }\,{\textrm{d}^2}{\boldsymbol{\xi} ^ + }}$. With the further assumptions that the velocity distribution function is an anisotropic Maxwellian (Creyssels et al. Reference Creyssels, Dupont, Ould El Moctar, Valance, Cantat, Jenkins, Pasini and Rasmussen2009) and that the vertical velocity $\xi _y^ +$ of the ascending particles is much larger than 1/C_D (the settling velocity), we derive a simple expression for ${Q_{yyy}}$ in Appendix B

(2.9)\begin{equation}{Q_{yyy}} = \sqrt {\frac{8}{\rm \pi}} {c^ + }{(T_y^ + )^{3/2}}.\end{equation}

Here, c ⁺ and $T_y^ + $ are the volume concentration and the mean square of the vertical velocity fluctuations of the ascending particles, respectively. In a steady state, c ⁺ and $T_y^ + $ are related to c and T_y through (see Appendix B)

(2.10)\begin{equation}c = {c^ + }\left( {2 + \frac{3}{4}{C_D}\sqrt {\frac{{2T_y^ + }}{\rm \pi}} } \right);\end{equation}

and

(2.11)\begin{equation}{T_y} = \frac{{4\sqrt{{\rm \pi} } {C_D}T_y^{+} + 4\sqrt {2T_y^ + } }}{{8\sqrt{{\rm \pi} } {C_D} + 3C_D^2\sqrt {2T_y^ + } }}.\end{equation}

As explicitly shown in Appendix B, the assumption that the vertical velocity of the ascending particles is much larger than the settling velocity is a good approximation for ${C_D}\sqrt {{T_y}} > 1$. We will see in § 4 that this condition is not satisfied only near the top of the saltation layer.

Using (2.7), (2.10) and (2.11) into (2.1) leads to

(2.12)\begin{equation}\frac{{\textrm{d}{c^ + }}}{{\textrm{d}y}} = - {c^ + }\frac{{8{C_D} + 3C_D^2\sqrt {2T_y^ +{/}{\rm \pi}} }}{{4{C_D}T_y^{+} + 4\sqrt {2T_y^ +{/}{\rm \pi}} }} - {c^ + }\frac{{2{C_D} + \sqrt {2/{\rm \pi} T_y^ + } }}{{2{C_D}T_y^{+} + 2\sqrt {2T_y^ +{/}{\rm \pi}} }}\frac{{\textrm{d}T_y^ + }}{{\textrm{d}y}},\end{equation}

that, combined with (2.8) and (2.9), permits an ordinary differential equation for $T_y^ + $ to be obtained

(2.13)\begin{equation}[\sqrt {8{\rm \pi}} {C_D}{(T_y^ + )^{3/2}} + 8T_y^ + ]\frac{{\textrm{d}T_y^ + }}{{\textrm{d}y}} = (6 - 4{\rm \pi} )C_D^2{(T_y^ + )^2} - 8T_y^ + .\end{equation}

The analytical solution of this, with the boundary condition $T_y^ + (y = 0) = T_{y,0}^ +$, for the vertical particle velocity fluctuations at the rigid boundary, is

(2.14)\begin{align} & 4\sqrt {2{{\rm \pi} (}2{\rm \pi} - 3\textrm{)}} {\tan ^{ - 1}}\left( {\dfrac{{\sqrt {2{\rm \pi} - 3} }}{2}{C_D}\sqrt {T_y^ + } } \right)\nonumber\\ & \qquad - 2(2{\rm \pi} - 3)\left\{ {\sqrt {2{\rm \pi}} {C_D}\sqrt {T_y^ + } + 2\log [4 + (2{\rm \pi} - 3)C_D^2T_y^ + ]} \right\}\nonumber\\ & \quad = {(3 - 2{\rm \pi} )^2}C_D^2y + 4\sqrt {2{{\rm \pi} (}2{\rm \pi} - 3\textrm{)}} {\tan ^{ - 1}}\left( {\dfrac{{\sqrt {2{\rm \pi} - 3} }}{2}{C_D}\sqrt {T_{y,0}^ + } } \right)\nonumber\\ & \qquad - 2(2{\rm \pi} - 3)\left\{ {\sqrt {2{\rm \pi}} {C_D}\sqrt {T_{y,0}^ + } + 2\log [4 + (2{\rm \pi} - 3)C_D^2T_{y,0}^ + ]} \right\}. \end{align}

The non-zero value of $T_{y,0}^ +$ is associated with the rebound velocity of the particles at a rigid, bumpy bed.

Under the same assumptions that were employed to derive equation (2.5), a simple dependence of $T_{y,0}^ +$ calculated at the bed on the impact angle θ is obtained in Appendix C

(2.15)\begin{equation}T_{y,0}^ + = {\alpha _T}\frac{1}{{C_D^2\theta }}\frac{{2{d_w}}}{{1 + {d_w}}},\end{equation}

where α_T is another strictly positive coefficient of order unity that is weakly dependent on the rebound properties of the particles and the flow regime of the fluid. We emphasize that the impact angle, and consequently $T_{y,0}^ +$, remains an unknown at this stage of the analysis. We will determine it in the next two sub-sections, when we describe quantities associated with the horizontal motion of the particles that are influenced by the flow regime of the fluid.

Equation (2.14) provides the analytical distribution of $T_y^ +$ along y, which, when employed in (2.11), gives the analytical distribution of T_y shown in figure 2(a) for different values of ${C_D}\sqrt {T_{y,0}^ + }$, corresponding to an impact angle θ in (2.15) in the range 1° to 20°. As shown in figure 2(b), the distribution of T_y is roughly linear and well approximated as

(2.16)\begin{equation}{T_y} = \frac{{h - y}}{h}{T_{y,0}},\end{equation}

where the depth of the saltation layer, h, is to be determined. The mean square of the particle velocity fluctuations at the bed, ${T_{y,0}}$, can be obtained from (2.11), when $T_{y,0}^ +$ is known.

Figure 2. (a) Normalized profiles of $T_y^ + $ (dot-dashed lines) and T_y (dashed lines) obtained from (2.14) and (2.11), respectively, when $C_D^2T_{y,0}^ + = 2.5$ (blue lines), $C_D^2T_{y,0}^ + = 10$ (orange lines) and $C_D^2T_{y,0}^ + = 100$ (purple lines). (b) Normalized profiles of T_y (dashed lines, same colour legend of figure 2a) and the linear distribution of (2.16) (solid black line).

Inserting (2.7) into (2.1), with (2.16), and integrating gives the following power-law distribution of the particle concentration:

(2.17)\begin{equation}c = {c_0}{\left( {\frac{{h - y}}{h}} \right)^{(h - {T_{y,0}})/{T_{y,0}}}},\end{equation}

where c ₀ is the particle concentration at the bed. Integrating equation (2.1) with (2.17) also gives

(2.18)\begin{equation}{\sigma _y} = {c_0}{T_{y,0}}{\left( {\frac{{h - y}}{h}} \right)^{h/{T_{y,0}}}}.\end{equation}

Upon introducing the hold-up, $M = \int_0^h {c\,\textrm{d}y}$, i.e. the particle mass per unit basal area, by integrating the concentration through the saltation layer, with σ_y = 0 at y = h, we obtain

(2.19)\begin{equation}{c_0} = \frac{M}{{{T_{y,0}}}}.\end{equation}

Integrating the energy balance, (2.8) with (2.18) and (2.19), provides

(2.20)\begin{equation}{Q_{yyy}} = 2{C_D}M{T_{y,0}}\frac{h}{{h + {T_{y,0}}}}\left[ {{{\left( {\frac{{h - y}}{h}} \right)}^{(h + {T_{y,0}})/{T_{y,0}}}} - 1} \right] + {Q_{yyy,0}},\end{equation}

where the value of Q_yyy at the bed can be obtained from the ratio of (2.9) and (2.7), with (2.10), (2.11) and (2.19), as

(2.21)\begin{equation}{Q_{yyy,0}} = \frac{{4{C_D}T_{y,0}^ + }}{{{C_D}\sqrt {2{\rm \pi} T_{y,0}^ + } + 2}}M.\end{equation}

Given that Q_yyy must vanish at the top of the saltation layer, the depth h of the saltation layer is determined from (2.20) as

(2.22)\begin{equation}h = \frac{{{Q_{yyy,0}}}}{{2{C_D}M{T_{y,0}} - {Q_{yyy,0}}}}{T_{y,0}}.\end{equation}

The governing equations, constitutive relations and analytical results described so far apply to saltation in both viscous and in turbulent shearing flows and can be calculated only after the determination of the impact angle θ, required in (2.5) and (2.15). To proceed, we must distinguish between the two regimes of the fluid shearing flow.

Before doing that, we support the assumption that particle collisions above the bed can be neglected. As suggested by Pasini & Jenkins (Reference Pasini and Jenkins2005), chances of mid-trajectory collisions are low if the mean free path, λ, the average distance travelled by a particle in between two successive collisions predicted by the kinetic theory of granular gases, is less than twice the height of the particle trajectory, here determined by gravity and fluid drag. If we use the expression for the mean free path in a dilute gas of Chapman & Cowling (Reference Chapman and Cowling1970), evaluate this at the bed and take h as the average height of the particle trajectories, we obtain

(2.23)\begin{equation}{\lambda _0} = \frac{1}{{6{c_0}\sqrt 2 }} \ge 2h.\end{equation}

Equation (2.23), used with (2.22), implies that, for given values of r, St and Sh, there is a maximum hold-up above which mid-trajectory collisions play a role. That is, because the mean free path is a decreasing function of the particle concentration, as the hold-up increases, collisions become more likely. However, there is also a limiting hold-up at which particles begin to deposit on the bed. Analysis of periodic trajectories over rigid beds (Jenkins & Valance Reference Jenkins and Valance2014) and discrete simulations of the present paper indicate that the concentration at the bed at which deposition begins is greater than that at which collisions begin to occur. Consequently, in saltation over such beds, collisions occur before deposition.

2.1. Saltation in viscous shearing flows

In the absence of turbulence, the expression for the fluid shear stress is (Valance & Berzi Reference Valance and Berzi2022)

(2.24)\begin{equation}S = \frac{{1 - c}}{{St}}\frac{{\textrm{d}U}}{{\textrm{d}y}} \simeq \frac{1}{{St}}\frac{{\textrm{d}U}}{{\textrm{d}y}},\end{equation}

where we have neglected the particle concentration with respect to unity. Equation (2.24) implies that, in the absence of particles and in the boundary layer approximation that we have employed, the fluid velocity distribution would be linear.

Using statistical mechanics arguments, the particle shear stress can be obtained as $s = - ({\rm \pi}/6)\int_{all\;{\xi ^ + }} {({\xi _x^{+}} - {\xi _x^{-}} ){\xi _y^{+}} {f^ + }\,{\textrm{d}^2}{\boldsymbol{\xi} ^ + }}$. To determine $\xi _x^ -$, we assume in Appendix A that the fluid velocity profile encountered by an ascending particle during its ballistic trajectory above a certain position y is locally linear; that is, it is linear in the region between y and the top of the trajectory. However, the slope of the linear profile changes with y due to the drag of the particles. This assumption and the assumptions on the form of the velocity distribution function ${f^ + }$ and ${C_D}\xi _y^ + \gg 1$, that we have already employed in the derivation of the expression for Q_yyy in (2.9), permit the derivation in Appendix B of a simple expression for the particle shear stress

(2.25)\begin{equation}s = \frac{3}{5}{C_D}{\sigma _y}\left( {U - u + \frac{2}{3}\frac{1}{{C_D^2}}\frac{{\textrm{d}U}}{{\textrm{d}y}}} \right).\end{equation}

Interestingly, the particle shear stress does not depend on the particle shear rate, as in Jenkins et al. (Reference Jenkins, Cantat and Valance2010), but only on the velocity difference and on the fluid shear rate. The physical reason is that the particles do not interact with each other, but only with the surrounding fluid.

As explained in Appendix C, for saltation in viscous flows the mean angle θ between the particle velocity and the horizontal before the impact with the bed is estimated in the limit of large vertical velocity of ascending particles in Appendix A as

(2.26)\begin{equation}\theta = \frac{{{C_D}}}{{St(Sh - {s_0})}},\end{equation}

where s ₀ is the particle shear stress at the bed.

Using (2.4) and (2.5) in (2.3a), with (2.26), we obtain a relationship between the drag coefficient and the particle shear stress at the bed

(2.27)\begin{equation}C_D^3 = \frac{{18}}{{St}}C_D^2 + \frac{{0.3}}{r}{\alpha _u}St(Sh - {s_0}).\end{equation}

The particle shear stress at the bed is obtained from (2.25), with (2.4)–(2.6), (2.19) and (2.24), as

(2.28)\begin{equation}{s_0} = \frac{{(2 - 3{\alpha _u})MStSh}}{{5{C_D} + (2 - 3{\alpha _u})MSt}}.\end{equation}

Equation (2.28) implies that for the particle shear stress at the bed to be positive, α_u must be less than 2/3. The system of (2.27) and (2.28) can be solved to determine the drag coefficient and the particle shear stress at the bed once the values of d_w, r, St, Sh and M are given. Then, the impact angle (2.26), the values of all the variables at the bed ((2.5), (2.15) and (2.21)), the depth of the saltation layer (2.22), and the analytical distributions of $T_y^ + $, T_y, c, σ_y and Q_yyy ((2.14), (2.16)–(2.18) and (2.20)) can also be calculated for saltation in viscous shearing flows.

Figure 3 shows the variation of the impact angle, the drag coefficient and the particle shear stress at the bed with the particle hold-up, as predicted by ((2.26)–(2.28)), for, e.g. saltation of 100 ${\rm \mu}$m basalt grains in viscous shearing flows on Venus, assuming that α_u = 0.60 (§ 4 provides more details on the choice of the parameters). Equation (2.23) is satisfied for the range of hold-up in figure 3. Notice that for almost the entire range of hold-up, the particle shear stress at the bed increases linearly with M, while θ and C_D are constant.

Figure 3. Predicted impact angle in degrees (solid line, –), drag coefficient (dashed line, ‐‐) and particle shear stress at the bed (dot-dashed line, ⋅-) as functions of the hold-up for saltation in viscous shearing flows when d_w = 1, St = 100, r = 50, and Sh = 0.05 (with α_u = 0.60).

Using (2.25) in (2.2), with (2.6), (2.18) and (2.24), gives a first-order, linear, non-homogeneous, differential equation for the particle shear stress

(2.29) \begin{align}\frac{{\textrm{d}s}}{{\textrm{d}y}} = - \left[ {\frac{5}{3}\frac{1}{{{T_{y,0}}}}{{\left( {\frac{{h - y}}{h}} \right)}^{ - 1}} + \frac{2}{3}\frac{{{c_0}St}}{{{C_D}}}{{\left( {\frac{{h - y}}{h}} \right)}^{(h - {T_{y,0}})/{T_{y,0}}}}} \right]s + \frac{2}{3}\frac{{{c_0}StSh}}{{{C_D}}}{\left( {\frac{{h - y}}{h}} \right)^{(h - {T_{y,0}})/{T_{y,0}}}}.\end{align}

The analytical solution of (2.29) is

(2.30) \begin{align} s & = Sh{\left({\dfrac{2}{3}\dfrac{{MSt}}{{{C_D}}}} \right)^{5/3}}\varGamma \left[ { - \dfrac{2}{3},\dfrac{2}{3}\dfrac{{MSt}}{{{C_D}}}{{\left( {\dfrac{{h - y}}{h}} \right)}^{h/{T_{y,0}}}}} \right]\nonumber\\ &\quad \times \exp \left\{ {\ln \left[ {{{\left( {\dfrac{{h - y}}{h}} \right)}^{(5/3)(h/{T_{y,0}})}}} \right] + \dfrac{2}{3}\dfrac{{MSt}}{{{C_D}}}{{\left( {\dfrac{{h - y}}{h}} \right)}^{h/{T_{y,0}}}}} \right\}\nonumber\\ & \quad - Sh{\left( {\dfrac{2}{3}\dfrac{{MSt}}{{{C_D}}}} \right)^{5/3}}\varGamma \left[ { - \dfrac{2}{3},\dfrac{2}{3}\dfrac{{MSt}}{{{C_D}}}} \right]\nonumber\\ &\quad \times \exp \left\{ {\ln \left[ {{{\left( {\dfrac{{h - y}}{h}} \right)}^{(5/3)(h/{T_{y,0}})}}} \right] + \dfrac{2}{3}\dfrac{{MSt}}{{{C_D}}}{{\left( {\dfrac{{h - y}}{h}} \right)}^{h/{T_{y,0}}}}} \right\}\nonumber\\ & \quad + {s_0} \exp \left\{ {\ln \left[ {{{\left( {\dfrac{{h - y}}{h}} \right)}^{(5/3)(h/{T_{y,0}})}}} \right] + \dfrac{2}{3}\dfrac{{MSt}}{{{C_D}}}{{\left( {\dfrac{{h - y}}{h}} \right)}^{h/{T_{y,0}}}} - \dfrac{2}{3}\dfrac{{MSt}}{{{C_D}}}} \right\}, \end{align}

in which $\varGamma (s,x) \equiv \int_x^\infty {{t^{s - 1}}\,{e^{ - t}}\,{\rm d}t}$ is the upper incomplete gamma function.

The fluid velocity profile is determined by integrating the constitutive relation (2.24) $\textrm{d}U/\textrm{d}y = St(Sh - s)$, with the particle shear stress given by (2.30). Unfortunately, there is no general analytical solution. Hence, we obtain the distribution of U by numerically solving this ordinary differential equation, with the no-slip condition, (2.4), at the bed. Once U is determined, the particle horizontal velocity is given by (2.25) as

(2.31)\begin{equation}u = U - \frac{5}{3}\frac{s}{{{C_D}{\sigma _y}}} + \frac{2}{3}\frac{{St(Sh - s)}}{{C_D^2}}.\end{equation}

We can now calculate the particle flux q per unit basal area through numerical integration as

(2.32)\begin{equation}q = \int_0^h {cu\,\textrm{d}y} .\end{equation}

In the special case of rarefied saltation, that is for $M \to 0$ and $Sh \gg s$, the analytical solution of (2.29) simplifies to

(2.33)\begin{equation}s = {s_0}{\left( {\frac{{h - y}}{h}} \right)^{(5/3)(h/{T_{y,0}})}} + \frac{{MStSh}}{{{C_D}}}\left[ {{{\left( {\frac{{h - y}}{h}} \right)}^{h/{T_{y,0}}}} - {{\left( {\frac{{h - y}}{h}} \right)}^{(5/3)(h/{T_{y,0}})}}} \right],\end{equation}

and the fluid velocity profile is simply

(2.34)\begin{equation}U = (StSh)y.\end{equation}

Then, from (2.18), (2.28), (2.31), (2.33) and (2.34), the particle horizontal velocity is

(2.35)\begin{equation}u = \left[ {(1 + {\alpha_u}){{\left( {\frac{{h - y}}{h}} \right)}^{(2/3)(h/{T_{y,0}})}} + C_D^2y - 1} \right]\frac{{StSh}}{{C_D^2}}.\end{equation}

Notice that, in the rarefied limit of saltation in viscous shearing flows, the particle and the fluid velocity profiles are independent of the hold-up. Using (2.17), (2.19) and (2.35) in (2.32), and integrating, gives the particle flux as

(2.36)\begin{equation}q = MStSh{T_{y,0}}\frac{h}{{h + {T_{y,0}}}}.\end{equation}

The above equations permit the determination of profiles of particle concentration, particle and fluid velocities and particle and fluid stresses in viscous shearing flow that will be compared with the results of DC numerical simulations in § 4.

2.2. Saltation in turbulent shearing flows

In the case of saltation in turbulent shearing flows, we model the fluid shear stress using the classical mixing length approach (Jenkins et al. Reference Jenkins, Cantat and Valance2010)

(2.37)\begin{equation}S = \frac{{1 - c}}{r}{\kappa ^2}{(y + {y_0})^2}{\left( {\frac{{\textrm{d}U}}{{\textrm{d}y}}} \right)^2} \simeq \frac{{{\kappa ^2}{{(y + {y_0})}^2}}}{r}{\left( {\frac{{\textrm{d}U}}{{\textrm{d}y}}} \right)^2},\end{equation}

where κ = 0.41 is von Kármán's constant and ${y_0} = \sqrt {r/Sh} /(9St) + ({k_s}/30)[1 - \exp ( - {k_s}St\sqrt {Sh/r} /26)]$ is an empirical expression for the origin of the logarithmic fluid velocity profile that encompasses hydrodynamically smooth, rough and transitional beds (Guo & Julien Reference Guo and Julien2007). The roughness length scale k_s is taken to be equal to d_w (Ho Reference Ho2012).

Equation (2.37) implies that, in the absence of particles and in the boundary layer approximation that we have employed, the fluid velocity distribution would be logarithmic. We assume, as in Pähtz et al. (Reference Pähtz, Liu, Xia, Hu, He and Tholen2021), that the fluid velocity profile encountered by an ascending particle during its ballistic trajectory is logarithmic even if other particles are present. This assumption, and the assumptions that we have already employed in the derivation of the expression for Q_yyy, (2.9), permits the derivation in Appendix B of an expression for the particle shear stress, $s = - ({\rm \pi}/6)\int_{all\;{\xi ^ + }} {(\xi _x^ + - \xi _x^ - )\xi _y^ + {f^ + }\,{\textrm{d}^2}{\xi ^ + }}$

(2.38) \begin{align} s &= \frac{3}{5}{C_D}{\sigma _y}\left\{ {U - u - \frac{{\sqrt {rS} }}{\kappa }\exp \left[ {{C_D^2}\frac{{\sqrt {rS} }}{\kappa }{{\left( {\frac{{\textrm{d}U}}{{\textrm{d}y}}} \right)}^{ - 1}}} \right]Ei\left[ { - {C_D^2}\frac{{\sqrt {rS} }}{\kappa }{{\left( {\frac{{\textrm{d}U}}{{\textrm{d}y}}} \right)}^{ - 1}}} \right]} \right\},\end{align}

in which $Ei(x) \equiv{-} \int_{ - x}^\infty {{e^{ - t}}{t^{ - 1}}\,{\rm d}t} $ is the exponential integral. As in the viscous case of (2.25), the particle shear stress in turbulent shearing flows does not depend on the particle shear rate. The saltating particles, when treated as a continuous medium, can experience shear stress even if the horizontal particle velocity is uniform across the flow domain.

As explained in Appendix C, the mean impact angle θ for saltation in turbulent flows in the limit of large vertical velocity of ascending particles is estimated as

(2.39)\begin{equation}\theta \simeq{-} {\left[ {{C_D}\frac{{\sqrt {r(Sh - {s_0})} }}{\kappa }\exp (C_D^2{y_0})Ei( - C_D^2{y_0})} \right]^{ - 1}}.\end{equation}

Integrating equation (2.2) between y = 0 and y = h, with vanishing particle shear stress at the top of the saltation layer gives ${s_0} = M{C_D}(\bar{U} - \bar{u})$ (Pähtz et al. Reference Pähtz, Liu, Xia, Hu, He and Tholen2021). Using this in (2.3b) gives the drag coefficient in terms of the particle shear stress at the bed

(2.40)\begin{equation}{C_D} = \frac{9}{{St}} + \sqrt {\frac{{81}}{{S{t^2}}} + \frac{{0.3}}{r}\frac{{{s_0}}}{M}} .\end{equation}

Equation (2.38) evaluated at the bed, with (2.5) and (2.39) gives the particle shear stress at the bed as a function of the impact angle

(2.41)\begin{equation}{s_0} = \frac{3}{5}\frac{M}{\theta }(1 - {\alpha _u}).\end{equation}

Equation (2.41) implies that, for the particle shear stress at the bed to be positive, α_u must be less than 1. Then, (2.39) with (2.40) and (2.41) results in an implicit equation for θ

(2.42) \begin{align} \dfrac{1}{\theta }{\left({\dfrac{9}{{St}} + \sqrt {\dfrac{{81}}{{S{t^2}}} + \dfrac{{0.9}}{{5r}}\dfrac{{1 - {\alpha_u}}}{\theta }} } \right)^{ - 1}} & \simeq{-} \dfrac{1}{\kappa }\sqrt {r\left( {Sh - \dfrac{{3M}}{5}\dfrac{{1 - {\alpha_u}}}{\theta }} \right)}\nonumber\\ &\quad \times \exp \left[ {{{\left( {\dfrac{9}{{St}} + \sqrt {\dfrac{{81}}{{S{t^2}}} + \dfrac{{0.9}}{{5r}}\dfrac{{1 - {\alpha_u}}}{\theta }} } \right)}^2}{y_0}} \right]\nonumber\\ & \quad \times Ei\left[ { - {{\left( {\dfrac{9}{{St}} + \sqrt {\dfrac{{81}}{{S{t^2}}} + \dfrac{{0.9}}{{5r}}\dfrac{{1 - {\alpha_u}}}{\theta }} } \right)}^2}{y_0}} \right].\end{align}

This can be solved to determine the impact angle. Once θ is known, s ₀ (2.41), C_D (2.40), the boundary values of the remaining variables at the bed ((2.5), (2.15) and (2.21)), the depth of the saltation layer (2.22), and the analytical distributions of $T_y^ + $, T_y, c, σ_y and Q_yyy ((2.14), (2.16)–(2.18) and (2.20)) can also be calculated in saltation in turbulent shearing flows.

Figure 4 shows the variation of the impact angle, the drag coefficient and the particle shear stress at the bed with the particle hold-up, as predicted by ((2.40)–(2.42)), for saltation of 240 ${\rm \mu}$m sand grains in a turbulent wind on Earth, with α_u = 0.85 (see § 4 for more details about the choice of the parameters). The range of the hold-ups in figure 4 is that for which (2.23) is satisfied. As for saltation in viscous shearing flows (figure 3), the particle shear stress at the bed increases linearly with M, while θ and C_D are almost constant.

Figure 4. Predicted impact angle in degrees (solid line, –), drag coefficient (dashed line, ‐‐) and particle shear stress at the bed (dot-dashed line, ⋅-) as functions of the hold-up for saltation in turbulent shearing flows with d_w = 1.5, St = 1681, r = 2208 and Sh = 0.04 (with α_u = 0.85).

Equation (2.2), with (2.38) and (2.6), and (2.37) can be written as a system of two ordinary differential equations

(2.43)\begin{equation}\frac{{\textrm{d}s}}{{\textrm{d}y}} = - {C_D}c\frac{{\sqrt {r(Sh - s)} }}{\kappa }\exp [C_D^2(y + {y_0})]Ei[ - C_D^2(y + {y_0})] - \frac{5}{3}\frac{s}{{{T_y}}};\end{equation}

and

(2.44)\begin{equation}\frac{{\textrm{d}U}}{{\textrm{d}y}} = \frac{{\sqrt {r(Sh - s)} }}{{\kappa (y + {y_0})}};\end{equation}

these can be numerically integrated, with the boundary conditions of (2.4) and (2.41), to obtain the distributions of the particle shear stress and the fluid horizontal velocity. Then, the particle horizontal velocity is given by (2.38), with (2.6), as

(2.45)\begin{equation}u = U - \frac{{\sqrt {r(Sh - s)} }}{\kappa }\exp [C_D^2(y + {y_0})]Ei[ - C_D^2(y + {y_0})] - \frac{5}{3}\frac{s}{{{C_D}{\sigma _y}}}.\end{equation}

Finally, the particle flux q per unit basal area is determined through numerical integration of (2.32), with the profiles of u and c obtained under turbulent conditions.

The above equations permit the determination of profiles of particle concentration, particle and fluid velocities and particle and fluid stresses in turbulent shearing flow that will be compared with the results of DC numerical simulations in § 4.

3. Scaling laws for rarefied saltation

The semi-analytical solutions to particle saltation in shearing flows described in the previous section permit simple asymptotic scalings to be obtained in the limit of rarefied saltation, $M \to 0$ and ${s_0} \ll Sh$, and: (i) only Stokes drag, ${C_D} = 18/St$, that is when the particle Reynolds number based upon the relative velocity between the particles and the fluid at the bed, ${u_0}St/r$, is less than unity; or (ii) only form drag, ${C_D} = 0.3{u_0}/r$, that is when ${10^3} < {u_0}St/r < {10^5}$. In the following, we will consider that at leading order, (2.11) and (2.22) imply that ${T_{y,0}} \propto h \propto T_{y,0}^{{ + ^{1/2}}}/{C_D}$.

We emphasize that there is no feedback of the particles on the fluid velocity profile, which, in rarefied saltation, is exactly linear in the viscous case, and logarithmic in the turbulent case. As a consequence, the particle hold-up only affects the particle concentration, is linearly proportional to it, and does not influence the particle velocity. Given that the mass flux involves the product of particle concentration and horizontal velocity, q must also be linearly related to the hold-up.

3.1. Rarefied saltation in viscous shearing flows with Stokes drag

In this case, the density ratio plays no role in the equations governing the saltation process. As mentioned, ${C_D} \propto S{t^{ - 1}}$. Then, with (2.26), we obtain that $\theta \propto S{t^{ - 2}}S{h^{ - 1}}$. With this, (2.5), (2.15), (2.19) and (2.36) imply the scalings for the various quantities that we report in table 1. In particular, we note that the scaling for the particle flux, $q \propto MS{t^4}S{h^{3/2}}$, was also derived in Valance & Berzi (Reference Valance and Berzi2022) using an approach in which all particles were assumed to follow the same PT, and not a distribution of trajectories as in the present work.

Table 1. Summary of the scaling laws for rarefied saltation.

4. Comparisons with numerical simulations and experiments

Here, we make comparisons between the predictions of the present theory and the results of DC numerical simulations of saltation of spheres over a rigid bumpy bed made of a layer of particles in close contact. The simulations are performed in a quasi-two-dimensional cell of streamwise length equal to 5120 particle diameters and transverse width equal to one particle diameter, with periodic boundary conditions in the streamwise direction. The cell is not bounded above. We checked that increasing the length in the streamwise direction up to ten times had no effect on the results. Data are available at https://doi.org/10.5281/zenodo.11264272 (Valance Reference Valance2024).

We solve Newton's equations of motion for the individual spherical particles under the influence of fluid drag, buoyancy, gravity and contact forces in collisions with the bed. Mid-trajectory collisions are forbidden. The fluid is treated as either a viscous or a turbulent flow, depending on the closure for the fluid shear stress, which is governed by a balance equation. The sum over all particles of drag and buoyancy enters the momentum balance for the fluid with a change in sign, thus ensuring the two-way coupling between the phases. In the numerical simulations, we assume that the fluid possesses only horizontal velocity and suppress the possibility of interparticle collisions above the bed.

The individual particles experience vertical and horizontal components of fluid drag based upon the local difference between the instantaneous velocity of the particle and the average velocity of the fluid. As in the theoretical treatment, we control the amount of particles in the simulations (the particle hold-up, M); the bumpiness of the rigid bed (the wall-particle diameter, d_w); the fluid viscosity (the inverse of the fall Stokes number, St); the fluid mass density (the inverse of the density ratio, r); and the intensity of the shearing flow (the Shields number, Sh). All measurements have been taken once a steady state is attained and, subsequently, time averaged. Such numerical simulations have already been used in the context of saltation in both turbulent (Durán et al. Reference Durán, Andreotti and Claudin2012; Pähtz et al. Reference Pähtz, Durán, Ho, Valance and Kok2015; Pähtz & Durán Reference Pähtz and Durán2020; Ralaiarisoa et al. Reference Ralaiarisoa, Besnard, Furieri, Dupont, Ould El Moctar, Naaim-Bouvet and Valance2020) and viscous (Valance & Berzi Reference Valance and Berzi2022) shearing flows. A more detailed description of the numerical simulations, including the contact parameters that we employ, can be found in Appendix D.

In a total of 50 simulations, we have numerically investigated the saltation process of four different types of solid particles in terrestrial and extra-terrestrial environment, as summarized in table 2, by changing the particle hold-up between 0.0004 and approximately 0.0300, that is, from the rarefied limit to the maximum hold-up for which the mid-trajectory collisions can be neglected and (2.23) is satisfied. Although the fluid regime on Mars and Venus is almost certainly turbulent, one can at least imagine performing experiments in pressurized wind tunnels, in which the turbulence is somehow suppressed, thus recovering the conditions reported in the first two rows of table 2. The conditions reported in the last two rows of table 2 are, instead, much closer to actual physical applications.

Table 2. Summary of the combination of parameters employed in the numerical simulations.

The flow conditions chosen serve to isolate and test the assumptions that we have made in building the theory: (i) an anisotropic Maxwellian velocity distribution for the particles, the vertical velocity of the ascending particles much larger than the settling velocity and a locally linear velocity for the fluid in the viscous regime, with Stokes drag (the first row in table 2); (ii) the additional assumption of drag coefficient uniform and equal to that evaluated at the bed in the viscous regime, with nonlinear drag (the second and third rows in table 2); and (iii) the fluid velocity profile encountered by ascending particles uniform and equal to the average of the logarithmic profile for the turbulent regime, with nonlinear drag (the fourth row in table 2).

The parameters of the turbulent case in table 2 match those of physical experiments of saltation on rigid, bumpy beds performed in a wind tunnel (Ho Reference Ho2012). For these, measurements of particle mass flux and profiles of particle concentration and horizontal particle and fluid velocities at different values of the particle hold-up are available. In the experiments, unlike the numerical simulations, the vertical velocity of the fluid surrounding the particles is non-zero, due to the no-slip condition on the particle surface, and this would permit a test of its influence on the results. However, as shown later, the predicted depth of the saltation layer in the absence of an upper bound and that measured in the numerical simulations is of the order of 4000 particle diameters; while the experiments were performed in a rectangular closed conduit, with a horizontal lid placed at approximately 1200 particle diameters above the rigid base. As a consequence, the top boundary conditions are different from those of the present numerical simulations and the semi-analytical treatment. For this reason, we do not show the profiles of the experiments. Although the order of magnitude of the profiles measured in the experiments is in good agreement with the predictions, their shape reveals the influence of the upper boundary. We postpone to a future work the solution of the appropriate two-point boundary value problem, with our proposed constitutive relations for the particle stresses and energy flux, and detailed comparisons against the experimental measurements.

Figures 5 and 6 show the comparisons between profiles of T_y, c, σ_y, s, u and U relative to selected values of the hold-up for saltation in viscous shearing flows with Stokes drag (the first column of plots in both figures), saltation in viscous shearing flows with nonlinear drag (the second column) and saltation in turbulent shearing flows with nonlinear drag (the third column). A similar agreement between the predictions and the simulations, not shown here for brevity, is obtained for all admissible values of M.

Figure 5. Profiles of (a–c) mean square of particle velocity fluctuations in the vertical direction, (d–f) particle concentration and (g–i) particle normal stress along y measured in numerical simulations of: saltation in viscous shearing flows with d_w = 1, St = 100, r = 150 000, Sh = 0.05 (first column) and M = 0.0008 (blue circles), M = 0.0033 (red circles), M = 0.0131 (purple circles); saltation in viscous shearing flows with d_w = 1, St = 100, r = 50, Sh = 0.05 (second column) and M = 0.0008 (blue squares), M = 0.0033 (squares), M = 0.0131 (purple squares); saltation in turbulent shearing flows with d_w = 1.5, St = 1681, r = 2208, Sh = 0.04 (third column) and M = 0.0008 (blue diamonds), M = 0.0033 (red diamonds), M = 0.0119 (purple diamonds). The solid lines are the predictions of the present theory.

Figure 6. Profiles of (a–c) particle shear stress, (d–f) fluid and (g–i) particle mean horizontal velocities measured in numerical simulations of: saltation in viscous shearing flows with d_w = 1, St = 100, r = 150 000, Sh = 0.05 (first column) and M = 0.0008 (blue circles), M = 0.0033 (red circles), M = 0.0131 (purple circles); saltation in viscous shearing flows with d_w = 1, St = 100, r = 50, Sh = 0.05 (second column) and M = 0.0008 (blue squares), M = 0.0033 (squares), M = 0.0131 (purple squares); saltation in turbulent shearing flows with d_w = 1.5, St = 1681, r = 2208, Sh = 0.04 (third column) and M = 0.0008 (blue diamonds), M = 0.0033 (red diamonds), M = 0.0119 (purple diamonds). The solid lines are the predictions of the present theory. The insets of figures 6(f) and 6(i) show the corresponding velocity profiles in semi-log scale.

In determining the semi-analytical solution that we have highlighted in § 2, we have employed α_u = 0.60 and α_T = 0.50 for saltation in viscous flows, and α_u = 0.85 and α_T = 0.18 for saltation in turbulent flows. We have fitted these values to exactly match the depth of the saltation layer, h, and the particle shear stress at the bed, s ₀, in one of the discrete simulations of saltation in viscous flows and one of the simulations of saltation in turbulent flows. Different values of these parameters would not alter the qualitative features of our semi-analytical solution, but impact the values of T_y, c and u at the bed, therefore causing a shift to the right or to the left in the relative profiles. As an overall estimate of their influence, we checked that a ten per cent change in either of the alphas would result in a ten per cent change in the mass flux.

Figures 5(a)–5(c) confirm that the intensity of the velocity fluctuations of the particles decreases linearly with the distance from the rigid bed, and that the model captures both the weak dependence of T_y on M for viscous saltation and Stokes drag (figure 5a), and its independence for nonlinear drag (figures 5b and 5c). The depth of the saltation layer is well predicted by the continuum model, but the slope of the linear decrease is underestimated with respect to the numerical simulations for nonlinear drag. Also, the non-monotonic behaviour of T_y near the rigid bed is not captured. These are the consequences of neglecting the dependence of the drag coefficient on the vertical direction. The agreement with the measurements in the numerical simulations would indeed improve if C_D is allowed to vary locally; however, this would prevent obtaining semi-analytical solutions of the governing equations, with little improvement on the results. The condition ${C_D}\sqrt {{T_y}} > 1$ that ensures that the assumption of large vertical velocities of ascending particles with respect to the settling velocity is sufficiently accurate, is satisfied in the flows of figures 5 and 6 everywhere, but for the upper quarter of the flow, where indeed the deviations between the simulations and the predictions are more significant.

In all cases, the continuum model satisfactorily reproduces the profiles of particle concentration and normal stress in the y-direction (figure 5d–i), which indeed follow the power-law distributions of (2.17) and (2.18). In the turbulent case, as anticipated, the particle concentration is rather uniform across the flow (figure 5f). The values of the particle concentration are also much smaller in the case of turbulent shearing flows (figure 5f), as a consequence of the greater agitation of the particles (figure 5c).

Figure 5 shows profiles of quantities associated with the vertical motion of the saltating particles, which is sensitive to the horizontal motion of the fluid only through the dependence of the drag coefficient on U. The horizontal motion of the fluid and, in particular, its flow regime, strongly affects the distribution of particle shear stress and particle and fluid horizontal velocity shown in figure 6.

The proposed constitutive relations for the particle shear stress ((2.25) and (2.38)) permit the qualitative and quantitative reproduction of the measurements in numerical simulations (figure 6a–c). In particular, the position and the magnitude of the peak in the particle shear stress are well captured. The physical reason for the presence of this maximum value of the particle shear stress above the bed is in the change of sign of the horizontal force exerted by the fluid on the particles. The particles slip at the rigid bed, while the fluid does not; hence, the particles drag the fluid near the bed, while the fluid drags the particles as the top of the saltation layer is approached. Interestingly, the particle shear stress can reach values up to 50 %–60 % of the total shear stress – the Shields parameter – for saltation in viscous shearing flows, but only up to 10 % of the total shear stress for saltation in turbulent shearing flows.

For saltation in viscous shearing flows (Figure 6d–h), the fluid and particle horizontal velocity profiles become progressively nonlinear as the particle hold-up increases; but the continuum model is capable of capturing this behaviour. For saltation in turbulent shearing flows, the fluid velocity profile is logarithmic, and the hold-up has little influence on it (figure 6f). The horizontal velocity of the particles is almost uniform vertically (figure 6i). This confirms our finding that, in the absence of particle interactions, the particle shear stress does not depend on the particle shear rate. The slight discrepancies between the predictions of the continuum model and the results of the numerical simulations are likely due to the assumption of a drag coefficient independent of y.

Finally, we compare the dependence of the particle flux on the particle hold-up in figure 7. The continuum model notably reproduces the measurements in numerical simulations in both the viscous and turbulent shearing flows. The experimental results on turbulent saltation in a wind tunnel (Ho Reference Ho2012) are in relatively good agreement with both the numerical simulations and the continuum model (figure 7b). The overestimate of the experimental mass flux is likely due to the additional resistance induced by the presence of the horizontal lid above the rigid bed. Figure 7(a) also assesses the validity of the scaling for q reported in table 1 for the case of viscous saltation and Stokes drag, in the rarefied limit.

Figure 7. (a) Particle flux against scaled particle hold-up measured in DC numerical simulations (open symbols) and predicted from the present theory (lines) in the case of saltation in viscous shearing flows with d_w = 1, r = 150 000 and: St = 60 and Sh = 0.05 (blue circles); St = 100 and Sh = 0.05 (orange circles); St = 100 and Sh = 0.075 (purple circles). The inset depicts the same data in terms of q against M. (b) Particle flux against particle hold-up measured in DC numerical simulations (open symbols) and in wind-tunnel experiments (filled symbols, after Ho Reference Ho2012) and predicted from the present theory (lines) in the case of saltation in viscous shearing flows with d_w = 1, St = 100 Sh = 0.05 and r = 2.5 (orange squares) and r = 50 (blue squares); and in the case of saltation in turbulent shearing flows with d_w = 1.5, St = 1681, r = 2208 and Sh = 0.025 (purple diamonds) and Sh = 0.04 (green diamonds).