domain: phys sde

text:
Diego Berzi

Alexandre Valance

James T Jenkins

Collisionless

Collisionless kinetic theory for saltation

Introduction

Saltation, that is the motion of solid particles in a gravitational field through successive jumps and rebounds from a base driven by a shearing flow, has been established as the main mode of transport of windblown sand [START_REF] Andreotti | A twospecies model of aeolian sand transport[END_REF][START_REF] Bagnold | The Physics of Blown Sand and Desert Dunes[END_REF][START_REF] Charru | Sand Ripples and Dunes[END_REF][START_REF] Owen | Saltation of uniform grains in air[END_REF][START_REF] Valance | The physics of Aeolian sand transport[END_REF] and is crucial in the dynamics of dunes [START_REF] Sauermann | A Continuum Saltation Model for Sand Dunes[END_REF]. Although originally conceived as the transport of sand grains by the wind on Earth, it has been recognized as significant also for the transport of sand or gravel in water on Earth [START_REF] Abbott | Saltation and Suspension Trajectories of Solid Grains in a Water Stream[END_REF][START_REF] Ancey | Saltating motion of a bead in a rapid water stream[END_REF][START_REF] Fernandez Luque | Erosion And Transport Of BedLoad Sediment[END_REF][START_REF] Niño | Experiments on saltation of sand in water[END_REF], basalt particles on Venus [START_REF] Greeley | Windblown sand on Venus Preliminary results of laboratory simulations[END_REF][START_REF] Iversen | Saltation threshold on Earth, Mars and Venus[END_REF] and on Mars [START_REF] Iversen | Saltation threshold on Mars: The effect of interparticle force, surface roughness, and low atmospheric density[END_REF][START_REF] Iversen | Saltation threshold on Earth, Mars and Venus[END_REF], and ice particles on Titan [START_REF] Burr | Higherthanpredicted saltation threshold wind speeds on Titan[END_REF]. While these studies considered the carrier fluid to be turbulent, saltation in viscous, nonturbulent fluids has also been investigated [START_REF] Charru | Instability of a bed of particles sheared by a viscous flow[END_REF][START_REF] Ouriemi | Sediment dynamics. Part 1. Bedload transport by laminar shearing flows[END_REF][START_REF] Seizilles | Crossstream diffusion in bedload transport[END_REF], given its practical relevance in limiting the transport capacity of oil in pipes [START_REF] Dall'acqua | Experimental results of pipeline dewatering through surfactant injection[END_REF][START_REF] Leporini | On the numerical simulation of sand transport in liquid and multiphase pipelines[END_REF].

There is a large body of mathematical models of saltation in the literature (see, e.g., the reviews of [START_REF] Kok | The physics of windblown sand and dust[END_REF]Pähtz et al. 2020;[START_REF] Valance | The physics of Aeolian sand transport[END_REF]. As highlighted by [START_REF] Valance | The physics of Aeolian sand transport[END_REF], the different models can be grouped into two main categories, Lagrangian or Eulerian approaches, based upon the point of view adopted in the description of the particle motion.

In many Lagrangian approaches [START_REF] Anderson | Simulation of Eolian Saltation[END_REF][START_REF] Creyssels | Saltating particles in a turbulent boundary layer: experiment and theory[END_REF][START_REF] Kok | A comprehensive numerical model of steady state saltation (COMSALT)[END_REF][START_REF] Werner | A steadystate model of windblown sand transport[END_REF]), an assigned distribution of particle ballistic trajectories is calculated by solving their differential momentum balances, governed by fluid drag and gravity, when immersed in a fluid with a prescribed velocity profile and for a given distribution of initial takeoff velocity at the basal boundary. The calculation is repeated once the fluid velocity profile and the distribution of the takeoff velocities are updated using a simple constitutive relation for the fluid shear stress (in the case of turbulent fluid, usually based on a mixing length approach) and a suitable set of boundary conditions that govern the impact of the particles at the base [START_REF] Beladjine | Collision process between an incident bead and a threedimensional granular packing[END_REF][START_REF] Crassous | Impact of a Projectile on a Granular Medium Described by a Collision Model[END_REF][START_REF] Oger | Discrete Element Method studies of the collision of one rapid sphere on 2D and 3D packings[END_REF], until a steady state is attained.

Simpler Lagrangian approaches in which the distribution of the particle trajectories is replaced by one [START_REF] Jenkins | Periodic trajectories in aeolian sand transport[END_REF] or two [START_REF] Andreotti | A twospecies model of aeolian sand transport[END_REF]) modes have also been recently proposed and successfully compared against experiments and numerical simulations. In particular, the one species models, also called Periodic Trajectory (PT) models, have been applied to saltation in both turbulent [START_REF] Berzi | Periodic saltation over hydrodynamically rough beds: aeolian to aquatic[END_REF] and viscous [START_REF] Valance | Particle saltation over rigid bumpy beds in viscous shearing flows[END_REF] shearing flows, for values of the ratio of the graintofluid mass densities encountered on Earth and other planetary bodies [START_REF] Berzi | The threshold for continuing saltation on Earth and other solar system bodies[END_REF]. PT models are simple enough to allow even for fully analytical solutions of steady and fully developed saltation, and permit the determination, with sufficient accuracy, of global quantities, such as the particle mass flux, as a function of the intensity of the shearing flows. On the other hand, many variables of interest, such as the distribution with height of the particle concentration and the total depth of the saltation layer, are poorly predicted. It is also not obvious how to extend this model to unsteady and/or boundaryvalued (inhomogeneous) problems.

Another family of Lagrangian approaches is based on the framework of Discrete Element Method (DEM, [START_REF] Cundall | A discrete numerical model for granular assemblies[END_REF], and solves the Newton's laws of motion for the individual grains that are allowed to collide with other particles and with the base, while the surrounding fluid is replaced by discrete forces (such as drag and buoyancy) acting on the particles themselves [START_REF] Tsuji | Discrete particle simulation of twodimensional fluidized bed[END_REF]. The forces transmitted by the fluid on the particles are then introduced, after a change of sign, in the fluid momentum balance, ensuring an instantaneous two way coupling. These DiscreteContinuum (DC) numerical simulations are a powerful tool, in that they greatly reduce the number of assumptions necessary to solve for the transport process. For instance, the dynamics of the impact of the particles with the base must not be modelled in advance, but is an output of the simulations. Likewise, the possibility of interparticle collisions above the base is naturally accounted for. On the other hand, the number of particles that is feasible to simulate is severely limited by computational power and nowhere near to the actual number of grains involved in realscale applications. DC simulations have been applied to saltation in turbulent [START_REF] Durán | Numerical simulation of turbulent sediment transport, from bed load to saltation[END_REF][START_REF] Pähtz | The fluctuation energy balance in nonsuspended fluidmediated particle transport[END_REF]Pähtz and Durán 2020b;[START_REF] Ralaiarisoa | Transition from Saltation to Collisional Regime in Windblown Sand[END_REF]) and viscous flows [START_REF] Valance | Particle saltation over rigid bumpy beds in viscous shearing flows[END_REF] and provide a large quantity of measurements, some of which are simply unattainable in physical experiments. Hence, they serve as severe tests of more sophisticated approaches.

Eulerian approaches, in which both the fluid and the saltating particles are treated as two superimposed continuum phases offer the most promising perspective for modelling large scale phenomena. In the model of [START_REF] Sauermann | A Continuum Saltation Model for Sand Dunes[END_REF], the motion of the two phases is depthaveraged over the saltation layer, and cannot therefore predict the distribution of the variables of interest with height. A key point of the model is that the wind profile in the saltation layer is determined by assuming an exponential distribution of the particle shear stress there. The same assumption is adopted in more recent Eulerian models [START_REF] Lämmel | A twospecies continuum model for aeolian sand transport[END_REF][START_REF] Pähtz | The apparent roughness of a sand surface blown by wind from an analytical model of saltation[END_REF]. Assuming an a priori distribution is crucially different from phrasing a constitutive relation for the particle shear stress, and obtaining the distribution from the usual momentum balance. [START_REF] Jenkins | Continuum model for steady, fully developed saltation above a horizontal particle bed[END_REF] were the first to propose a constitutive relation for the particle shear stress based on averaging the equations governing the trajectory of the single particles and substituting averaging of products with products of averaging. They also equated the particle pressure to the product of the particle concentration and the granular temperature, the mean square of the particle velocity fluctuations. They assumed, as in [START_REF] Creyssels | Saltating particles in a turbulent boundary layer: experiment and theory[END_REF], that the latter is uniformly distributed with the distance from the base, which holds only if the vertical drag on the particles is negligible. From the uniform distribution of the granular temperature, they obtained the exponential decay of the particle concentration and distributions of particle and fluid mean horizontal velocities that agreed reasonably well with experiments [START_REF] Creyssels | Saltating particles in a turbulent boundary layer: experiment and theory[END_REF]. The model was later extended to deal with unsteadiness and inhomogeneities [START_REF] Jenkins | Twophase continuum theory for windblown sand[END_REF]. Interestingly, the constitutive relation for the particle shear stress of [START_REF] Jenkins | Continuum model for steady, fully developed saltation above a horizontal particle bed[END_REF] coincides with the dilute and collisionless limit of the expression derived by [START_REF] Garzó | Enskog kinetic theory for monodisperse gas solid flows[END_REF] by solving the Enskog kinetic theory for monodisperse gassolid flows. This is an indication that both constitutive relations share the same assumption, namely that they only account for the influence of the fluid drag and not gravity nor buoyancy on the particle path.

Here, we formulate a collisionless kinetic theory for saltating particles over an horizontal bed. We distinguish between ascending and descending particles and assume that the velocity distribution function of the ascending particles is an anisotropic Maxwellian, as in [START_REF] Creyssels | Saltating particles in a turbulent boundary layer: experiment and theory[END_REF]. Then, we employ approximate analytical expressions for the particle trajectories in viscous and turbulent shearing flows to determine the transfers of horizontal and vertical momentum, and kinetic energy associated with the particle vertical motion across an horizontal surface. Once averaged over all possible particle trajectories using the velocity distribution function, we obtain approximate analytical constitutive relations for the particle stresses and the energy flux of fluctuation kinetic energy.

As suggested by [START_REF] Pasini | Aeolian transport with collisional suspension[END_REF], collisions above the base can be ignored if the mean free path of kinetic theory is larger than the length of the ballistic trajectory. Given that the mean free path strongly decreases if the particle concentration increases [START_REF] Chapman | The mathematical theory of nonuniform gases The mathematical theory of nonuniform gases[END_REF], and that the particle concentration near erodible beds (composed of particles at rest identical to those in saltation) is large [START_REF] Tholen | Anomalous Scaling of Aeolian Sand Transport Reveals Coupling to Bed Rheology[END_REF], we argue that, strictly speaking, collisionless saltation is possible only if it takes place over rigid beds.

We phrase and solve the balances of fluid and particle momenta and particle fluctuation energy for a steady, fullydeveloped saltation over an horizontal, rigid, bumpy bed, driven by either a viscous or a turbulent shearing flow, in the boundary layer approximation, with appropriate boundary conditions. As in all previous models of saltation, we assume that the fluid velocity is solely in the horizontal direction. This permits to uncouple the determination of quantities associated with the particle vertical motion, such as intensity of the vertical velocity fluctuations, particle concentration and normal stress, from those that involve the particle horizontal motion and are, therefore affected by the flow regime of the fluid, such as particle shear stress and particle and fluid mean horizontal velocity. Depending on the flow regime of the fluid, only one or two differential equations must be solved numerically to determine the vertical profiles of the corresponding quantities. The rest of the profiles, including the depth of the saltation layer, are obtained analytically.

After deriving scaling laws for, e.g., the particle mass flux in some special limits, we successfully test the results of the theory against quasi2D DC simulations of saltation in viscous and turbulent shearing flows that were carried out suppressing the possibility of midtrajectory collisions for a wide range of particleto fluid mass density ratio, ranging from Mars to Earth environments, at different values of the intensity of the shearing flows, the amount of particles in the system and the viscosity of the fluid. The relation between the particle mass flux and the amount of particles in the system in windtunnel experiments [START_REF] Ho | Etude expérimentale du transport de particules dans une couche limite turbulente[END_REF]) is also satisfactorily captured.

In Section II, we present the constitutive relations, the balance equations and the boundary conditions that we employ to obtain semianalytical solutions to steady and fully developed saltation over rigid, bumpy beds. In Section III, we show how to obtain scaling laws for various quantities in the limit of rarefied saltation in both viscous and turbulent flows, at least when the drag coefficient reduces to its asymptotic expressions. Comparisons against DC simulations and experiments are detailed in Section IV. Finally, in Section V, we conclude with a summary of the main findings and an outline of future endeavors.

II. Governing equations and semianalytical solutions

The saltating particles are identical spheres of diameter d and mass density  s . A shearing flow of a fluid of mass density  f and molecular viscosity  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. 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 dw, in close contact to each other, over a flat plate. Hence, the ratio of dw to d is a natural measure of the bed roughness. We characterize the particles through the fall particle Reynolds number

( ) 3/ 2 R 1 / / f f g r rd   = -
, 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(r1)/r. Then, lengths, velocities and stresses are expressed in units of d, [g(r1)d/r] 1/2 and  s g(r 1)d/r, respectively. The balance equations for the x and ymomenta of the particles, under steady and fullydeveloped conditions, are:

d ; d y c y  = - (1) ( ) d , d D s cC U u y = - - (2)
where: y  and s are the particle normal stress in the ydirection and the particle shear stress, respectively; c is the particle volume concentration; and CD 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 that CD is independent of y and equal to its average over all the spherical particles evaluated at the bed,

0 0 18 0.3 , St D C U u r = + - (3)
where St = rR is the fall Stokes number, and U0 and u0 are the fluid and particle horizontal velocities evaluated at the bed (here and in what follows, the subscript 0 refers to quantities evaluated at the bed). We employ the term Stokes drag to refer to situations in which 18 St

D C

, form drag for situations in which

0 0 0.3 D C U u r -
and nonlinear drag for the generic case in which both Stokes and form drag must be accounted for. A more appropriate expression for the average CD should also involve the strength of the particle velocity fluctuations [START_REF] Jenkins | Collisional sheet flows of sediment driven by a turbulent fluid[END_REF]. 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. For flows over rigid beds, we take the noslip boundary condition for the fluid,

0 0. U = (4)
If saltating particles with a large value of the coefficient of sliding friction impact a rigid bed, made bumpy by gluing identical spheres in close contact on a flat plate, at an average small angle  with respect to the horizontal (Figure 1), then, using the analysis of [START_REF] Lämmel | Grainscale modeling and splash parametrization for aeolian sand transport[END_REF]

in Appendix C, 0 1 , u D u C   = (5)
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.

In the boundary layer approximation, the sum of the fluid and particle shear stress, S and s respectively, is constant and equal to the farfield fluid shear stress, which, in dimensionless terms, is the Shields parameter:

Sh. s S + = (6)
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 ydirection is simply given by the average vertical flux of ymomentum, ,

y y cT  = (7)
where y T is the mean square of the vertical velocity fluctuations of the particles.

The distribution of the additional hydrodynamic field, Ty, along y is governed by the balance of kinetic energy associated with the particle vertical motion, that is, the yycomponent of the particle secondmoment tensor, which, in a steady and fullydeveloped flow, and in the absence of midtrajectory collisions, reduces to (Saha andAlam 2016, 2017)

d 2 0. d yyy D y Q C y  - - = (8)
Here, yyy Q is the yyyelement of the flux of second moment, Q , a third rank tensor, and 2 D y C  is the energy 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. With the further assumptions that the velocity distribution of the ascending particles is an anisotropic Maxwellian [START_REF] Creyssels | Saltating particles in a turbulent boundary layer: experiment and theory[END_REF]) and that the vertical velocity of the ascending particles is much larger than 1/CD (the settling velocity), in Appendix B we derive a simple expression for yyy Q :

( ) 3/2 8 . yyy y Q c T  + + = (9)
where c + and Ty + are the volume concentration and the mean square of the vertical velocity fluctuations of the ascending particles, respectively. At least under steady state conditions, c + and Ty + are related to c and Ty through (see Appendix B):

3 2 ; 4 y D T c c C  + +      = +     (
  + + + + +   + = - -     (13)
whose analytical solution, with the boundary condition

( ) ,0 D y y D D y D y D D y y D D y D y C T T C iC T iC T C C y T T C iC T iC T                - + + + + - + + + +   - -       - - + - - + + -       = - + - -       - - + - - + + -     (14)
Under the same assumptions that we employ to derive equation ( 5), a simple dependence of ,0 y T + calculated at the bed on the impact angle  can be obtained (Appendix C) in the form:

,0 2 2 1 , 1 w y T D w d T C d   + = + ( 15
)
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 ,0 y T + , remains an unknown at this stage of the analysis. We will show how to determine it in the next two subsections, when we describe quantities associated with the horizontal motion of the particles, and, therefore, influenced by the flow regime of the fluid. Equation ( 14) provides the analytical distribution of Ty + along y that, inserted into equation ( 11), would also give the analytical distribution of Ty. As shown in Figure 2, the latter is roughly linear and well approximated as ,0 ,

y y h y T T h - = (16)
where the depth of the saltation layer, h, is to be determined. The mean square of the particle velocity fluctuations at the bed, Ty,0, can be obtained from equation (11), if ,0 y T + is known. Inserting equation [START_REF] Guo | Buffer law and transitional roughness effect in turbulent openchannel flows[END_REF] into equation ( 1), with equation ( 16), and integrating gives the following powerlaw distribution of the particle concentration, ,0 ,0 0 ,

y y h T T h y c c h - -   =     ( 17
)
where c0 is the particle concentration at the bed. Integrating equation (1) with equation ( 17) also gives

,0 0 ,0 . y h T y y h y c T h  -   =     (18)
After introducing the particle mass holdup, 0 h M cdy =  , i.e., the particle mass per unit basal area, by integrating equation (1) over the saltation layer, with y = 0 at y = h, we obtain

0 ,0 . y M c T = (19)
Integrating the energy balance (equation 8) with equations ( 18) and ( 19), provides

,0 ,0 ,0 ,0 ,0 2 1 , y y h T T yyy D y yyy y h h y Q C MT Q h T h +   -     = -+     +       (20)
where the value of Qyyy at the bed can be obtained from the ratio of equations ( 9) and ( 7), with equations (10), ( 11) and ( 19), as

,0 ,0 ,0 4 . 2 2 D y yyy D y C T Q M C T  + + = + (21)
Given that Qyyy must vanish at the top of the saltation layer, the depth h is determined from equation (20) as ,0 ,0 ,0 ,0 . 2

yyy y D y yyy Q h T C MT Q = - (22)
The governing equations, constitutive relations and analytical results described so far apply, in general, to both saltation in viscous and in turbulent shearing flows and can be calculated only after the determination of the impact angle . To proceed, we now must distinguish between the two regimes of the shearing flow.

Before doing that, we support the assumption regarding particle collisions above the bed. As suggested by [START_REF] Pasini | Aeolian transport with collisional suspension[END_REF], chances of midtrajectory 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 [START_REF] Chapman | The mathematical theory of nonuniform gases The mathematical theory of nonuniform gases[END_REF], evaluate this at the bed, and take h as the average height of the particle trajectories, we obtain

0 0 1 2 . 6 2 h c  =  (23)
Given that, from equation ( 22),

,0 y h T , equation ( 23) implies that there is a maximum mass holdup M (for given values of r, St and Sh) above which midtrajectory collisions play a role. The mean free path is a decreasing function of the particle concentration. If the saltation process were to take place over an erodible -a densely packed assembly of particles identical to those in saltation that can be mobilized due to impacts rather than a rigid bed, and if the volume concentration is continuous across the interface with the erodible bed, then the theoretical mean free path calculated at the top of the erodible bed, where 0.6 c , would be much less than the mean free path evaluated at the rigid bed, where 1 c

. Hence, we expect that the mass holdup at which saltating particles over a rigid bed start to deposit and form an erodible substrate is larger than the maximum mass holdup for which midtrajectory collisions can be ignored.

Saltation in viscous shearing flows

In the absence of turbulence, the expression for the fluid shear stress is [START_REF] Valance | Particle saltation over rigid bumpy beds in viscous shearing flows[END_REF])

1 d 1 d , St d St d c U U S y y - = (24)
where we have neglected the particle concentration with respect to unity. Equation ( 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. To solve for the particle trajectory, 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 at least locally linear, in the sense that is linear in the region comprised between y and the top of the trajectory. However, we let the slope of the linear profile to change with y. This assumption and the assumptions that we have already employed in the derivation of the expression for Qyyy in equation ( 9) permit the derivation in Appendix B of a simple expression for the particle shear stress,

2 3 2 1 d . 5 3 d D y D U s C U u C y    = -+     (25)
Interestingly, the particle shear stress does not depend on the particle shear rate, as in [START_REF] Jenkins | Continuum model for steady, fully developed saltation above a horizontal particle bed[END_REF], 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, the angle  between the velocity of the particles before the impact with the bed and the horizontal for saltation in viscous flows is estimated as

( ) 0 , St Sh D C s  = - ( 26
)
where s0 is the particle shear stress at the bed.

Using equations ( 4) and ( 5) in equation ( 3), with equation ( 26), we obtain a relationship between the drag coefficient and the particle shear stress at the bed, ( )

3 2 0 18 0.3 St Sh . St D D u C C s r  = + - (27)
The particle shear stress at the bed is obtained from equation ( 25), with equations ( 4), ( 5), ( 6), ( 19) and ( 24), as

( ) ( ) 0 2 3 StSh , 5 2 3 St u D u M s C M   - = + - (28)
The system of equations ( 27) and ( 28) can be solved to determine the drag coefficient and the particle shear stress at the bed once the values of dw, r, St, Sh and M are given. Then, the impact angle (equation 26), the values of all the variables at the bed (equations 5, 15 and 21), the depth of the saltation layer (equation 22), and the analytical distributions of Ty + , Ty, c , y and Qyyy (equations 14, 1618 and 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 mass holdup, as predicted by equations ( 2628), for, e.g., saltation of 100 m basalt grains in viscous shearing flows on Venus, with u = 0.6 (see Section IV for more details about the choice of the parameters). The range of the mass holdup in Figure 3 is that for which equation ( 23) is satisfied. Notice that for almost the entire range of mass holdup, the particle shear stress at the bed linearly increases with Using equation ( 25) in equation ( 2), with equations ( 6), ( 18) and ( 24), gives a first order, linear, non homogeneous, differential equation for the particle shear stress,

,0 ,0 ,0 ,0 1 0 0 ,0 St StSh d 5 1 2 2 . d 3 3 3 y y y y h T h T T T y D D c c s h y h y h y s y T h C h C h - - -   - - -         = - + +                   (29)
The analytical solution of equation ( 29) is:

,0 ,0 ,0 ,0 5 5/3 3 5 5/3 3 2 St 2 2 St 2 St Sh , exp ln 3 3 3 3 2 St 2 2 St 2 St Sh , exp ln 3 3 3 3 y y y y h h h T T T D D D h T D D D M M h y h y M h y s C C h h C h M M h y M h y C C h C h         - - -             =  - +                                   - -       -  - +                    ,0 ,0 ,0 5 3 0 2 St 2 St exp ln , 3 3 y y y h T h h T T D D h y M h y M s h C h C              - -         + + -                 (30)
in which  is the incomplete gamma function.

The fluid velocity profile is determined by integrating the constitutive relation (equation 24)

( )

d d St Sh U y s =
-, with the particle shear stress given by equation (30). Unfortunately, there is no general analytical solution. Hence, we obtain the distribution of U by numerically solving this ordinary differential equation, with the noslip condition, equation ( 4), at the bed. Once U is determined, the particle horizontal velocity is given by equation ( 25) as ( )

2 St Sh 5 2 . 3 3 D y D s s u U C C  - = - + (31)
We can now calculate the particle mass flux q per unit basal area through numerical integration as

0 d . h q cu y =  (32)
In the special case of rarefied saltation, that is for 0 M → and Sh s , the analytical solution of equation ( 29) simplifies to ,0 ,0 ,0

5 5 3 3 0 StSh , y y y h h h T T T D h y M h y h y s s h C h h   - - -         = + -                 (33)
and the fluid velocity profile is simply ( )

StSh . U y = (34)
Then, from equations ( 28), ( 31), ( 33) and ( 34), and the expression for s0, the particle horizontal velocity is

,0 2 3 2 2 5 2 3 StSh 1 . 3 y h T u D D h y u C y h C    -+ -     = + -         (35)
Notice that in the rarefied limit of saltation in viscous shearing flows, the particle and the fluid velocity profiles are independent of the mass holdup. Using equations ( 17), ( 19) and ( 35) in equation ( 32), and integrating, gives ,0 ,0

StSh .

y y h q M T h T = + (36)

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 [START_REF] Jenkins | Continuum model for steady, fully developed saltation above a horizontal particle bed[END_REF],

( ) ( ) 2 2 2 2 2 0 2 0 1 , y y c dU dU S y y r dy r dy   +     - = +         (37)
where  = 0.41 is Von Kármán's constant and 0 1 St Sh 1 exp 9St Sh 30 26

s s k r y k r     = + - -            
is an expression for the origin of the logarithmic fluid velocity profile that encompasses hydrodynamically smooth, rough and transitional beds [START_REF] Guo | Buffer law and transitional roughness effect in turbulent openchannel flows[END_REF]. The roughness length scale ks is taken to be ks = dw [START_REF] Ho | Etude expérimentale du transport de particules dans une couche limite turbulente[END_REF]. Equation (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 that the fluid velocity profile encountered by an ascending particle during its ballistic trajectory above a certain position y is uniform and equal to the average of the logarithmic profile in the region comprised between y and the top of the trajectory. This average value of U changes with y. This assumption, and the assumptions that we have already employed in the derivation of the expression for Qyyy, equation ( 9), permits the derivation in Appendix B of an expression for the particle shear stress,

3 3 ln 1 . 5 4 y D y D T rS dU s C U u C rS dy    +         = -+ +         (38)
As in the viscous case of equation ( 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 in the flow domain.

As explained in Appendix C, the impact angle  for saltation in turbulent flows is estimated as

( ) 1 ,0 0 0 Sh 1 0.85 ln . y D D T r s C C y   - +     -             (39)
Using equation ( 5) in equation ( 3) gives the drag coefficient in terms of the impact angle,

2 9 1 18 1.2 1 . St 2 St D u C r     = + +     (40)
Equation ( 38) evaluated at the bed, with equations ( 5) and (39) also gives the particle shear stress at the bed as a function of the impact angle, ( )

0 3 0.9 , 5 u M s   = - (41)
Then, equation ( 39) with equations ( 15), ( 40) and ( 41) results in a transcendental equation for , ( )

1 2 2 2 0 1 9 1 18 1.2 1 0.85 3 1 Sh 0.9 St 2 St 5 2 1 9 1 18 1.2 1 ln , St 2 St 1 u u w T u w r M r d y r d           - -         + + - -                        + +       +         (42)
that can be solved to determine the impact angle. Once  is known, CD (equation 40), s0 (equation 41), the boundary values of the remaining variables at the bed (equations 5, 15 and 21), the depth of the saltation layer (equation 22), and the analytical distributions of Ty + , Ty, c , y and Qyyy (equations 14, 1618 and 20) can also be calculated for 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 mass holdup, as predicted by equations ( 4042), for, e.g., saltation of 240 m sand grains in a turbulent wind on Earth, with u = 0.8 and T = 1 (see Section IV for more details about the choice of the parameters). The range of the mass holdup in Figure 4 is that for which equation ( 23) is satisfied. As for saltation in viscous shearing flows (Figure 3), the particle shear stress at the bed linearly increases with Equation ( 2), with equations ( 38) and ( 6), and equation ( 37) can be written as a system of two ordinary differential equations:

( )

0 Sh d 3 1 5 ln 1 ; d 4 3 , y D D y T r s s s C c y C y y T  +   -   = + -   +   (43)
and

( ) ( ) 0 Sh d ; d r s U y y y  - = + (44)
these can be numerically integrated, with the boundary conditions of equations ( 4) and ( 41), to obtain the distributions of the particle shear stress and the fluid horizontal velocity. Then, the particle horizontal velocity is given by equation ( 38), with equation ( 6), as

( ) 0 Sh 3 1 5 ln 1 . 4 3 y D D y T r s s u U C y y C   +   -   = + + -   +   (45)
Finally, the particle mass flux q per unit basal area is determined through numerical integration of equation ( 32), with the profiles of u and c obtained under turbulent conditions.

III. Scaling laws for rarefied saltation

The semianalytical solutions to particle saltation in shearing flows described in the previous section permit simple scalings to be obtained, at least in the extreme cases of rarefied saltation, that is for 0 M → and In the following, we will consider that at leading order, equations ( 11) and ( 22) imply that

1/2 ,0 ,0 y y D T h T C +   .
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 mass holdup only affects the particle concentration, and it 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 mass holdup.

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,

1 St D C - 
. Then, with equation ( 26), we obtain that

2 1 St Sh  - - 
. With this, equations ( 5), ( 15), ( 19) and

(36) imply the scalings for the various quantities that we report in Table 1. In particular, we notice that the scaling for the mass flux,

4 3/2
St Sh q M  , was also derived in [START_REF] Valance | Particle saltation over rigid bumpy beds in viscous shearing flows[END_REF] by using an approach in which all particles were assumed to follow the same periodic trajectory, and not a distribution of trajectories as in the present work.

Rarefied saltation in viscous shearing flows with form drag

In the case of form drag,

1 0 D C u r - 
. Then, we expect the scaling laws to involve also the density ratio. With equation ( 5), we obtain

1/2 1/2 0 u r  - 
, and, with equations ( 26) and ( 15),

1/3 2/3 2/3 St Sh r  - - -  and ,0 y T r + 
. Hence, equations ( 5), ( 15), ( 19) and ( 36) imply the scalings for the various quantities that we report in Table 1.

Rarefied saltation in turbulent shearing flows with form drag

In this limiting case, the fall Stokes number, St, plays no role in the equations governing the saltation process and cannot be involved in the scalings. With

1 0 D C u r - 
and equation ( 5), we obtain

1/2 1/2 0 u r  - 
, and, with equations ( 41) and ( 15), at leading order,

1 Sh  -  and ,0 y T r +  . Hence, 1/2 1/2 0 Sh u r  , 1/2 1/2 Sh D C r - 
, and

1/2 1/2 ,0 ,0 Sh y y D h T T C r + -   
. Equation ( 19), then, provides

1 1/2 0 Sh c Mr -  . As
shown later, the concentration and particle velocity profiles in the case of turbulent saltation are approximately uniform along y. Therefore, we expect

1/2 1/2 0 0 Sh q c u h M r  
. We summarize the scaling laws for rarefied turbulent saltation with form drag in Table 1.

St - 1/3 1/3 1/3 St Sh r - 1/2 1/2 Sh r - Impact angle,  2 1 St Sh - - 2/3 2/3 1/3 St Sh r - - - 1 Sh - Particle slip velocity, u0 3 St Sh 1/3 1/3 2/3 St Sh r 1/2 1/2

Sh r

Depth of saltation layer, h

3 1/2 St Sh 1/3 1/3 5/6 St Sh r - - 1/2 Sh r - Particle concentration at the bed, c0 3 1/2 St Sh M -- 1/3 1/3 5/6 St Sh M r - 1/2 1 Sh M r - Particle mass flux, q 4 3/2 St Sh M 2/3 2/3 5/6 St Sh M r 1/2 1/2 Sh M r

IV. Comparisons with numerical simulations and experiments

Here, we make comparisons between the predictions of the present theory and the results of discrete continuum numerical simulations of saltation of spheres over a rigid bumpy bed made of a layer of particles identical to those saltating in close contact. The simulations are performed in a quasi2D 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 upperbounded. We solve Newton's equations of motion for the individual spherical particles under the influence of fluid drag, buoyancy, gravity and contact forces when they collide with the bed (midtrajectory collisions are forbidden), while treating the fluid as either a viscous or a turbulent continuum, depending on the closure for the fluid shear stress, for which balance equations are phrased. Once changed in sign, the sum over all particles of drag and buoyancy enter the momentum balance for the fluid, thus ensuring the twoway coupling between the phases. In the numerical simulations, we assume that the fluid only possesses horizontal velocity and we suppress the possibility of interparticle collisions above the bed. The individual particles experience a fluid drag based upon the local difference between the instantaneous velocity of the particle and the average velocity of the fluid. As in our theoretical treatment, we control the amount of particles in the simulations -the particle mass holdup, M, the bumpiness of the rigid bed -the wallparticle diameter, dw, the viscosity -the inverse of the fall Stokes number, St, the mass density -the inverse of the density ratio, r, of the fluid, and the intensity of the shearing flow -the Shields number, Sh. In what follows, all measurements have been taken once a steady state is attained and, subsequently, timeaveraged. This type of numerical simulations have already been used in the context of saltation in both turbulent [START_REF] Durán | Numerical simulation of turbulent sediment transport, from bed load to saltation[END_REF][START_REF] Pähtz | The fluctuation energy balance in nonsuspended fluidmediated particle transport[END_REF]Pähtz and Durán 2020a;[START_REF] Ralaiarisoa | Transition from Saltation to Collisional Regime in Windblown Sand[END_REF]) and viscous [START_REF] Valance | Particle saltation over rigid bumpy beds in viscous shearing flows[END_REF] 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 extraterrestrial environment, as summarized in Table 2, by changing the mass holdup between 0.0004 and about 0.03 that is from the rarefied limit to the maximum mass hold up for which the midtrajectory collisions can be neglected and equation ( 23) is satisfied. Although the fluid regime on Mars and Venus is almost certainly turbulent, one can at least imagine performing experiments in pressurized windtunnels, 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.

The chosen flow conditions serve the purpose of isolating and testing the assumptions that we have made in building our theory: (i) anisotropic Maxwellian velocity distribution for the particles, vertical velocity of the ascending particles much larger than the settling velocity and locally linear velocity for the fluid in the case of viscous regime and Stokes drag -first row in Table 2; (ii) additional assumption of drag coefficient uniform and equal to that evaluated at the bed in the case of viscous regime and nonlinear drag -second and third row in Table 2; and (iii) fluid velocity profile encountered by an ascending particle uniform and equal to the average of the logarithmic profile in the case of turbulent regime and nonlinear drag -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 windtunnel [START_REF] Ho | Etude expérimentale du transport de particules dans une couche limite turbulente[END_REF]; for these, measurements of particle mass flux and profiles of particle concentration and horizontal particle and fluid velocities at different values of the particle mass holdup are available. In the experiments, unlike the numerical simulations, the vertical velocity of the fluid surrounding the particles is nonzero, due to the noslip condition on the particle surface, and this would permit a test of its influence on the results. However, as shown later, the depth of the saltation layer in the absence of an upper bound predicted by our model and measured in the numerical simulations is of the order of 4000 particle diameters. The experiments were performed in a rectangular closed conduit, with an horizontal lid placed at about 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 semianalytical treatment. Although the order of magnitude of the profiles measured in the experiments is in good agreement with our model, their shape reveals the influence of the upper boundary. We postpone to a future work the solution of the appropriate twopoint boundary value problem, with our proposed constitutive relations for the particle stresses and energy flux, and detailed comparisons against the experimental measurements. Table 2. Summary of the combination of parameters employed in the numerical simulations. Figures 5 and6 show the comparisons between profiles of Ty, c, y, s, u and U relative to selected values of the mass holdup (similar agreement, not shown here for brevity, is obtained for all admissible values of M) for saltation in viscous shearing flows with Stokes drag -first column of plots in both figures, saltation in viscous shearing flows with nonlinear drag -second column, and saltation in turbulent shearing flows with nonlinear drag -third column. In determining the semianalytical solution that we have highlighted in Section II, we have employed u = 0.6 and T = 0.5 for saltation in viscous flows, and u = 0.8 and T = 1 for saltation in turbulent flows. Different values of these parameters would not alter the qualitative features of our semianalytical solution, but only slightly impact the values of Ty, c and u at the bed, therefore causing a shift to the right or to the left in the relative profiles.

Figures 5(ac) confirm that the intensity of the velocity fluctuations of the particles decreases linearly with the distance from the rigid bed, and the model captures both the weak dependence of Ty on M for viscous saltation and Stokes drag (Figure 5a), and its independence for nonlinear drag (Figures 5b andc). 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 nonmonotonic behaviour of Ty 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 CD is allowed to vary locally; however, this prevents our obtaining semi analytical solutions of the governing equations, with little gain on the results. In all cases, the continuum model satisfactorily reproduces the profiles of particle concentration and normal stress in the ydirection (Figures 5di), which indeed follow the powerlaw distributions of equations ( 17) and ( 18). In the turbulent case, as anticipated, the particle concentration is rather uniform along y (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 only profiles of quantities associated with the vertical motion of the saltating particles, which is entirely uncoupled from the horizontal motion of the fluid in our continuum model and in the numerical simulations, were it not for 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.

Our proposed constitutive relations for the particle shear stress (equations 25 and 38) permit the qualitatively and quantitatively reproduction of the measurements in numerical simulations (Figures 6ac).

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 5060 % 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.

As the mass holdup increases, the fluid and particle horizontal velocity profiles become progressively nonlinear for saltation in viscous shearing flows (Figures 6d, e, g andh), but the continuum model is capable of capturing this behaviour. For saltation in turbulent shearing flows, the fluid velocity profile is logarithmic, and the mass holdup 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 a seemingly minor consequence of the lack of perfect agreement in terms of particle shear stress (Figure 6c). Finally, we compare the dependence of the particle mass flux on the particle mass holdup 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 windtunnel [START_REF] Ho | Etude expérimentale du transport de particules dans une couche limite turbulente[END_REF] are in sufficiently 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, that we do not account for. Figure 7a 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.

V. Conclusions

We have derived constitutive relations for the particle stresses and flux of particle kinetic energy associated with the fluctuating vertical velocities that apply to saltation of particles in viscous and turbulent shearing flows, in the absence of collisions above the bed. To do this, we have employed the averaging methods of statistical mechanics based on an anisotropic Maxwellian velocity distribution function for the ascending particles and approximate analytical expressions for the particle trajectories under the influence of fluid drag, gravity and buoyancy. Given that we neglect the possibility of particleparticle interactions, we have obtained the perhaps unexpected result that the particle shear stress does depend on the local relative velocity between the grains and the carried fluid and the fluid shear rate, but not on the particle shear rate. We have combined the constitutive relations for the particle phase and wellknown expressions for the fluid shear stress with momentum and energy balances to obtain semianalytical solutions of steady, fullydeveloped saltation over horizontal, rigid beds. We have employed boundary conditions appropriated for particles rebounding at a bumpy base with no upperbounds, and assumed that the carrier fluid has only mean horizontal velocity. Hence, the regime of the carrier fluid, viscous or turbulent, does not affect the profiles of particle concentration and normal stress, which follow a powerlaw decrease with height, and mean square of the vertical velocity fluctuations, linearly decreasing with the distance from the bed. These results are in contrast with previous theoretical and experimental works in which the granular temperature was assumed to be uniform and the concentration was shown to decrease exponentially with height.

We have confirmed our findings through extensive comparisons with discretecontinuum simulations in both viscous and turbulent regime. Thus, we are inclined to blame midfluid collisions for the above mentioned qualitative difference in the distribution of particle concentration. We have also shown that the predictions of the theory in terms of profiles of particle shear stress, and particle and fluid mean horizontal velocity are in excellent agreement with the numerical simulations. In natural units of particle mass density, diameter and reduced gravity, we have determined that the saltation process is controlled by the mass of particles in the system, the intensity of the shearing flow, the fluid mass density and viscosity and the bumpiness of the rigid base.

We have successfully tested the dependence of the horizontal mass flux per unit basal area of the particles as a function of the mass holdup obtained with our theory against discretecontinuum numerical simulations and windtunnel experiments, for a large range of the control parameters. We have also determined simple scaling laws for the dependence of various quantities of interest on the control variables in the special limits of rarefied viscous and turbulent saltation, and solely Stokes or form drag on the particles.

In a future, we plan to apply our theory to twopoint boundary value problems, to, e.g., mimic available experiments of saltation in enclosed windtunnels. We also wish to extend the present work to deal with unsteady and /or developing flows. More importantly, the inclusion of midfluid collisions, that is the transition from a collisionless to a collisional kinetic theory of saltation, is a crucial future step, especially in view of modelling transport phenomena over erodible beds, of more interest for geophysical and planetary science applications. Declaration of interests. The authors report no conflict of interest.

Appendix A. Approximate analytical trajectories for the saltating particles

We characterize the drag exerted on the particles through a drag coefficient, CD, that we take to be independent of y, as in [START_REF] Pasini | Aeolian transport with collisional suspension[END_REF], to obtain analytical expressions for the particle trajectories.

Integrating the vertical momentum balance for the particle:

d 1 0, d y D y C t + + =   (A1)
with y the vertical component of the particle velocity, we obtain the vertical velocity at any time t after the particle reaches a certain ascending velocity y  + at a certain location yr, ( ) ( )

1 exp ; exp D y D y D D C C t C C t   + + - = (A2)
and, with another integration, the position y as a function of time, ( )

( ) 2 1 1 exp . D y D D r D C C t C t y y C  + + - - -     -= (A3)
From equation (A3) we obtain an implicit expression for the time that the particle spends at y ≥ yr ( )

1 1 exp . D y f D f D C t C t C  + +   = - -   (A4)
Then, from equation (A2), with t = tf, the downward vertical velocity of the particle at y = yr is .

y y f t   - + = - (A5)
To permit subsequent analytical integrations, we propose the following explicit expression for the time of flight,

( ) ( )

1 1 exp , f y D y D D t C C C   + + = + - - (
1 1 exp . y D y D D C C C   - + = - + - (A7)
In the limit of large D y C  + , equation (A7) reduces to 1 .

y D C  -= - (A8)
We next integrate the horizontal particle momentum balances, and distinguish between saltation in viscous and in turbulent shearing flows.

A1. Saltation in viscous shearing flows

In the case of saltation in viscous shearing flows, we assume that the fluid velocity profile is locally linear (that is for y ≥ yr), so that the horizontal momentum balance for the particles is:

( ) d d , d d x D r x U C U y y t y     = + - -     (A9)
with x the horizontal component of the particle velocity, and U and d d

U y

the fluid horizontal velocity and shear rate at the reference level. In the following integrations, we treat d d

U y

for saltation in viscous shearing flows as if it were constant. Then, after inserting equation (A3) into equation (A9), and integrating with the

initial condition x(t = 0) = x + , ( ) ( ) ( ) ( ) ( ) ( )
2 2 1 d 1 exp exp d 1 exp d exp . d D y x D D D D D D x D D C U U C t C t C t y C C t C t U U C t y C    + + + = + - - - -     -+ - - + - - (A10)
The horizontal velocity after the time tf is, then, with equations (A4) and (A6), (

)

2 1 d 1 1 exp 1 exp exp d 1 exp , 1 x D y D y D y D y D y D D y x D y D y U U C C C C C C y C U C C           - + + + + + + + + +     = + + - - - - - -     + - + - + (A11)
In the limit of large D y C  + , equation (A11) gives

2 1 d . d x D U U C y  -= + (A12)
Equation ( A10) can be integrated with the boundary condition x(t = 0) = xi to obtain the horizontal displacement from the initial position at the reference level as a function of time:

( ) ( ) ( ) ( ) ( ) 2 2 2 1 d 2 2 exp d 1 d 1 1 1 exp 1 exp . d 2 D y i D D D D D D x D D D D D C U x x Ut C t t t C y C C C t U t C t U C t C y C C C   + + +     = + + - + + -           - -- - + + - - -         (A13)
The horizontal displacement  after the time tf is

( ) ( ) ( ) ( ) ( ) 2 2 2 1 d 2 2 exp d 1 d 1 1 1 exp 1 exp . d 2 D y f D f f f D D D D f f D f x D f D D D D C U Ut C t t t C y C C C t U t C t U C t C y C C C   + + +      = + - + + -             - -- - + + - - -         (A14)
In the limit of large D y

C  + , 1 f y D t C  + - = +
and equation (A14) reduces to

2 1 d . 2 d x y y D D U U U C y C    + + + -  = + + (A15)

A2. Saltation in turbulent shearing flows

In the case of saltation in turbulent shearing flows, we assume that the fluid velocity is equal to its average, U , in the region between the reference location yr and the top of the trajectory, yM, determined from equation (A3) at the time t for which y in equation (A2) vanishes. Then, the horizontal momentum balance for the particles reads:

( )

d d x D x C U t   = - (A16)
where, with the expression for the turbulent fluid shear stress, ( ) ( )

0 0 0 2 0 2 0 0 1 ln 1 d 1 ln 1 1 ln 1 1 1 ln 1 1 ln 1 1 0.85 ln M r y r M r r y r M r M r r D y D y r D D r D y D y y D r y y rS U U y y y y y y y y y rS U y y y y C C y y rS U C C y y C C rS U C y y          + + + + +   - = + +   - +         + - = + + + -       - +             - + +       = + + + -     + - +          + + +  d 1 0.85 ln 1 , d y D rS U U y C rS   +     = + +             (A17)
with U and d d U y the fluid horizontal velocity and shear rate at the reference level . In equation (A17), we have taken the limit of large D y C  + to approximate the average fluid velocity. The factor 0.85 in front of the logarithm allows to better fit the trajectories of discretecontinuum simulations of saltation in turbulent shearing flows.

Then, after integrating equation (A16) with the initial condition x(t = 0) = x + , ( ) ( )

1 exp exp . x D x D U C t C t   + = - - + -     (A18)
The horizontal velocity after the time tf is, then, with equations (A5) and (A7), ( ) ( )

exp 1 exp . 1 1 D y x x D y D y D y C U C C C       + + - + + +   -   = - + - + +     (A19)
In the limit of large D y C  + , equation (A19) gives

. x U  -= (A20)
Equation ( A18) can be integrated with the boundary condition x(t = 0) = xi to obtain the horizontal displacement from the initial position at the reference level as a function of time:

( )

1 1 exp . i x D D x x Ut U C t C  +   = + + - - -       (A21)
The horizontal displacement  after the time tf is ( )

1 1 exp . f x D f D Ut U C t C  +      = + - - -     (A22)
In the limit of large D y

C  + , 1 f y D t C  + - = +
and equation (A22) reduces to

.

x y D U C   + +  = + (A23) all all all f f f c c c        + - + + + -- - = = +    ξ ξ ξ ξ ξ ξ (B1)
where f is the velocity distribution of particles at that height and 6c/ is the number density. As in [START_REF] Creyssels | Saltating particles in a turbulent boundary layer: experiment and theory[END_REF], we distinguish between ascending or descending particles and introduce velocity distributions for both species. We assume that the velocity distribution of the ascending particles is an anisotropic Maxwellian:

( )

2 2 6 exp exp , 2 2 x y x y x y u c f T T T T     + + + + + + + +   -     = - -           (B2)
where:

u c f d T T T T c f T C T T c C c               - - + + + + + + + + -- - + + + + - + + + + + - - + + + + -   -     - - = -              - =        =      ξ ξ ξ ξ (B4)
that is a relation between the concentrations of the two species of particles. As shown later, we found that equation (B4) does not correctly capture the right limit as Ty + tends to zero, and also slightly underestimates c . An almost perfect agreement can be obtained, instead, with

3 . 4 y D T c c c C  + - + +  = + (B5)
From this and the definition of c,

3 2 . 4 y D T c c c c C  + + - +      = + = +     (B6)
Taking the limit for large D y C  + permits equations (B4)(B6) to be obtained without assuming a velocity distribution for the descending particles. We show the comparisons between equations (B5) and (B6) and the results of numerical simulations on saltation of particles in a viscous shearing flow experiencing Stokes drag in figure B1.

f f f f J c c c c f f J c             + - + + + + + + -- - + + + -- - + + + -- - + = + = +   = +        ξ ξ ξ ξ ξ ξ ξ ξ ξ ξ (B7)
where J is the Jacobian of the transformation from

+ ξ into - ξ . When
f f J f c c           + + + + + - -- - + + - + + +     = + = -       ξ ξ ξ ξ (B8)
where we made use of the fact that, at steady state, 0

u c C T T T T T c T c T c C T C T T                + + + + + + + + + + + + + - + + + + + + + + + + + + + +   -       = + - -                   = + - = +              (B10)
The limit for large D y C  + permits a simple expression for y (and for the other constitutive relations) to be obtained. This approximation is valid whenever the Maxwellian vertical velocity distribution of the ascending particles is wide, that is for large Ty + , so that the contribution of the large y The difference between the approximate equation (B10) and the more accurate equation ( B11) is shown in figure B2 when the drag coefficient is equal to 0.18, and we notice that the mismatch becomes relevant only at small Ty + , as expected. The agreement with the measurements in the numerical simulations of saltation of particles in a viscous shearing flow experiencing Stokes drag is remarkable.

u c C C C T T T T c T C C C T T T C T c T c T                 + + + + + + + + + + + + + + - + + + + + + + + + + + + + + + +   -       = + - - - -                   = + - - -             = +        2 erfc . 2 D y C T +           (B11)
Q f       + + + - + + + + - = -   (B14)
u c Q C T T T T T c T c T C T C T T               + + + + + + + + + + + + + - + + + + + + + + + + + + +   -       = - - -                     = - - = -                  (B15)
An undesirable, nonphysical consequence of equation ( B15) is that Qyyy vanishes at a finite, nonzero value of Ty + . However, in the same limit of large Figure B4 indicates that the predictions of equations (B15) and (B16) are in excellent agreement with the measurements in numerical simulations of saltation of particles in a viscous shearing flow experiencing Stokes drag, but only equation (B16) has the correct limit for vanishing Ty + .

The yyelement of the dissipation tensor due to fluid drag in the balance of particle second moment ( )

0 d d . 6 x x y x y s f       + + + - + + + + - = - -   (B18)
To proceed, we must distinguish between saltation in viscous and turbulent shearing flows, because of the influence of the fluid velocity on the particle horizontal velocity.

B1. Saltation in viscous shearing flows

In the limit of large D y

u c U s u u U T C dy T T T c U T U u d C y T T T T U c U u C y              + + + + + + + + + + + + + - + + + + + + + + + + +   -       = - - -+ -- -                   = -+ -             = -+        (B19)
Equation (B19) indicates that there is no influence of Tx + on the shear stress, because of the assumed Gaussian distribution of x + that causes the term involving (x + u) in the integrand of equation (B15) to disappear. However, [START_REF] Valance | Particle saltation over rigid bumpy beds in viscous shearing flows[END_REF] showed that the PDF of the downward velocity of the particles impacting the bed (dominated by the horizontal component) is actually nonsymmetric. Given the perfect agreement with respect to the other constitutive relations obtained in this Appendix, which do not involve the horizontal velocity of the ascending particles, we are inclined to ascribe to the nonsymmetric distribution of x + the fact that the shear stress measured in simulations of saltation particles in a viscous shearing flow experiencing Stokes drag (figure B6a) is not linear in the velocity difference, Uu, as implied by equation (B19). At small Ty + , equation (B10) in equation (B19) actually indicates that

2 1 d d D y D U s C U u C y     -+     . We found that 2 3 2 1 d 5 3 d D y D U s C U u C y    = -+     (B20)
permits to reproduce the behaviour of the particle shear stress for saltation in viscous shearing flows better than equation (B19) in almost the entire flow domain (figure B6b).

B2. Saltation in turbulent shearing flows

In the limit of large D y C  + , using equations (A17), (A20) and (B2) into equation (B14), and integrating, we obtain the following expression, in the case of saltation in turbulent shearing flows, f is the dimensionless contact force between particles p and q; I = /10 is the moment of inertia of a sphere; and , p q n is the unit vector along the contact direction.

          = - - -+ -- +                           -    -           = + + -                 
+ + + + + + + +     -                          -+ + -+ +                  (
The normal component fc,n of the contact force is modelled by a linear spring dashpot, so that fc,n = (knδ+γnvn), where kn is the spring stiffness, δ the overlap between the compliant spheres, γn the viscous damping coefficient and vn the normal component of the relative translational particle velocities. The negative of the ratio between the normal relative velocity before and after the collision is the coefficient of normal restitution en. If the values of en and kn are prescribed, γn is deduced from the following relation: 𝛾 𝑛 = (𝜋/6)√12𝑘 𝑛 /(1 + 𝜋 2 / ln(𝑒 𝑛 ) 2 ). The tangential component fc,t of the contact force is described via a Coulomb friction model regularized through a viscous damping: fc,t = -min(μfc,n, γtvt)sign (vt), where μ is the Coulomb friction coefficient, vt the relative slip velocity at contact and γt the tangential viscous damping coefficient. The values chosen for the parameters are: kn = p/6•10 7 , γn = γt, en = 0.88, and μ = 0.5.

The fluid motion is solved by an Eulerian description based on Eq. ( 24) or Eq. ( 37) for laminar and turbulent flows, respectively. The vertical component of the fluid velocity is assumed to be zero, so that only the horizontal momentum balance is required and reads where the infinite upper bound of the integral means that all moving grains located above y must be accounted for. Once the vertical profile S(y) of the fluid shear stress is determined, the horizontal fluid velocity profile can then be obtained from the integration Eq. ( 24) or Eq. ( 37) for laminar and turbulent shearing flows, respectively, with the noslip boundary condition U = 0 at y = 0.

The simulated system is quasitwodimensional with a streamwise length equal to 5120 particle diameters and a transverse length equal to one diameter. We use spherical particles with a polydispersity of ±10 % and adjust their number in the system to obtain the desired value of the mass holdup M. Periodic boundary conditions are employed in the streamwise direction. The domain is not upper bounded, while the lower boundary is composed of a layer of immobile particles of diameter dw in close contact (rigid, bumpy bed; see Figure 1). As mentioned earlier in the paper, we suppress the possibility of particle-particle collisions above the bed.

Operatively, at every time step, we integrate equations (D1) and (D2) for every particle in the system and determine its new velocity and position. We then use this information to update the profiles of fluid shear stress and horizontal velocity via equations (4), and (24) for viscous flows or equation (37) for turbulent flows and proceed with the next time step until we reach a steady state, that is, when the horizontal mass flux averaged over a window of 100 unit time steps is stationary. Initially, the fluid profile is taken to be linear or logarithmic for viscous and turbulent shearing flows, respectively, and corresponds to the unperturbed profile in the absence of particles. The particles are initially displayed on a horizontal row located at ten diameters from the rigid, bumpy bed, with a constant interparticle distance equal to two diamters and zero initial velocity. Importantly, the final state is independent of the initial conditions as long as the number of particles in the flow does not surpass its transport capacity.

Fig. 1 .

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

Fig. 2 .

 2 Fig.2. Normalized profiles of Ty + (dotdashed line) and Ty (dashed line) obtained from equation (14) and (11), respectively, when CD = 0.18. The solid line represents the approximate linear distribution of Ty (equation 16).

Fig. 3 .

 3 Fig. 3. Predicted impact angle (in degrees, solid line), drag coefficient (dashed line) and particle shear stress at the bed (dotdashed line) as functions of the mass holdup for saltation in viscous shearing flows when dw = 1, St = 100, r = 50, and Sh = 0.05 (with u = 0.6).Using equation (25) in equation (2), with equations (6), (18) and (24), gives a first order, linear, non homogeneous, differential equation for the particle shear stress,

Fig. 4 .

 4 Fig. 4. Predicted impact angle (in degrees, solid line), drag coefficient (dashed line) and particle shear stress at the bed (dotdashed line) as functions of the mass holdup for saltation in turbulent shearing flows with dw = 1.5, St = 1681, r = 2208, and Sh = 0.04 (with u = 0.8 and T= 1).Equation (2), with equations (38) and (6), and equation (37) can be written as a system of two ordinary differential equations:

0

  Sh s, when there is only Stokes drag -18 St D C = , when the particle Reynolds number based upon the relative velocity between the particles and the fluid at the bed, 0 St u r , is less than unity, or form drag -

Fig. 5 .

 5 Fig. 5. Profiles of (ac) mean square of particle velocity fluctuations in the vertical direction, (df) particle concentration and (gi) particle normal stress along y measured in numerical simulations of: saltation in viscous shearing flows with dw = 1, St = 100, r = 150000, Sh = 0.05 and M = 0.0008 (blue circles), M = 0.0033 (red circles), M = 0.0131 (purple circles); saltation in viscous shearing flows with dw = 1, St = 100, r = 50, Sh = 0.05 and M = 0.0008 (blue squares), M = 0.0033 (squares), M = 0.0131 (purple squares); saltation in turbulent shearing flows with dw = 1.5, St = 1681, r = 2208, Sh = 0.04 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.

Fig. 6 .

 6 Fig. 6. Profiles of (ac) particle shear stress, (df) fluid and (gi) particle mean horizontal velocities measured in numerical simulations of: saltation in viscous shearing flows with dw = 1, St = 100, r = 150000, Sh = 0.05 and M = 0.0008 (blue circles), M = 0.0033 (red circles), M = 0.0131 (purple circles); saltation in viscous shearing flows with dw = 1, St = 100, r = 50, Sh = 0.05 and M = 0.0008 (blue squares), M = 0.0033 (squares), M = 0.0131 (purple squares); saltation in turbulent shearing flows with dw = 1.5, St = 1681, r = 2208, Sh = 0.04 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.

Fig. 7 .

 7 Fig. 7. (a) Particle mass flux against scaled particle mass holdup measured in discretecontinuum numerical simulations (open symbols) and predicted from the present theory (lines) in the case of saltation in viscous shearing flows with dw = 1, r = 150000 and: St = 60 and Sh = 0.05 (blue circles); St = 100 and Sh = 0.05 (orange circles); St = 100 and

  A6) which has the right limits at both small D y C  + , i.e.,

  mean squares of the horizontal velocity fluctuations of the ascending particles. The concentration of the ascending particles, c + , is equal to take the limit for large D y C  + of y  -(equation A9) and, inserting equation (B2) into equation (

Fig. B1 .

 B1 Fig. B1. Concentration ratios c/c + (circles) and c /c + (squares) as functions of Ty + measured in numerical simulations of saltation in a viscous shearing flow with dw = 1, St = 100, r = 150000, Sh = 0.05 (so that the Stokes drag is dominant and CD = 18/St = 0.18) and M = 0.0131. The dashed and solid lines are the predictions of equations (B5) and (B6), respectively.

y

  = and, therefore, from equation (B3), if the trajectories are independent of each other (that is for small mass holdups), stress in the ydirection, in the absence of particle interaction, is simply given by the average vertical flux of ymomentum, taking the limit for large D y C  + of y  -(equation A9) and, inserting equation (B2) into equation (B9), and integrating, we obtain



  + is significant. Conversely, it becomes less accurate near the top of the saltating layer, where we expect small values of Ty + and a subsequent narrower distribution of y  + around zero. If we use the more accurate equation (A7) in equation (B9), we would obtain the following constitutive relation for the particle normal stress:

Fig. B2 .

 B2 Fig. B2. Ratio of particle normal stress over particle concentration of ascending particles as a function of Ty + measured in numerical simulations of saltation in a viscous shearing flow with dw = 1, St = 100, r = 150000, Sh = 0.05 (so that the Stokes drag is dominant and CD = 18/St = 0.18) and M = 0.0131. The solid and dashed lines are the predictions of equations (B10) and (B11), respectively.

  agreement between equations (B12) and (B13) and the measurements in the numerical simulations of saltation of particles in a viscous shearing flow experiencing Stokes drag in figureB3.

Fig. B3 .

 B3 Fig. B3. Ty (circles) and Ty (squares) as functions of Ty + measured in numerical simulations of saltation in a viscous shearing flow dw = 1, St = 100, r = 150000, Sh = 0.05 (so that the Stokes drag is dominant and CD = 18/St = 0.18) and M = 0.0131. The solid and dashed lines are the predictions of equations (B12) and (B13), respectively.

  Fig. B4. Ratio of energy flux Qyyy over particle concentration of ascending particles as a function of Ty + measured in numerical simulations of saltation in a viscous shearing flow with dw = 1, St = 100, r = 150000, Sh = 0.05 (so that the Stokes drag is dominant and CD = 18/St = 0.18) and M = 0.0131. The dashed and solid lines are the predictions of equations (B15) and (B16), respectively.

Fig. B5 .

 B5 Fig. B5. Negative of the derivative of Qyyy over particle concentration of ascending particles as a function of Ty + measured in numerical simulations of saltation in a viscous shearing flow with dw = 1, St = 100, r = 150000, Sh = 0.05 (so that the Stokes drag is dominant and CD = 18/St = 0.18) and M = 0.0131. The solid line is the predictions of Dyy/c + from equations (B17).

  Alam 2017), so that, with the definition of y and equation (B10), indicates that Dyy should be equal to the negative of the derivative of Qyyy. simulations of saltation of particles in a viscous shearing flow experiencing Stokes drag as a function of Ty + do indeed match the predictions of / yy D c + from equation (B17).Finally, the particle shear stress is equal to the negative of the average vertical flux of xmomentum,

Fig. B6 .

 B6 Fig.B6. Ratio of (a) particle shear stress s and the product of particle concentration of ascending particles and the square root of Ty + and (b) particle shear stress and the particle normal stress as functions of the velocity difference, Uu, measured in numerical simulations of saltation in a viscous shearing flow with dw = 1, St = 100, r = 150000, Sh = 0.05 (so that the Stokes drag is dominant and CD = 18/St = 0.18) and M = 0.0004. The solid lines are the predictions of equations (B19) and (B20), respectively.

  of the logarithm to carry out the analytical integration, and, we have kept the same dependency on the drag coefficient and the normal stress as in the case of saltation in viscous shearing flows (equation B20). The factor 3/4 in front of the logarithm allows to better fit the results of discretecontinuum simulations of saltation in turbulent shearing flows. of proportionality T  weakly dependent on the coefficients of normal and tangential restitution and, perhaps, on the flow regime of the fluid.[START_REF] Lämmel | Grainscale modeling and splash parametrization for aeolian sand transport[END_REF] also obtained

  = 𝑐 < ∑ 𝐶 𝐷 (𝑈 -𝜉 𝑥 ) 𝑝𝜖[𝑦; 𝑦+𝑑𝑦] >/ < ∑ 1 𝑝𝜖[𝑦; 𝑦+𝑑𝑦] >,with the angular brackets denoting ensemble averaging, represents the xcomponent of the average volume force exerted on the fluid by the particles whose centers are located in the horizontal slice between y and y+dy. The integration of equation

Table 1 .

 1 Summary of the scaling laws for rarefied saltation.

		Viscous saltation	Viscous saltation	Turbulent saltation
		with Stokes drag	with form drag	with form drag
	Drag coefficient, CD	1

  and p ω are the translational and angular velocity vectors of particle p, respectively; x e and y e are the horizontal and vertical unit vectors, respectively;

									ω					n					f	(D2)
	where p ξ ( p drag D C U   = - f	p x	)	e	x	-		p y y e	 	is the
	dimensionless drag force exerted by the fluid on the pparticle and , p q c
			( ) ( ) ( ) ( ) 2 ,0 ,0 2 ,0 ,0 y x y x     + + --+ +	2 2		=	   	t e	-	(	t e	2	2 e t e n -	2	)		2 1 + d d w w	   	2	,	(C6)
	so that, with equations (A8) and (A12), and ( ) ,0 x  +	2					,0 + x		2	=	( ) 0 u +	2	,
	u	0 +		1 D C	2	  	1	+	2 1 	      	t e	-	(	2 e e 2 t n t e -	2	)		2 1 + d d w w	   	2	-	,0 + y T	.

(C7)

Then, we can determine the total particle slip velocity at the bed as =

Journal of Fluid Mechanics

Appendix B. Constitutive relations

We define the average of a certain quantity  at a certain distance y from the bed as: Appendix C. Boundary conditions for shallow impacts at a rigid, bumpy bed

After [START_REF] Lämmel | Grainscale modeling and splash parametrization for aeolian sand transport[END_REF], in the case of angle of impact  much less than /2 (shallow impact) of particles with large coefficient of friction, the negative of the average ratio of the vertical particle velocity after and before the impact is

where en and et are the coefficients of normal and tangential restitution, respectively, and, here and in what follows, we use the subscript 0 to indicate quantities evaluated at the bed.

In the limit of large takeoff velocity, equations (A8), (A12) and ( 24) give

for saltation in viscous shearing flows.

In the case of saltation in turbulent shearing flows, instead, in the limit of large takeoff velocity, equations (A8), (A17) and (A20) give

With equations (B5) and (B6), and 0 u -given by equation (A12), we obtain, at leading order, The discretecontinuum (DC) simulations is based on the combination of a discrete element method for the particle dynamics coupled to a continuum description of hydrodynamics, as developed in [START_REF] Durán | Numerical simulation of turbulent sediment transport, from bed load to saltation[END_REF], [START_REF] Ralaiarisoa | Transition from Saltation to Collisional Regime in Windblown Sand[END_REF] and [START_REF] Valance | Particle saltation over rigid bumpy beds in viscous shearing flows[END_REF]. The particle motion is described by a Lagrangian approach according to which the particle labelled p obeys the following dimensionless equations: