Universal friction law at granular solid-gas transition explains scaling of sediment transport load with excess fluid shear stress
UUniversal friction law at granular solid-gas transition explains scaling of sedimenttransport load with excess fluid shear stress
Thomas P¨ahtz , ∗ and Orencio Dur´an
1. Institute of Port, Coastal and Offshore Engineering,Ocean College, Zhejiang University, 310058 Hangzhou, China2. State Key Laboratory of Satellite Ocean Environment Dynamics,Second Institute of Oceanography, 310012 Hangzhou, China3. Department of Ocean Engineering, Texas A&M University, College Station, Texas 77843-3136, USA
A key interest in geomorphology is to predict how the shear stress τ exerted by a turbulent flowof air or liquid onto an erodible sediment bed affects the transport load M ˜ g (i.e., the submergedweight of transported nonsuspended sediment per unit area) and its average velocity when exceedingthe sediment transport threshold τ t . Most transport rate predictions in the literature are based onthe scaling M ˜ g ∝ τ − τ t , the physical origin of which, however, has remained controversial. Herewe test the universality and study the origin of this scaling law using particle-scale simulationsof nonsuspended sediment transport driven by a large range of Newtonian fluids. We find thatthe scaling coefficient is a universal approximate constant and can be understood as an inversegranular friction coefficient (i.e., the ratio between granular shear stress and normal-bed pressure)evaluated at the base of the transport layer (i.e., the effective elevation of energetic particle-bedrebounds). Usually, the granular flow at this base is gaslike and rapidly turns into the solidlikegranular bed underneath: a liquidlike regime does not necessarily exist, which is accentuated bya nonlocal granular flow rheology in both the transport layer and bed. Hence, this transitionfundamentally differs from the solid-liquid transition (i.e., yielding) in dense granular flows eventhough both transitions are described by a friction law. Combining this result with recent insightsinto the nature of τ t , we conclude that the transport load scaling is a signature of a steady reboundstate and unrelated to entrainment of bed sediment. PACS numbers: 45.70.-n, 47.55.Kf, 92.40.Gc ∗ a r X i v : . [ c ond - m a t . s o f t ] O c t I. INTRODUCTION
The transport of sediment mediated by the turbulent shearing flow of a Newtonian fluid over an erodible granular bedis responsible for the evolution of fluid-sheared surfaces composed of loose sediment, such as river and ocean beds, andwind-blown sand surfaces on Earth and other planets, provided that the sediment is not kept suspended by the fluidturbulence [1–14]. Nonsuspended sediment transport thus constitutes one of the most important geomorphologicalprocesses in which granular particles collectively move like a continuum flow, and predicting the associated sedimenttransport rate Q (i.e., the total particle momentum in the flow direction per unit bed area) and flow threshold τ t (i.e.,the value of the fluid shear stress τ below which sediment transport ceases) are considered central problems in Earthand planetary geomorphology [1–14]. Here we provide the theoretical base necessary to understand the scaling of Q and τ t and, by doing so, show that and why nonsuspended sediment transport constitutes a class of granular flowswith unique properties, such as a nonlocal granular flow rheology even relatively far from the flow threshold. A. The scaling of the transport rate of nonsuspended sediment
Numerous experimental and theoretical studies (e.g., Refs. [15–67]) have measured or derived analytical expressionsfor the transport rate Q as a function of particle and environmental parameters, such as the particle (fluid) density ρ p ( ρ f ), kinematic fluid viscosity ν , characteristic particle diameter d , gravitational constant g , and τ and τ t . Mostof the theoretical derivations are based on, or can be reformulated in the spirit of, Bagnold’s [18–20] pioneering ideas.Defining a Cartesian coordinate system x = ( x, y, z ), where x is in the flow direction, z in the direction normal to thebed oriented upwards, and y in the lateral direction, Bagnold assumed that there is a well-defined interface z = z r between granular bed ( z < z r ) and transport layer ( z > z r ), which we henceforth call the “Bagnold interface,” withthe following properties (Fig. 1):1. The transport rate Q r above z r well approximates the total transport rate Q (i.e., z r cannot be too far away fromthe actual granular bed). Hence, one can separate Q into the mass M = ρ p (cid:82) ∞ z r φ d z of particles located above z r per unit bed area, where φ is the particle volume fraction (i.e., the fraction of space covered by particles), andthe average horizontal velocity v x with which particles located above z r move: v x ≡ Q r /M (cid:39) Q/M .2. The ratio µ ≡ − P zx /P zz between the particle shear stress − P zx and normal-bed pressure P zz , where P ij is theparticle stress tensor, at z r does not significantly depend on the fluid shear stress τ : µ b ≡ µ ( z r ) (cid:54) = f ( τ ).3. The ratio − P zx ( z r ) /τ between particle and fluid shear stress increases from nearly zero at low transport stages( τ /τ t − (cid:28)
1) to nearly unity at large transport stages ( τ /τ t − (cid:29) − P zx ( z r ) = τ − τ t and − P zx ( z r ) = √ τ ( √ τ − √ τ t ). Note that the former expression is usuallyattributed to Owen [55] (“Owen’s second hypothesis” [68]) in the aeolian transport literature [10–12] even thoughBagnold [18] was its originator and also applied it to aeolian transport. Q > z / Q -4-2024 E l e v a t i o n , ( z − z r ) / d µ -4-2024 0 0.4 0.8 1.2 − P zx / τ -4-2024 (c)(b)(a) Arrows : Increasing τ FIG. 1.
Visualization of Bagnold interface properties.
Vertical profiles of (a) the fraction Q >z /Q of sediment transportoccurring above elevation z , (b) the friction coefficient µ , and (c) the ratio − P zx /τ between the particle shear stress − P zx andfluid shear stress τ . The solid lines correspond to data obtained from direct sediment transport simulations (see Sec. II) fortwo representative cases: turbulent bedload (turquoise) and saltation transport (brown). The black, dashed lines mark theBagnold interface z = z r . Combining these three properties and using the vertical momentum balance P (cid:48) zz (cid:39) − ρ p φ ˜ g of steady, homogeneoussediment transport [69], where the prime denotes the derivative d / d z and ˜ g = (1 − ρ f /ρ p ) g the buoyancy-reducedvalue of g , then yields Q (cid:39) µ − b ˜ g − ( τ − τ t ) v x if − P zx ( z r ) = τ − τ t ,Q (cid:39) µ − b ˜ g − √ τ ( √ τ − √ τ t ) v x if − P zx ( z r ) = √ τ ( √ τ − √ τ t ) . (1)Indeed, the functional behaviors in Eq. (1) resemble the vast majority of theoretical and experimental threshold shearstress-based expressions for the transport load M ˜ g (cid:39) Q ˜ g/v x and transport rate Q in the literature, which differ onlyin their prediction of v x . For example, experiments of nonsuspended sediment transport driven by turbulent streamsof liquid (turbulent “bedload”) suggest that v x is linear in (cid:112) τ /ρ f [29–32], whereas experiments of nonsuspendedsediment transport driven by turbulent streams of air (turbulent “saltation”) suggest that v x is constant with τ [48–50]. The capability of Eq. (1) to reproduce experimental data is indirect evidence that the Bagnold interfaceexists for these conditions. However, there are a number of unsolved problems, even inconsistencies, regarding thegenerality and physical origin of the Bagnold interface that currently prevent us from understanding and predictingthe scaling laws of nonsuspended sediment transport for arbitrary conditions and from integrating nonsuspendedsediment transport within the framework of granular flow rheology. B. Open questions
1. Existence of the Bagnold interface
Natural granular beds are locally very heterogeneous and undergo continuous rearrangements during sedimenttransport, which renders the definition of a bed-transport-layer interface difficult. For steady, homogeneous transportconditions, four different definitions have been proposed in the literature: the elevation at which the friction coefficient µ exhibits a certain constant value [70], the elevation at which the particle volume fraction φ exhibits a certain constantportion of the bed packing fraction φ b [51], the elevation at which the particle shear rate ˙ γ exhibits a certain constantportion of its maximal value [32], and the elevation at which the production rate P zz ˙ γ of cross-correlation fluctuationenergy is maximal [71, 72]. However, whether any of these interfaces is the Bagnold interface and whether the Bagnoldinterface even exists for nonsuspended sediment transport in arbitrary environments remain unclear.In this study, we provide answers to the following questions: • Does the Bagnold interface exist in general settings? • If so, is there a general definition of the Bagnold interface?
2. Physical origin of friction law
Property 2 of the Bagnold interface represents a macroscopic, dynamic friction law, analogous to Coulomb frictiondescribing the sliding of an object down an inclined plane, where the constant dynamic bed friction coefficient µ b is the analog to the ratio between the horizontal and normal force acting on the sliding object. In the context ofdense ( φ (cid:38) .
4) granular flows and suspensions, it is well established that a constant dynamic friction coefficient (theyield stress ratio) characterizes the transition between solidlike and liquidlike flow behavior [73–92]. Here liquidlikebehavior refers to dense flows that obey a local rheology (i.e., µ depends only on a single local quantity, such as φ ),while solidlike behavior refers to both quasistatic and creeping flows (not to be confused with Bagnold’s term “surfacecreep” [1]). Quasistatic flows are associated with very small, reversible deformations of dense packed granular systems,while creeping flows are associated with an exponential relaxation of the particle shear rate ˙ γ between quasistaticand liquidlike flows [90–96] and characterized by a nonlocal granular flow rheology [90–92]. Based on the fact thata friction law characterizes the solid-liquid transition, it has been very common to argue that the Bagnold interfaceseparates a solidlike granular bed from a liquidlike transport layer on its top and that µ b is the yield stress ratio[21–28], which is in the spirit of Bagnold’s original reasoning [18–20]. However, this interpretation is inconsistent withProperty 3 of the Bagnold interface, which predicts that the particle shear stress − P zx ( z r ), and thus the particlevolume fraction φ ( z r ) [69], becomes very small when the fluid shear stress approaches the flow threshold ( τ → τ t ). Itis further inconsistent with the fact that the Bagnold interface is also found in highly simplified numerical sedimenttransport simulations that do not resolve particle interactions [25, 62].An alternative interpretation of the friction law came from studies on saltation transport [60–65, 97]. They suggestedthat µ b is an effective restitution coefficient characterizing an approximately constant ratio between the averagehorizontal momentum loss and vertical momentum gain of particles rebounding at the Bagnold interface. However,this interpretation has never been tested against experiments or numerical particle-scale simulations of sedimenttransport, and it is unclear how it can be generalized to the bedload transport regime, in which transported particlesexperience long-lasting contacts with the granular bed and each other [98].In this study, we provide answers to the following questions: • What is the physical origin of the friction law at the Bagnold interface? • Is this origin in some way associated with the rheology of dense granular flows and suspensions?
3. Universality of friction law
For the purpose of understanding the scaling laws of nonsuspended sediment transport in arbitrary environments, itis crucial to know how much the dynamic bed friction coefficient µ b at the Bagnold interface varies with environmentalparameters other than τ . Currently, the literature suggests that the friction coefficient µ at elevations near the bedsurface, and thus near the Bagnold interface, strongly depends on the fluid driving transport (reported values rangefrom 0 . . µ have been limited to systems that only crudely represent natural nonsuspended sediment transport, such as themotion of externally fed particles along rigid beds [36, 99, 100], or µ has been estimated as τ /P zz [101], which makessense only for intense transport conditions due to Property 3.In this study we provide an answer to the following question: • How much does the dynamic friction coefficient µ b at the Bagnold interface vary with environmental parameters? C. Organization of this paper
The method that we use to answer the open questions outlined above, direct sediment transport simulations withthe model of Ref. [51], is briefly introduced in Sec. II. Section III then puts forward our definition of the bed-transport-layer interface as the effective elevation at which the most energetic transported particles rebound when colliding withbed surface particles and shows that this interface is the Bagnold interface. It also shows that the friction law atthe Bagnold interface is, indeed, universal. Section IV links this finding, for the vast majority of sediment transportregimes, to a steady transport state in which transported particles continuously rebound at the bed surface and showsthat alternative explanations associated with the rheology of dense granular flows and suspensions in general fail dueto the absence of a liquidlike flow regime. Finally, Sec. V summarizes the main conclusions that can be drawn fromour results and discusses our results in the context of sediment transport modeling.
II. NUMERICAL SIMULATIONS
In this section, we describe the numerical model (Sec. II A), the simulated sediment transport conditions (Sec. II B),and how we use the simulation data to compute relevant physical quantities (Sec. II C).
A. Numerical model description
The numerical model of sediment transport in a Newtonian fluid of Ref. [51] belongs to a new generation ofsophisticated grain-scale models of sediment transport [11, 51, 69–72, 85, 98, 102–125] and has been shown to reproducemany observations concerning viscous and turbulent nonsuspended sediment transport in air and water [11, 51, 71,72, 105], and bedform formation [106]. It couples a discrete element method for the particle motion with a continuumReynolds-averaged description of hydrodynamics, which means that it neglects turbulent fluctuations around the meanturbulent flow. It simulates the translational and rotational dynamics of ≈ ,
000 spheres, including >
10 layers ofbed particles (more than sufficient to completely dissipate the energy of particles impacting the bed surface), withdiameters d p evenly distributed within two sizes (0 . d and 1 . d ) in a quasi-2-D, vertically infinite domain of length1181 d . Periodic boundary conditions are imposed along the flow direction, while the bottommost layer of particlesis glued to a bottom wall. The particle contact model considers normal repulsion (restitution coefficient e ), energydissipation, and tangential friction, where the magnitude of the tangential friction force relative to the normal contactforce is limited through a Coulomb friction criterion (contact friction coefficient µ c = 0 . e (Sec. II B). We refer thereader to the original publication [51] for further details (note that we recently corrected slight inaccuracies in theoriginal model [71]). B. Simulated sediment transport conditions
Using the numerical model, we simulate steady, homogeneous sediment transport for a particle-fluid-density ratio s ≡ ρ p /ρ f within the range s ∈ [1 . , ≡ (cid:112) ( s − gd /ν within the range Ga ∈ [0 . , e = 0 .
9. For small density ratio ( s ≤ . e = 0 .
01 because e can become very small for small Stokes numbers due to lubrication forces [126–128]. For each set of s , Ga, and e , we vary the dimensionless fluid shear stress (“Shields number”) Θ = τ / [( ρ p − ρ f ) gd ]in regular intervals above its threshold value Θ t = τ t / [( ρ p − ρ f ) gd ], which we obtain from extrapolation to vanishingtransport [72]. The simulated conditions cover four major, and very distinct, natural transport regimes, which dependon the transport layer thickness and the thickness of the viscous sublayer of the turbulent boundary layer [72]: viscousbedload transport, such as the transport of sand by oil; turbulent bedload transport, such as the transport of gravel bywater; viscous saltation transport, such as the transport of sand by wind on Mars; and turbulent saltation transport,such as the transport of sand by wind on Earth. They also cover 5 orders of magnitude of the ‘impact number’Im ≡ √ s + 0 . (cid:39) √ s Ga, which characterizes the mode of entrainment of bed sediment under threshold conditions[71]: Im (cid:38)
20 when entrainment by particle-bed impacts dominates entrainment by the mean turbulent flow, Im (cid:46) (cid:46) Im (cid:46) C. Computation of local averages and particle stresses
We use the simulation data to compute local averages of particle properties and the particle stress tensor, which isexplained in the following.
1. Local, mass-weighted time average and particle volume fraction
We compute the local, mass-weighted time average (cid:104) A (cid:105) of a particle quantity A through [69] (cid:104) A (cid:105) = 1∆ φ (cid:88) n V np A n δ ( z − z n ) T , (2) φ = 1∆ (cid:88) n V np δ ( z − z n ) T , (3)where ∆ = 1181 d is the simulation area, φ is the local particle volume fraction, z n ( V np = πd n p /
6) is the elevation(volume) of particle n , δ the δ distribution, and · T = T (cid:82) T · d t denotes the time average over a sufficiently long time T .The δ kernels have been coarse grained through spatial averaging over a discretization box of size 1181 d × d × ∆ z , where∆ z varies between 0 . d in dense and dilute flow regions ( φ (cid:38) .
1) and larger values in rarefied regions. Henceforth,the δ symbol should thus be interpreted as the associated coarse-graining function.
2. Particle stress tensor
The particle stress tensor P ij is composed of a kinetic contribution due to the transport of momentum betweencontacts (superscript ‘t’) and a contact contribution (superscript ‘c’) and computed through [69] P ij = P tij + P cij , (4a) P tij = ρ p φ (cid:104) c i c j (cid:105) , (4b) P cij = 12∆ (cid:88) mn F mnj ( x mi − x ni ) K ( z, z m , z n ) T , (4c)where K = (cid:82) δ { z − [( z m − z n )˜ s + z n ] } d˜ s , c = v − (cid:104) v (cid:105) is the fluctuation velocity, and F mn the contact force applied byparticle n on particle m ( F mm = 0). We confirmed that these definitions are consistent with the steady momentumbalance P (cid:48) zi = ρ p φ (cid:104) a i (cid:105) [69], where a is the particle acceleration due to noncontact forces. III. EXISTENCE OF THE BAGNOLD INTERFACE IN ARBITRARY ENVIRONMENTS
In Sec. III A, we first put forward our definition of the bed-transport-layer interface. In Sec. III B, we then showwith data from our direct transport simulations that this definition, in contrast to common alternative definitions,obeys the properties of the Bagnold interface (except for a slight restriction regarding Property 3) with a universallyapproximately constant bed friction coefficient µ b . A. Definition of the bed-transport-layer interface
In order to motivate a definition of the bed-transport-layer interface that results in the Bagnold interface, we exploitthe fact that numerical studies that represent the granular bed surface by a rigid bottom wall found that this wallobeys Properties 1-3 of the Bagnold interface [36, 62]. This finding suggests that an appropriate definition should havecharacteristics that mimic those of particle rebounds at rigid boundaries. One such characteristic is the productionof particle velocity fluctuations. For example, gravity-driven granular flows down an inclined, rigid base exhibit amaximum of the granular temperature (cid:104) c (cid:105) near this base [129]. The probable reason is that such rigid boundariesinduce strong correlations between the velocities of descending particles before rebound and ascending particles afterrebound.In steady sediment transport, the mass balance dictates (cid:104) v z (cid:105) = 0 [69], which can be achieved only if rebounds oftransported particles at the granular bed partially convert horizontal momentum of descending particles into verticalmomentum of ascending particles (i.e., negative correlation). Similar to gravity-driven granular flows, this constraintimplies that particle-bed rebounds are a strong source of the negative cross-correlation fluctuation energy density − ρ p φ (cid:104) c z c x (cid:105) .The balances of − ρ p φ (cid:104) c z c x (cid:105) and of the actual fluctuation energy density ρ p φ (cid:104) c (cid:105) can be derived rigorously fromNewton’s axioms. For steady sediment transport ( ∂/∂x = ∂/∂y = ∂/∂t = 0), they read [69] (Einsteinian summation) − q (cid:48) z ( xz ) = 12 P zz ˙ γ + Γ drag( xz ) + Γ coll( xz ) , (5a) q (cid:48) zii = − P zx ˙ γ − Γ drag ii − Γ coll ii , (5b)respectively, where the parentheses denote the symmetrization in the indices [ A ( ij ) = ( A ij + A ji )]. Furthermore, q ijk = ρ p φ (cid:104) c i c j c k (cid:105) + (cid:80) mn F mnj c k ( x mi − x ni ) K ( z, z m , z n ) T is the flux tensor of fluctuation energy, ˙ γ = (cid:104) v x (cid:105) (cid:48) the par-ticle shear rate, Γ drag ij = − ρ p φ (cid:104) a i c j (cid:105) the drag dissipation rate tensor, and Γ coll ij = − (cid:80) mn F mni ( v mj − v nj ) δ ( z − z m ) T the collisional dissipation rate tensor. In Eq. (5b), − P zx ˙ γ corresponds to the production rate and Γ drag ii and Γ coll ii to the dissipation rate of ρ p φ (cid:104) c (cid:105) by fluid drag and collisions, respectively. In Eq. (5a), P zz ˙ γ corresponds to theproduction rate and − Γ drag( xz ) and − Γ coll( xz ) to the dissipation rate of − ρ p φ (cid:104) c z c x (cid:105) by fluid drag and collisions, respectively.Hence, if we identify the bed-transport-layer interface as the average elevation of energetic particle-bed rebounds anduse that such rebounds are a strong source of − ρ p φ (cid:104) c z c x (cid:105) , it makes sense to define this interface through a maximumof the local production rate of − ρ p φ (cid:104) c z c x (cid:105) : max( P zz ˙ γ ) = [ P zz ˙ γ ]( z r ) , (6)which is exactly the definition that we applied in two recent studies [71, 72]. Figures 2(a) and 2(b) show exemplaryvertical profiles relative to z r of P zz ˙ γ for (a) weak and (b) intense viscous and turbulent bedload transport andturbulent saltation transport, where the bedload cases have been simulated using two different restitution coefficientsto mimic the minimal ( e = 0 .
9) and nearly maximal ( e = 0 .
01) effect that lubrication forces can possibly have. Itcan be seen that the value of e does not significantly affect these profiles. As we will see later, the influence of e onbedload transport properties is very small in general, consistent with previous studies [70–72, 122, 130].The interface z = z r defined by Eq. (6) shares some similarities with the region in which the production rate offluctuation energy is nearly balanced by the collisional energy dissipation rate: − P zx ˙ γ (cid:39) Γ coll ii . For turbulent bedloadtransport, it has been speculated that this region is a distinct granular layer (the “dense algebraic layer”) with athickness of several particle diameters d and that the bottom of this layer corresponds to the bed-transport-layerinterface [39, 40]. However, Figs. 2(c) and 2(d) show for the same cases as before that the thickness of the region in -5 -4 -3 -2 -1 P zz ˙ γ / [ ρ p ˜g √ ˜gd ] -4-2024 ( z − z r ) / d -4 -3 -2 -1 P zz ˙ γ / [ ρ p ˜g √ ˜gd ] -4-2024 s = . , Ga = . , e = . = . , Ga = . , e = . = . , Ga = , e = . = . , Ga = , e = . -3 -2 -1 − P zx ˙ γ / Γ collii -4-2024 ( z − z r ) / d -3 -2 -1 − P zx ˙ γ / Γ collii -4-2024 s = , Ga = , e = . (a)(c) Intense transportWeak transportWeak transport Intense transport (d)(b)
FIG. 2.
Exemplary vertical profiles of quantities associated with the (cross-correlation) fluctuation energybalance. (a, b) Vertical profiles relative to the rebound location z r of P zz ˙ γ/ ( ρ p ˜ g √ ˜ gd ). (c, d) Vertical profiles of − P zx ˙ γ/ Γ coll ii .Symbols correspond to viscous bedload transport [ s = 2 .
65, Ga = 0 .
1, Θ / Θ t = (1 . , . s = 2 . / Θ t = (2 . , . s = 2000, Ga = 5, Θ / Θ t = (2 . , / Θ t are shown in (a) and (c) and those with larger values of Θ / Θ t in (b) and(d). For the bedload transport conditions, the restitution coefficient has been varied to mimic the minimal ( e = 0 .
9) and nearlymaximal ( e = 0 .
01) possible effect of lubrication forces. which − P zx ˙ γ/ Γ coll ii (cid:39) (cid:28) d ), especially for bedload transport, regardless of whether transportis weak or intense. In order words, the dense algebraic layer usually does not exist. One of the reasons may be thefact that drag dissipation (Γ drag ii ), which has been neglected in Refs. [39, 40], actually dominates collisional dissipation(Γ coll ii ) in bedload transport [Fig. 1(b) in Ref. [69], which is based on the same numerical model]. B. Test of interface definition against data from our direct transport simulations
Figures 3 and 4 show that the interface z = z r defined by Eq. (6) obeys Properties 1-3 of the Bagnold interfacefor most simulated conditions. In fact, the numerical data support that most transport (80 − z r [Figs. 3(a) and 4(a)], that the bed friction coefficient µ b does not change much with τ [Figs. 3(b) and 4(b)], andthat the expression − P zx ( z r ) = τ − τ t is approximately obeyed for conditions with √ s Ga (cid:39)
10 [Figs. 3(c) and 4(c)].Furthermore, µ b varies overall between about 0 . . τ [Fig. 3(b)],which is surprisingly small given the large variability of the simulated conditions. That is, µ b can be consideredan approximate universal constant for the purpose of sediment transport modeling, which is, indeed, what we didin a recent study [72]. In contrast, interfaces defined through a constant value of φ/φ b [line-connected symbols inFig. 4(b)], through a constant value of µ [line-connected symbols in Fig. 4(c)], or through other definitions proposedin the literature (not shown) in general do not fulfill the requirements of the Bagnold interface.Conditions with √ s Ga (cid:46)
10 deviate from Property 3 [Figs. 3(c) and 4(c)], the reason for which can be seen inFig. 4(d). It shows that the local fluid shear stress τ f = τ + P zx at z r is near the flow threshold τ t at low transportstages and remains constant or decreases with increasing Θ, consistent with Property 3. However, once a criticalvalue Θ ≈ . τ f ( z r ) begins to increase and enters a regime in which it becomes proportional to Θ t τ .This proportionality causes − P zx ( z r ) /τ to approach a limiting value at large transport stages that is smaller thanthe value unity required by Property 3, with larger values of the flow threshold Shields number Θ t corresponding tolarger deviations. In fact, the sediment transport regime that exhibits the largest values of the flow threshold forcohesionless particles [max(Θ t ) ≈ .
2] is viscous bedload transport, which is characterized by comparably small valuesof √ s Ga [72]. -2 √ sGa − ( τ − τ t ) / P z x ( z r ) -2 Impact number, √ sGa B e d f r i c t i o n c o e ffi c i e n t , µ b -2 √ sGa T r a n s p o r t r a t e r a t i o , Q r / Q s = . , Ga = . = . , Ga = . = . , Ga = = . , Ga = = . , Ga = = . , Ga = = . , Ga = = . , Ga = = , Ga = . = , Ga = . = , Ga = = , Ga = = , Ga = = , Ga = = , Ga = = , Ga = = , Ga = . = , Ga = . = , Ga = = , Ga = = , Ga = = , Ga = = , Ga = = , Ga = = . , Ga = (c)(b)(a) Filled symbols : e = . : e = . FIG. 3.
Test of Bagnold interface properties.
Test of (a) Property 1, (b) Property 2, and (c) Property 3 of the Bagnoldinterface with data from our direct transport simulations for various combinations of the particle-fluid-density ratio s , Galileonumber Ga, Shields number Θ, and thus impact number √ s Ga. For conditions with s ≤ .
65 (corresponding to bedloadtransport), the restitution coefficient has been varied to mimic the minimal ( e = 0 .
9) and nearly maximal ( e = 0 .
01) possibleeffect of lubrication forces. The vertical bars indicate the range of values the quantities cover with varying Θ above about 2Θ t .This lower limit is imposed to separate the random variability due to bad statistics when Θ is close to Θ t [e.g., see Fig. 4(c)]from the actual variability. Indications that the Bagnold interface properties are obeyed: (a) the sediment transport rate ratio Q r /Q is near unity, (b) the bed friction coefficient µ b is approximately constant with Θ (relatively small vertical bars), and (c)the quantity − ( τ − τ t ) /P zx ( z r ) is near unity. IV. PHYSICAL ORIGIN OF FRICTION LAW
As explained in Sec. I B 2, there have been two interpretations of the friction law (Property 2) in the literature.In Sec. IV A, we show that the first interpretation based on the rheology of dense granular flows and suspensions ingeneral is inconsistent with data from our direct transport simulations. In particular, we present strong evidence forthe absence of a liquidlike flow regime at low transport stages. In Sec. IV B, we show that the second interpretationassociated with particle rebounds at the bed surface is consistent with the simulation data for most conditions. Inparticular, we explain why this kinematic interpretation also applies to bedload transport, in which the particledynamics are dominated by long-lasting intergranular contacts rather than particle kinematics.
A. Dense rheology interpretation of friction law
Figure 5(a) shows that the particle volume fraction φ ( z r ) at the Bagnold interface, obtained from our direct transportsimulations, increases with the Shields number Θ until it approaches at large Θ a constant maximal value that dependson whether the simulated condition corresponds to bedload ( φ maxbedl (cid:39) .
45) or saltation transport ( φ maxsalt (cid:39) . φ (cid:38) .
4, particularly when considering that the values of φ ( z r ) are near 10 − for some simulated conditionsand could possibly be even lower for conditions more extreme than those simulated. However, conditions correspondingto sufficiently intense bedload transport [e.g., conditions with √ s Ga ≤ (cid:38) t ; see ellipse in Fig. 5(a)] posea notable exception as φ ( z r ) (cid:38) .
4. For these conditions, the dense rheology interpretation of the friction law may,indeed, be consistent with the simulation data. Rescaled Shields number, Θ / Θ t Q r / Q s = . , Ga = . , e = . = . , Ga = . , e = . Rescaled Shields number, Θ / Θ t µ b s = , Ga = , e = . = , Ga = , e = . , z r ≡ z | φ = . φ b Rescaled Shields number, Θ / Θ t − ( τ − τ t ) / P z x ( z r ) s = . , Ga = , e = . = . , Ga = , e = . = . , Ga = , e = . , z r ≡ z | µ = . -2 -1 Shields number, Θ τ f ( z r ) / τ t (a) (b)(d)(c) FIG. 4.
Exemplary trends of quantities associated with Bagnold interface properties. (a) Sediment transport rateratio Q r /Q , (b) bed friction coefficient µ b , and (c) − ( τ − τ t ) /P zx ( z r ) versus rescaled Shields number Θ / Θ t . (d) Rescaled surfacefluid shear stress τ f ( z r ) /τ t versus Shields number Θ. The interface z = z r is calculated by Eq. (6) if not otherwise stated inthe legends. Symbols correspond to viscous bedload transport ( s = 2 .
65, Ga = 0 . s = 2 . s = 2000, Ga = 5). For the bedload transport conditions, the restitution coefficienthas been varied to mimic the minimal ( e = 0 .
9) and nearly maximal ( e = 0 .
01) possible effect of lubrication forces.
Absence of liquidlike granular flow regime
The simulation data indicate that a liquidlike granular flow regime does not necessarily exist. For example, Fig. 5(b)shows for saltation transport with sufficiently low Θ / Θ t (brown, dashed lines) that the local friction coefficient µ canremain well below the yield stress ratio µ s (cid:39) .
277 [81] within the dense flow region ( φ (cid:38) . φ (cid:39) .
58) to gaslike( φ (cid:46) .
4) values is, regardless of the transport regime, very thin ( < d ) at low transport stages (Fig. 4 in Ref. [51],which is based on the same numerical model). In this transient zone and slightly beyond, the average particle velocity (cid:104) v x (cid:105) and thus the particle shear rate ˙ γ obey an exponential relaxation behavior (Fig. 7 in Ref. [51]), and the Bagnoldinterface ( z = z r ) is located within this relaxation zone [Fig. 2(a) in Ref. [71], which is based on the same numericalmodel]. Hence, one may interpret the Bagnold interface as the base of the gaslike transport layer.Furthermore, an exponential relaxation of ˙ γ is reminiscent of granular creeping [91, 92, 95], which is associatedwith a nonlocal rheology [91–94]. In fact, if the rheology was local, µ would solely depend on the particle volumefraction φ or, alternatively, on the dimensionless number that characterizes the rapidness of the granular shearingmotion relative to particle rearrangement processes: the viscoinertial number [81–84] K = (cid:113) ( ρ p d ˙ γ + 2 ρ f ν ˙ γ ) /P zz ≡ (cid:112) I + 2 J. (7)The viscoinertial number K reconciles inertial granular flows, characterized by the inertial number I = ˙ γd/ (cid:112) P zz /ρ p ,with viscous suspensions, characterized by the viscous number J = ρ f ν ˙ γ/P zz . However, a data collapse of µ ( φ )and µ ( K ) is found only when Θ is sufficiently far from the flow threshold Θ t (consistent with Ref. [85]), where“sufficiently” usually refers to relatively intense transport conditions, as shown in Figs. 5(b) and 5(c) for two casesthat are exemplary for turbulent bedload (turquoise lines) and saltation transport (brown lines).Put together, the fact that µ < µ s within the dense flow region, the very thin creepinglike transient zone fromquasistatic to gaslike particle volume fractions, and the absence of a local and thus liquidlike rheology are strongevidence for a granular solid-gas transition around the Bagnold interface, where the solidlike and gaslike regime areconnected by the creepinglike zone. Note that a granular solid-gas transition and the absence of a liquidlike granularflow regime at low transport stages are rather unusual in the context of granular flows and suspensions. To ourknowledge, they have previously been found only in viscous bedload transport experiments [86]. Further note that0 Rescaled Shields number, Θ / Θ t -3 -2 -1 P a r t i c l e v o l u m e f r a c t i o n , φ ( z r ) φ maxsalt ≃ . φ maxbedl ≃ . -4 -3 -2 -1 Particle volume fraction, φ -2 -1 F r i c t i o n c o e ffi c i e n t , µ Θ < Θ t Θ ≥ Θ t Θ < . Θ t Θ ≥ . Θ t µ ≃ µ s Viscoinertial number, K F r i c t i o n c o e ffi c i e n t , µ (b) s = . , Ga = = , Ga = : Increasing Θ (c)(a) Open symbols : e = . : e = . √ sGa ≤ FIG. 5.
Failure of dense rheology interpretation. (a) Particle volume fraction φ ( z r ) at the Bagnold interface versusShields number Θ. (b, c) Friction coefficient µ versus (b) particle volume fraction and (c) viscoinertial number K . Symbolsin (a) correspond to data from our direct transport simulations for various combinations of the particle-fluid-density ratio s ,Galileo number Ga, and Θ. For symbol legend, see Fig. 3. For conditions with s ≤ .
65 (corresponding to bedload transport),the restitution coefficient has been varied to mimic the minimal ( e = 0 .
9) and nearly maximal ( e = 0 .
01) possible effect oflubrication forces. The turquoise and brown lines in (b) and (c) correspond to the conditions ( s, Ga , e ) = (2 . , , .
9) and( s, Ga , e ) = (2000 , , . the absence of a liquidlike rheology at low transport stages implies that two-phase flow models of sediment transportthat are based on local rheology models [85, 131] can be applied only to sufficiently intense transport conditions. Very viscous bedload transport
For conditions corresponding to very viscous bedload transport ( √ s Ga ≤ (cid:46) t ). In fact, for Θ (cid:38) t , boththe friction coefficient µ [Figs. 3(b) and 4(b)] and particle volume fraction φ [ellipse in Fig. 5(a)] are approximatelyconstant at z r , which is consistent with a local rheology µ ( φ ) around the Bagnold interface (i.e., liquidlike flow behaviordue to φ (cid:38) . J ( z r ) for Θ (cid:38) t , which is consistent with a localrheology µ ( J ). Consistently, Figs. 6(b) and 6(c) show exemplary for the case ( s, Ga , e ) = (2 . , . , .
01) that thesimulation data of the effective friction coefficient τ /P zz collapse as a function of J for sufficiently large Θ / Θ t , whereasthis local rheology behavior is disobeyed for small Θ / Θ t . This finding and the shape of the profiles of [ τ /P zz ]( J )shown in Figs. 6(b) and 6(c) are in qualitative agreement with recent viscous bedload transport measurements (cf.Fig. 9 in Ref. [86]).We now show that the approximate constancy of J ( z r ) for sufficiently large Θ / Θ t can be inferred from the definitionof the Bagnold interface [Eq. (6)] applied to viscous conditions. First, using µ = − P zx /P zz and the fact that the localviscous fluid shear stress can be expressed as τ f = τ + P zx = ρ f ν (1 − φ ) u (cid:48) x [51, 71], where u x is the mean horizontalfluid velocity, we obtain from Eq. (6) that the following condition must be obeyed at the Bagnold interface ( z = z r ):( P zz ˙ γ ) (cid:48) = P zz ˙ γ (cid:48) − µ (cid:48) P zx ˙ γ − ρ f νµ ˙ γ [(1 − φ ) u (cid:48) x ] (cid:48) = 0 . (8)Second, we neglect spatial changes of the particle volume fraction φ because it is close to the packing fraction in densesystems, and thus we also neglect spatial changes of µ as they are of the same order [81]. Using these approximationsand the shear rate definition ˙ γ = (cid:104) v x (cid:105) (cid:48) in Eq. (8), we approximately obtain J ( z r ) ≈ [ (cid:104) v x (cid:105) (cid:48)(cid:48) /u (cid:48)(cid:48) x ]( z r ) µ b [1 − φ ( z r )] . (9)1 Rescaled Shields number, Θ / Θ rt -2 -1 V i s c o u s nu m b e r a t i n t e r f a c e , J ( z r ) s = . , Ga = . ( √ sGa ≤ ) s = . , Ga = . ( √ sGa ≤ ) s = , Ga = . ( √ sGa ≤ ) -7 -5 -3 -1 Viscous number, J E ff e c t i v e f r i c t i o n c o e ffi c i e n t , τ / P zz Θ / Θ t = . Θ / Θ t = . Θ / Θ t = . Θ / Θ t = . Θ / Θ t = . -7 -5 -3 -1 Viscous number, J -2 -1 E ff e c t i v e f r i c t i o n c o e ffi c i e n t , τ / P zz Open symbols : e = . : e = . (a) s = . , Ga = . , e = .
01 s = . , Ga = . , e = . (b) (c) FIG. 6.
Dense rheology interpretation for very viscous bedload transport. (a) Viscous number J ( z r ) at the Bagnoldinterface versus rescaled Shields number Θ / Θ t . Symbols correspond to data from our direct transport simulations for thosecombinations of the particle-fluid-density ratio s , Galileo number Ga, and Shields number Θ that obey √ s Ga ≤
1. The twoconditions ( s, Ga , e ) = (2 . , , .
9) (turbulent bedload transport, turquoise circles) and ( s, Ga , e ) = (2000 , , .
9) (turbulentsaltation transport, brown triangles) from Figs. 5(b) and 5(c) are also shown for comparison. (b, c) Effective friction coefficient τ /P zz versus J for the case ( s, Ga , e ) = (2 . , . , .
01) and several Θ / Θ t in (b) log-linear and (c) log-log scale. The quantity [ (cid:104) v x (cid:105) (cid:48)(cid:48) /u (cid:48)(cid:48) x ]( z r ) is expected to exhibit an approximately constant value smaller than unity as the particlevelocity profile (cid:104) v x (cid:105) ( z ) is strongly coupled to the flow velocity profile u x ( z ) when the bed is fully mobile (i.e., liquidlike)due to a strong viscous drag forcing [71], which explains the approximate constancy of J ( z r ) for sufficiently large Θ / Θ t (Fig. 6a). Hence, for conditions corresponding to very viscous bedload transport ( √ s Ga ≤
1) sufficiently far from theflow threshold (Θ (cid:38) t ), µ b ≈ const can be explained in the context of dense granular flows and suspensions. B. Rebound interpretation of friction law
The gaslike transport layer is composed of particles that hop, slide, and/or roll along a solidlike granular bed at lowtransport stages or a liquidlike granular bed at large transport stages [Figs. 5(b) and 5(c)]. Except for very viscousbedload transport (which is therefore excluded from the following considerations), the hopping motion is significantand usually even dominates above the Bagnold interface ( z > z r ) [72]. Now we argue that a steady transport statein which particles hop along a granular bed (Fig. 7) causes the kinetic friction coefficient µ t ≡ − P tzx /P tzz to beapproximately constant at z r : µ tb ≡ µ t ( z r ) ≈ const. FIG. 7.
Sketch of the trajectory of a particle hopping along a granular bed.
Driven by the flow, a transported particle(blue) hops along the solidlike or liquidlike granular bed (yellow particles). Instants of particle contacts are colored deep blue,and the ones for which the center of mass of the transported particle is above the Bagnold interface ( z > z r ) are numberedconsecutively (for illustrating the mathematical derivation in the Appendix). Constant kinetic friction coefficient
First, defining the average (cid:104) A (cid:105) ↑ ( ↓ ) = φ (cid:104) AH [+( − ) v z ] (cid:105) /φ ↑ ( ↓ ) of a quantity A over ascending (descending) particles,where H the Heaviside function and φ ↑ ( ↓ ) = φ (cid:104) H [+( − ) v z ] (cid:105) the volume fraction of ascending (descending) particles,we approximately obtain φ (cid:104) v z v i (cid:105) = φ ↑ (cid:104) v z v i (cid:105) ↑ + φ ↓ (cid:104) v z v i (cid:105) ↓ ≈ φ ↑ (cid:104) v z (cid:105) ↑ (cid:104) v i (cid:105) ↑ + φ ↓ (cid:104) v z (cid:105) ↓ (cid:104) v i (cid:105) ↓ = φ ↑ (cid:104) v z (cid:105) ↑ ( (cid:104) v i (cid:105) ↑ − (cid:104) v i (cid:105) ↓ ) , (10)where we neglected velocity correlations and used the steady-state mass balance φ (cid:104) v z (cid:105) = φ ↑ (cid:104) v z (cid:105) ↑ + φ ↓ (cid:104) v z (cid:105) ↓ = 0 [69].Further using the definition of the kinetic stresses [Eq. (4b)] and (cid:104) c z c i (cid:105) = (cid:104) v z v i (cid:105) (which follows from (cid:104) v z (cid:105) = 0), wethen obtain from Eq. (10) µ t = − P tzx P tzz = − φ (cid:104) c z c x (cid:105) φ (cid:104) c z (cid:105) = − φ (cid:104) v z v x (cid:105) φ (cid:104) v z (cid:105) ≈ (cid:104) v x (cid:105) ↓ − (cid:104) v x (cid:105) ↑ (cid:104) v z (cid:105) ↑ − (cid:104) v z (cid:105) ↓ . (11)As the Bagnold interface is the effective elevation of energetic particles rebounding at the bed surface (Sec. III A),Eq. (11) implies that µ tb is a measure for the ratio between the average horizontal momentum loss [ ∝ ( (cid:104) v x (cid:105) ↓ −(cid:104) v x (cid:105) ↑ )( z r )]and vertical momentum gain [ ∝ ( (cid:104) v z (cid:105) ↑ − (cid:104) v z (cid:105) ↓ )( z r )] of hopping particles rebounding at the bed surface.Second, provided that the influence of fluid drag on the vertical motion of hopping particles can be neglected (thisprecondition is indirectly verified by the fact that the final result is consistent with data from our direct transportsimulations), a steady hopping motion requires (cid:104) v z (cid:105) ↑ ( z r ) ≈ −(cid:104) v z (cid:105) ↓ ( z r ) due to energy conservation. On average,only an approximately constant impact angle α i = − arctan[ (cid:104) v z (cid:105) ↓ / (cid:104) v x (cid:105) ↓ ]( z r ), resulting in an approximately constantrebound angle α r = arctan[ (cid:104) v z (cid:105) ↑ / (cid:104) v x (cid:105) ↑ ]( z r ), can ensure this constraint [64, 65, 97], which combined implies µ tb ≈ const. Approximate equality of friction coefficients
Until here our reasoning is largely in line with previous studies [60–65, 97]. These studies now concluded µ b ≈ constfrom µ tb ≈ const, which is consistent with our direct transport simulations, as shown in Fig. 8(a). In fact, it canbe seen that µ tb /µ b is relatively close unity for most simulated conditions, except for very viscous bedload transportconditions ( √ s Ga ≤
1) with Θ / Θ t (cid:38)
2. However, exactly for these conditions, µ b ≈ const has been explained from thelocal rheology of dense viscous suspension (Sec. IV A). Interestingly, conditions with √ s Ga ≤ / Θ t (cid:46) µ tb /µ b that are again relatively close to unity, as shown for an exemplary case in the inset of Fig. 8(a). Thissuggests that the rebound interpretation of µ b ≈ const explained in this section may actually apply to very viscousbedload transport at low transport stages even though the hopping motion is dominated by particles sliding androlling along the granular bed [72].Figure 8(b) shows that the contact friction coefficient µ c ≡ − P czx /P czz is relatively close to µ b for all simulatedconditions. Furthermore, Fig. 8(c) tests the hypothesis of previous studies [60–65, 97] that P tij ( z r ) ≈ P ij ( z r ) is thereason why µ tb ≈ µ b . It can be seen that, while this reasoning works well for saltation transport conditions, it doesnot hold for bedload transport conditions because µ tb ≈ µ b despite P tij ( z r ) (cid:28) P ij ( z r ).In the Appendix, we derive µ tb ≈ µ cb ≈ µ b from first physical principles. In summary, this derivation mainly exploitsthat the granular transport layer is gaslike, which means that collisions between particles located above the Bagnoldinterface are predominantly binary. This property allows us to write the contact stress tensor component P czi ( z r ) asthe total impulse per unit bed area per unit time generated by collisions between particles transported above theBagnold interface with bed particles below the Bagnold interface. Using Eq. (10) (which is based on the steady-statemass balance) and that the Bagnold interface ( z = z r ) is the effective elevation of energetic particle-bed rebounds(Sec. III A), it can then be shown that each such collision approximately generates the impulse equivalent per unitbed area per unit time of the associated kinetic stress tensor component P tzi ( z r ), which implies P czi ( z r ) ≈ R ↑ z r P tzi ( z r ),where R ↑ z r is the average number of such collisions per crossing of the Bagnold interface from below. As R ↑ z r is thesame for i = x and i = z , it eventually follows µ cb ≈ µ tb and thus µ tb ≈ µ cb ≈ µ b . V. DISCUSSION AND CONCLUSIONS
We have used numerical simulations that couple the discrete element method for the particle motion with a con-tinuum Reynolds-averaged description of hydrodynamics to study the physical origin and universality of theoreticalthreshold shear stress-based models of the rate of nonsuspended sediment transport for a large range of Newtonianfluids driving transport, including viscous and turbulent liquids and air. The vast majority of such models are based3 -1 Impact number, √ sGa F r i c t i o n c o e ffi c i e n t r a t i o , µ t b / µ b -1 Impact number, √ sGa F r i c t i o n c o e ffi c i e n t r a t i o , µ c b / µ b -1 Impact number, √ sGa P r e ss u r e r a t i o , [ P t zz / P zz ] ( z r ) Θ / Θ t (b) Open symbols : e = . : e = . (a) (c) FIG. 8.
Approximate equality of friction coefficients.
Friction ratios (a) µ tb /µ b and (b) µ cb /µ b and pressure ratio (c)[ P tzz /P zz ]( z r ) versus impact number √ s Ga. Symbols correspond to data from our direct transport simulations for variouscombinations of the particle-fluid-density ratio s , Galileo number Ga, and Shields number Θ. For symbol legend, see Fig. 3.For conditions with s ≤ .
65 (corresponding to bedload transport), the restitution coefficient has been varied to mimic theminimal ( e = 0 .
9) and nearly maximal ( e = 0 .
01) possible effect of lubrication forces. The vertical bars indicate the range ofvalues the quantities cover with varying Θ above about 2Θ t . This lower limit is imposed to separate the random variability dueto bad statistics when Θ is close to Θ t [e.g., see Fig. 4(c)] from the actual variability. Inset of (a): friction ratio µ tb /µ b versusrescaled Shields number Θ / Θ t for very viscous bedload transport ( s = 2 .
65, Ga = 0 . e = 0 . on, or can be reformulated in the spirit of, Bagnold’s [18–20] assumption that there is a well-defined interface be-tween granular bed and transport layer, which we have called the “Bagnold interface”, with certain special properties(Properties 1-3 in the Introduction). From our study, we have gained the following insights:1. Our simulations support the hypothesis that the Bagnold interface corresponds to the effective elevation atwhich the most energetic particles rebound, which can be mathematically defined through a maximum of thelocal production rate of cross-correlation fluctuation energy [Eq. (6)].2. Our simulations indicate that, in general, the transition between the solidlike granular bed and gaslike granulartransport layer occurs through a very thin granular creepinglike zone, which contains the Bagnold interface andwhich is associated with a nonlocal granular flow rheology. A local rheology, which is required for liquidlikebehavior, is usually found only for relatively intense transport conditions [Figs. 5(b) and 5(c)]. The absence ofa liquidlike rheology at low transport stages implies that two-phase flow models of sediment transport that arebased on local rheology models [85, 131] can be applied only to sufficiently intense transport conditions.3. As the majority of sediment transport is gaslike, the transport rate above the Bagnold interface well approximatesthe overall transport rate, as supported by our simulations [Figs. 3(a) and 4(a)] and demanded by Property 1.4. Our simulations indicate that the ratio between the particle shear stress and normal-bed pressure at the Bagnoldinterface, the bed friction coefficient µ b , varies between about 0 . . µ b is insensitive to the fluid shear stress τ , as demanded byProperty 2. The physical origin of this universal approximate invariance of µ b has been physically linked to asteady transport state in which particles continuously rebound at the bed surface (Fig. 7).5. Very viscous bedload transport ( √ s Ga (cid:46)
1) not too far above the flow threshold (Θ (cid:38) t ) poses a notableexception: our simulations indicate that the granular flow around the Bagnold interface is liquidlike [Figs. 5(a)and 6], and the friction law has been physically linked to the local rheology of dense viscous suspensions.6. As the friction law is obeyed at the base of the gaslike transport layer, µ b fundamentally differs from the constantyield stress ratio associated with the solid-liquid transition in dense granular flows and suspensions. This findingchallenges a large number of studies [18–28] according to which µ b is the yield stress ratio.7. Our simulations indicate that the local fluid shear stress τ f ( z r ) at the Bagnold interface reduces to a value nearthe flow threshold τ t at low transport stages and remains constant or decreases with increasing Shields number4Θ, consistent with Property 3. However, once a critical value Θ ≈ . τ f ( z r ) begins to increase andenters a regime in which it becomes proportional to Θ t τ . This behavior results in a deviation from Property 3for sufficiently viscous bedload transport ( √ s Ga (cid:46) τ f ( z r ) reduces to the smallest value that just allows entrain-ment of bed sediment (by the splash caused by particle-bed impacts and/or by the action of fluid forces), which isassumed to be near τ t [18–28, 54–63]. However, according to our recent study [72], τ t is not an entrainment thresholdbut rather a rebound threshold: the minimal fluid shear stress needed to compensate the average energy loss ofrebounding particles by the fluid drag acceleration during particle trajectories. That is, τ f ( z r ) reduces to the smallestvalue that just allows a long-lasting rebound motion. This interpretation (which was originally proposed by Bagnold[1] for turbulent saltation transport but later discarded) is independent of whether the bed is rigid or erodible andconsistent with our finding that µ b ≈ const is linked to a steady rebound state rather than the constant yield stressratio at the granular solid-liquid transition. In fact, based on this rebound picture, we have proposed a universalanalytical flow threshold model [72], which uses µ b = 0 .
63 (the simulation mean) and which predicts τ t for arbitraryenvironmental conditions in simultaneous agreement with available measurements in air and viscous and turbulentliquids despite not being fitted to any kind of experimental data. That is, the only ingredient that remains missingfor a universal scaling law predicting the rate of nonsuspended sediment transport [i.e., a version of Eq. (1) that isapplicable to arbitrary environmental conditions] is a universal scaling law for the average particle velocity v x in theflow direction. So far, we have succeeded in deriving an expression for v x for sufficiently low Θ / Θ t [72], and we arecurrently working on a generalization to arbitrarily large Θ / Θ t . Finally, we would like to emphasize that bed sedi-ment entrainment, even though it does not seem to affect the functional structure of the scaling laws of nonsuspendedsediment transport, is still required to sustain the equilibrium state described by such laws [72]. ACKNOWLEDGMENTS
We acknowledge support from a grant from the National Natural Science Foundation of China (No. 11750410687).
Appendix: Physical derivation of equality of friction coefficients
First, we use the steady momentum balance with respect to contact forces: − P c (cid:48) zi = ρ p φ (cid:104) a ci (cid:105) [69], where a c is theparticle acceleration due to contact forces ( F cm = (cid:80) n F mn ). Integrating this balance over elevations z > z r yields P czi ( z r ) = 1 T ∆ (cid:88) n T (cid:90) F cni H ( z n − z r )d t, (A.1)where we used (cid:82) ∞ z r δ ( z − z n )d z = H ( z n − z r ) and Eq. (2). Above the Bagnold interface ( z > z r ), the granular flowis gaslike [Fig. 5(a)], implying that particle contacts between hopping particles mainly occur during binary collisions.Because a binary contact between a particle m and a particle n does not contribute to Eq. (A.1) due to F cm + F cn = 0,the contacts contributing to Eq. (A.1) are predominantly particle-bed rebounds (colored deep blue in Fig. 7). Theterm (cid:82) T F cni H ( z n − z r )d t thus describes the total impulse gained by particle n in time T during those particle-bedrebounds in which its center of mass is located above the Bagnold interface ( z n > z r ). Consecutively numbering suchparticle-bed rebounds by r n = 1 , , ..., R nT (Fig. 7), where R nT is the total number of rebounds of particle n that occurin time T above z r , and denoting the velocity change caused by each rebound as δv r n i , which implies that ρ p V np δv r n i is the gained impulse at each rebound, we obtain from Eq. (A.1) P czi ( z r ) (cid:39) T ∆ (cid:88) n R nT (cid:88) r n =1 ρ p V np δv r n i = ρ p δv ri T ∆ (cid:88) n R nT V np , (A.2)where δv ri is the average of δv r n i over all particles and particle-bed rebounds above z r . Now we separate R nT intothe number of instants n,T ↑ z r particle n crosses the Bagnold interface from below in time T and the average number R n ↑ z r of rebounds of particle n per such crossing that occur above z r : R nT = R n ↑ z r n,T ↑ z r . Furthermore, as the Bagnoldinterface is the effective elevation of energetic particle-bed rebounds (Sec. III A), we approximate δv r n i by the averagevelocity gain at z r : δv r n i ≈ (cid:104) v z (cid:105) ↑ ( z r ) − (cid:104) v z (cid:105) ↓ ( z r ). Combining these mathematical manipulations and using Eqs. (4b)and (10), and the fact that the vertical upward-flux [ φ ↑ (cid:104) v z (cid:105) ↑ ]( z r ) of particles through the Bagnold interface equals5the total particle volume (cid:80) n n,T ↑ z r V n that crosses the Bagnold interface from below per unit bed area ∆ per unittime T , we approximately obtain from Eq. (A.2) P czi ( z r ) ≈ ρ p [ (cid:104) v z (cid:105) ↑ − (cid:104) v z (cid:105) ↓ ]( z r ) T ∆ (cid:88) n R n ↑ z r n,T ↑ z r V np = R ↑ z r ρ p [ φ ↑ (cid:104) v z (cid:105) ↑ ]( z r )[ (cid:104) v z (cid:105) ↑ − (cid:104) v z (cid:105) ↓ ]( z r ) ≈ R ↑ z r P tzi ( z r ) , (A.3)where R ↑ z r is the average number of particle-bed rebounds above z r per crossing of the Bagnold interface from below.Equation (A.3) means that the contact contribution P czi ( z r ) to the stress tensor P zi ( z r ) is approximately proportionalto the kinetic contribution P tzi ( z r ), where the proportionality factor R ↑ z r is the same for i = x and i = z . Hence,Eq. (A.3) implies µ tb ≈ µ cb ≈ µ b . [1] R. A. Bagnold, The Physics of Blown Sand and Desert Dunes (Methuen, New York, 1941).[2] M. Yalin,
Mechanics of Sediment Transport (Pergamon Press, Oxford, 1977).[3] W. H. Graf,
Hydraulics of Sediment Transport (Water Resources Publications, Littleton, CO, 1984).[4] L. C. van Rijn,
Principles of Sediment Transport in Rivers, Estuaries and Coastal Seas (Aqua, Amsterdam, 1993).[5] P. Y. Julien,
Erosion and Sedimentation (Cambridge University Press, Cambridge, 1998).[6] M. H. Garcia,
Sedimentation Engineering: Processes, Measurements, Modeling, and Practice (American Society of CivilEngineers, Reston, VA, 2007).[7] M. C. Bourke, N. Lancaster, L. K. Fenton, E. J. R. Parteli, J. R. Zimbelman, and J. Radebaugh, “Extraterrestrial dunes:An introduction to the special issue on planetary dune systems,” Geomorphology , 1–14 (2010).[8] K. Pye and H. Tsoar,
Aeolian Sand and Sand Dunes (Springer, Berlin, 2009).[9] X. Zheng,
Mechanics of Wind-Blown Sand Movements (Springer, Berlin, 2009).[10] Y. Shao,
Physics and Modelling of Wind Erosion (Kluwer, Dordrecht, 2008).[11] O. Dur´an, P. Claudin, and B. Andreotti, “On aeolian transport: Grain-scale interactions, dynamical mechanisms andscaling laws,” Aeolian Research , 243–270 (2011).[12] J. F. Kok, E. J. R. Parteli, T. I. Michaels, and D. Bou Karam, “The physics of wind-blown sand and dust,” Reports onProgress in Physics , 106901 (2012).[13] K. R. Rasmussen, A. Valance, and J. Merrison, “Laboratory studies of aeolian sediment transport processes on planetarysurfaces,” Geomorphology , 74–94 (2015).[14] A. Valance, K. R. Rasmussen, A. Ould El Moctar, and P. Dupont, “The physics of aeolian sand transport,” ComptesRendus Physique , 105–117 (2015).[15] E. Meyer-Peter and R. M¨uller, “Formulas for bedload transport,” in Proceedings of the 2nd Meeting of the InternationalAssociation for Hydraulic Structures Research (IAHR, Stockholm, 1948) pp. 39–64.[16] H. A. Einstein,
The Bed-Load Function for Sediment Transportation in Open Channel Flows (United States Departmentof Agriculture, Washington, DC, 1950).[17] M. S. Yalin, “An expression for bedload transportation,” Journal of the Hydraulics Division , 221–250 (1963).[18] R. A. Bagnold, “The flow of cohesionless grains in fluid,” Philosophical Transactions of the Royal Society London A ,235–297 (1956).[19] R. A. Bagnold, “An approach to the sediment transport problem from general physics,” in US Geological Survey Profes-sional Paper 422-I (U.S. Government Printing Office, Washington, DC, 1966).[20] R. A. Bagnold, “The nature of saltation and “bed-load” transport in water,” Proceedings of the Royal Society LondonSeries A , 473–504 (1973).[21] K. Ashida and M. Michiue, “Study on hydraulic resistance and bedload transport rate in alluvial streams,” Proceedingsof the Japan Society of Civil Engineers , 59–69 (1972).[22] F. Engelund and J. Fredsøe, “A sediment transport model for straight alluvial channels,” Nordic Hydrology , 293–306(1976).[23] A. Kovacs and G. Parker, “A new vectorial bedload formulation and its application to the time evolution of straight riverchannels,” Journal of Fluid Mechanics , 153–183 (1994).[24] Y. Nino and M. Garcia, “Gravel saltation 2. Modeling,” Water Resources Research , 1915–1924 (1994).[25] Y. Nino and M. Garcia, “Using Lagrangian particle saltation observations for bedload sediment transport modelling,”Hydrological Processes , 1197–1218 (1998).[26] G. Seminara, L. Solari, and G. Parker, “Bed load at low Shields stress on arbitrarily sloping beds: Failure of the Bagnoldhypothesis,” Water Resources Research , 1249 (2002).[27] G. Parker, G. Seminara, and L. Solari, “Bed load at low shields stress on arbitrarily sloping beds: Alternative entrainmentformulation,” Water Resources Research , 1183 (2003).[28] A. D. Abrahams and P. Gao, “A bed-load transport model for rough turbulent open-channel flows on plain beds,” EarthSurface Processes and Landforms , 910–928 (2006).[29] R. Fernandez Luque and R. van Beek, “Erosion and transport of bedload sediment,” Journal of Hydraulic Research ,127–144 (1976).[30] G. M. Smart, “Sediment transport formula for steep channels,” Journal of Hydraulic Engineering , 267–276 (1984). [31] E. Lajeunesse, L. Malverti, and F. Charru, “Bed load transport in turbulent flow at the grain scale: Experiments andmodeling,” Journal of Geophysical Research , F04001 (2010).[32] H. Capart and L. Fraccarollo, “Transport layer structure in intense bedload,” Geophysical Research Letters , L20402(2011).[33] J. J. J. Doorschot and M. Lehning, “Equilibrium saltation: Mass fluxes, aerodynamic entrainment, and dependence ongrain properties,” Boundary-Layer Meteorology , 111–130 (2002).[34] D. M. Hanes and A. J. Bowen, “A granular-fluid model for steady intense bed-load transport,” Journal of GeophysicalResearch , 9149–9158 (1985).[35] Y. Nino, M. Garcia, and L. Ayala, “Gravel saltation 1. Experiments,” Water Resources Research , 1907–1914 (1994).[36] Y. Nino and M. Garcia, “Experiments on saltation of sand in water,” Journal of Hydraulic Engineering , 1014–1025(1998).[37] F. Charru and H. Mouilleron-Arnould, “Viscous resuspension,” Journal of Fluid Mechanics , 303–323 (2002).[38] F. Charru, “Selection of the ripple length on a granular bed sheared by a liquid flow,” Physics of Fluids , 121508 (2006).[39] D. Berzi, “Analytical solution of collisional sheet flows,” Journal of Hydraulic Engineering , 1200–1207 (2011).[40] D. Berzi, “Transport formula for collisional sheet flows with turbulent suspension,” Journal of Hydraulic Engineering , 359–363 (2013).[41] F. Charru, J. Bouteloup, T. Bonometti, and L. Lacaze, “Sediment transport and bedforms: a numerical study oftwo-phase viscous shear flow,” Meccanica , 3055–3065 (2016).[42] R. Maurin, J. Chauchat, and P. Frey, “Revisiting slope influence in turbulent bedload transport: consequences for verticalflow structure and transport rate scaling,” Journal of Fluid Mechanics , 135–156 (2018).[43] J. E. Ungar and P. K. Haff, “Steady state saltation in air,” Sedimentology , 289–299 (1987).[44] M. P. Almeida, J. S. Andrade, and H. J. Herrmann, “Aeolian transport layer,” Physical Review Letters , 018001(2006).[45] M. P. Almeida, J. S. Andrade, and H. J. Herrmann, “Aeolian transport of sand,” The European Physical Journal E ,195–200 (2007).[46] M. P. Almeida, E. J. R. Parteli, J. S. Andrade, and H. J. Herrmann, “Giant saltation on mars,” Proceedings of theNational Academy of Science , 6222–6226 (2008).[47] A. Recking, P. Frey, A. Paquier, P. Belleudy, and J. Y. Champagne, “Bed-load transport flume experiments on steepslopes,” Journal of Hydraulic Engineering , 1302–1310 (2008).[48] M. Creyssels, P. Dupont, A. Ould El Moctar, A. Valance, I. Cantat, J. T. Jenkins, J. M. Pasini, and K. R. Rasmussen,“Saltating particles in a turbulent boundary layer: experiment and theory,” Journal of Fluid Mechanics , 47–74 (2009).[49] T. D. Ho, A. Valance, P. Dupont, and A. Ould El Moctar, “Scaling laws in aeolian sand transport,” Physical ReviewLetters , 094501 (2011).[50] R. L. Martin and J. F. Kok, “Wind-invariant saltation heights imply linear scaling of aeolian saltation flux with shearstress,” Science Advances , e1602569 (2017).[51] O. Dur´an, B. Andreotti, and P. Claudin, “Numerical simulation of turbulent sediment transport, from bed load tosaltation,” Physics of Fluids , 103306 (2012).[52] P. Aussillous, J. Chauchat, M. Pailha, M. M´edale, and ´E. Guazzelli, “Investigation of the mobile granular layer in bedloadtransport by laminar shearing flows,” Journal of Fluid Mechanics , 594–615 (2013).[53] S. Z. Ali and S. Dey, “Origin of the scaling laws of sediment transport,” Proceedings of the Royal Society London SeriesA , 20160785 (2017).[54] R. Kawamura, “Study of sand movement by wind,” in Translated (1965) as University of California Hydraulics EngineeringLaboratory Report , 5 (University of California, Berkeley, 1951).[55] P. R. Owen, “Saltation of uniform grains in air,” Journal of Fluid Mechanics , 225–242 (1964).[56] R. J. Kind, “A critical examination of the requirements for model simulation of wind-induced erosion/deposition phe-nomena such as snow drifting,” Atmospheric Environment , 219–227 (1976).[57] K. Lettau and H. H. Lettau, “Exploring the world’s driest climate,” in IES Report , Vol. 101 (University of Wisconsin,Madison, 1978) pp. 110–147.[58] M. Sørensen, “An analytic model of wind-blown sand transport,” Acta Mechanica Supplementum , 67–81 (1991).[59] M. Sørensen, “On the rate of aeolian sand transport,” Geomorphology , 53–62 (2004).[60] G. Sauermann, K. Kroy, and H. J. Herrmann, “A continuum saltation model for sand dunes,” Physical Review E ,031305 (2001).[61] O. Dur´an and H. J. Herrmann, “Modelling of saturated sand flux,” Journal of Statistical Mechanics , P07011 (2006).[62] T. P¨ahtz, J. F. Kok, and H. J. Herrmann, “The apparent roughness of a sand surface blown by wind from an analyticalmodel of saltation,” New Journal of Physics , 043035 (2012).[63] M. L¨ammel, D. Rings, and K. Kroy, “A two-species continuum model for aeolian sand transport,” New Journal of Physics , 093037 (2012).[64] J. T. Jenkins and A. Valance, “Periodic trajectories in aeolian sand transport,” Physics of Fluids , 073301 (2014).[65] D. Berzi, J. T. Jenkins, and A. Valance, “Periodic saltation over hydrodynamically rough beds: aeolian to aquatic,”Journal of Fluid Mechanics , 190–209 (2016).[66] H. J. Huang, T. L. Bo, and X. J. Zheng, “Numerical modeling of wind-blown sand on Mars,” The European PhysicalJournal E , 80 (2014).[67] P. Wang and X. Zheng, “Unsteady saltation on mars,” Icarus , 161–166 (2015). [68] B. Walter, S. Horender, C. Voegeli, and M. Lehning, “Experimental assessment of Owen’s second hypothesis on surfaceshear stress induced by a fluid during sediment saltation,” Geophysical Research Letters , 6298–6305 (2014).[69] T. P¨ahtz, O. Dur´an, T.-D. Ho, A. Valance, and J. F. Kok, “The fluctuation energy balance in non-suspended fluid-mediated particle transport,” Physics of Fluids , 013303 (2015), see supplementary material for theoretical derivations.[70] R. Maurin, J. Chauchat, B. Chareyre, and P. Frey, “A minimal coupled fluid-discrete element model for bedload trans-port,” Physics of Fluids , 113302 (2015).[71] T. P¨ahtz and O. Dur´an, “Fluid forces or impacts: What governs the entrainment of soil particles in sediment transportmediated by a Newtonian fluid?” Physical Review Fluids , 074303 (2017).[72] T. P¨ahtz and O. Dur´an, “The cessation threshold of nonsuspended sediment transport across aeolian and fluvial environ-ments,” Journal of Geophysical Research: Earth Surface , 1638–1666 (2018).[73] S. Courrech du Pont, P. Gondret, B. Perrin, and M. Rabaud, “Granular avalanches in fluids,” Physical Review Letters , 044301 (2003).[74] GDR MiDi, “On dense granular flows,” The European Physical Journal E , 341–365 (2004).[75] C. Cassar, M. Nicolas, and O. Pouliquen, “Submarine granular flows down inclined planes,” Physics of Fluids , 103301(2005).[76] P. Jop, Y. Forterre, and O. Pouliquen, “A constitutive law for dense granular flows,” Nature , 727–730 (2006).[77] Y. Forterre and O. Pouliquen, “Flows of dense granular media,” Annual Review of Fluid Mechanics , 1–24 (2008).[78] B. Andreotti, Y. Forterre, and O. Pouliquen, Granular Media: Between Fluid and Solid (Cambridge University Press,Cambridge, 2013).[79] P. Jop, “Rheological properties of dense granular flows,” Comptes Rendus Physique , 62–72 (2015).[80] F. Boyer, ´E. Guazzelli, and O. Pouliquen, “Unifying suspension and granular rheology,” Physical Review Letters ,188301 (2011).[81] M. Trulsson, B. Andreotti, and Philippe Claudin, “Transition from the viscous to inertial regime in dense suspensions,”Physical Review Letters , 118305 (2012).[82] C. Ness and J. Sun, “Flow regime transitions in dense non-Brownian suspensions: Rheology, microstructural characteri-zation, and constitutive modeling,” Physical Review E , 012201 (2015).[83] C. Ness and J. Sun, “Shear thickening regimes of dense non-Brownian suspensions,” Soft Matter , 914–924 (2016).[84] L. Amarsid, J.-Y. Delenne, P. Mutabaruka, Y. Monerie, F. Perales, and F. Radjai, “Viscoinertial regime of immersedgranular flows,” Physical Review E , 012901 (2017).[85] R. Maurin, J. Chauchat, and P. Frey, “Dense granular flow rheology in turbulent bedload transport,” Journal of FluidMechanics , 490–512 (2016).[86] M. Houssais, C. P. Ortiz, D. J. Durian, and D. J. Jerolmack, “Rheology of sediment transported by a laminar flow,”Physical Review E , 062609 (2016).[87] M. Houssais and D. J. Jerolmack, “Toward a unifying constitutive relation for sediment transport across environments,”Geomorphology , 251–264 (2017).[88] R. Delannay, A. Valance, A. Mangeney, O. Roche, and P. Richard, “Granular and particle-laden flows: from laboratoryexperiments to field observations,” Journal of Physics D: Applied Physics , 053001 (2017).[89] S. Roy, S. Luding, and T. Weinhart, “A general(ized) local rheology for wet granular materials,” New Journal of Physics , 043014 (2017).[90] K. Kamrin and G. Koval, “Nonlocal constitutive relation for steady granular flow,” Physical Review Letters , 178301(2012).[91] M. Bouzid, M. Trulsson, P. Claudin, E. Cl´ement, and B. Andreotti, “Nonlocal rheology of granular flows across yieldconditions,” Physical Review Letters , 238301 (2013).[92] M. Bouzid, A. Izzet, M. Trulsson, E. Cl´ement, P. Claudin, and B. Andreotti, “Non-local rheology in dense granular flows– Revisiting the concept of fluidity,” The European Physics Journal E , 125 (2015).[93] K. Nichol, A. Zanin, R. Bastien, E. Wandersman, and M. van Hecke, “Flow-induced agitations create a granular fluid,”Physical Review Letters , 078302 (2010).[94] K. A. Reddy, Y. Forterre, and O. Pouliquen, “Evidence of mechanically activated processes in slow granular flows,”Physical Review Letters , 108301 (2011).[95] M. Houssais, C. P. Ortiz, D. J. Durian, and D. J. Jerolmack, “Onset of sediment transport is a continuous transitiondriven by fluid shear and granular creep,” Nature Communications , 6527 (2015).[96] B. Allen and A. Kudrolli, “Granular bed consolidation, creep, and armoring under subcritical fluid flow,” Physical ReviewFluids , 074305 (2018).[97] D. Berzi, A. Valance, and J. T. Jenkins, “The threshold for continuing saltation on Earth and other solar system bodies,”Journal of Geophysical Research: Earth Surface , 1374–1388 (2017).[98] M. W. Schmeeckle, “Numerical simulation of turbulence and sediment transport of medium sand,” Journal of GeophysicalResearch: Earth Surface , 1240–1262 (2014).[99] J. R. D. Francis, “Experiments on the motion of solitary grains along the bed of a water-stream,” Philosophical Trans-actions of the Royal Society London A , 443–471 (1973).[100] J. E. Abbott and J. R. D. Francis, “Saltation and suspension trajectories of solid grains in a water stream,” PhilosophicalTransactions of the Royal Society of London A , 225–254 (1977).[101] D. M. Hanes and D. L. Inman, “Experimental evaluation of a dynamic yield criterion for granular fluid flows,” Journalof Geophysical Research , 3670–3674 (1985). [102] M. V. Carneiro, T. P¨ahtz, and H. J. Herrmann, “Jump at the onset of saltation,” Physical Review Letters , 098001(2011).[103] M. V. Carneiro, N. A. M. Ara´ujo, T. P¨ahtz, and H. J. Herrmann, “Midair collisions enhance saltation,” Physical ReviewLetters , 058001 (2013).[104] C. Ji, A. Munjiza, E. Avital, J. Ma, and J. J. R. Williams, “Direct numerical simulation of sediment entrainment inturbulent channel flow,” Physics of Fluids , 056601 (2013).[105] O. Dur´an, B. Andreotti, and P. Claudin, “Turbulent and viscous sediment transport - a numerical study,” Advances inGeosciences , 73–80 (2014).[106] O. Dur´an, P. Claudin, and B. Andreotti, “Direct numerical simulations of aeolian sand ripples,” Proceedings of theNational Academy of Science , 15665–15668 (2014).[107] A. G. Kidanemariam and M. Uhlmann, “Direct numerical simulation of pattern formation in subaqueous sediment,”Journal of Fluid Mechanics , R2 (2014).[108] A. G. Kidanemariam and M. Uhlmann, “Interface-resolved direct numerical simulation of the erosion of a sediment bedsheared by laminar channel flow,” International Journal of Multiphase Flow , 174–188 (2014).[109] A. G. Kidanemariam and M. Uhlmann, “Formation of sediment patterns in channel flow: minimal unstable systems andtheir temporal evolution,” Journal of Fluid Mechanics , 716–743 (2017).[110] B. Vowinckel, T. Kempe, and J. Fr¨ohlich, “Fluid-particle interaction in turbulent open channel flow with fully-resolvedmobile beds,” Advances in Water Resources , 32–44 (2014).[111] B. Vowinckel, R. Jain, T. Kempe, and J. Fr¨ohlich, “Entrainment of single particles in a turbulent open-channel flow: anumerical study,” Journal of Hydraulic Research , 158–171 (2016).[112] S. K. Arolla and O. Desjardins, “Transport modeling of sedimenting particles in a turbulent pipe flow using Euler-Lagrangelarge eddy simulation,” International Journal of Multiphase Flow , 1–11 (2015).[113] T. P¨ahtz, A. Omeradˇzi´c, M. V. Carneiro, N. A. M. Ara´ujo, and H. J. Herrmann, “Discrete element method simulationsof the saturation of aeolian sand transport,” Geophysical Research Letters , 1153–1170 (2015).[114] M. V. Carneiro, K. R. Rasmussen, and H. J. Herrmann, “Bursts in discontinuous aeolian saltation,” Scientific Reports , 11109 (2015).[115] A. H. Clark, M. D. Shattuck, N. T. Ouellette, and C. S. O’Hern, “Onset and cessation of motion in hydrodynamicallysheared granular beds,” Physical Review E , 042202 (2015).[116] A. H. Clark, M. D. Shattuck, N. T. Ouellette, and C. S. O’Hern, “Role of grain dynamics in determining the onset ofsediment transport,” Physical Review Fluids , 034305 (2017).[117] J. J. Derksen, “Simulations of granular bed erosion due to a mildly turbulent shear flow,” Journal of Hydraulic Research , 622–632 (2015).[118] J. R. Finn and M. Li, “Regimes of sediment-turbulence interaction and guidelines for simulating the multiphase bottomboundary layer,” International Journal of Multiphase Flow , 278–283 (2016).[119] J. R. Finn, M. Li, and S. V. Apte, “Particle based modelling and simulation of natural sand dynamics in the wave bottomboundary layer,” Journal of Fluid Mechanics , 340–385 (2016).[120] R. Sun and H. Xiao, “SediFoam: A general-purpose, open-source CFDDEM solver for particle-laden flow with emphasison sediment transport,” Computers & Geosciences , 207–219 (2016).[121] H. A. Elghannay and D. K. Tafti, “Les-dem simulations of sediment transport,” International Journal of Sediment Research(2017), 10.1016/j.ijsrc.2017.09.006.[122] H. A. Elghannay and D. K. Tafti, “Sensitivity of numerical parameters on DEM predictions of sediment transport,”Particulate Science and Technology , 438–446 (2017).[123] C. Gonz´alez, D. H. Richter, D. Bolster, S. Bateman, J. Calantoni, and C. Escauriaza, “Characterization of bedloadintermittency near the threshold of motion using a Lagrangian sediment transport model,” Environmental Fluid Mechanics , 111–137 (2017).[124] Z. Cheng, J. Chauchat, T. J. Hsu, and J. Calantonic, “Eddy interaction model for turbulent suspension in Reynolds-averaged Euler-Lagrange simulations of steady sheet flow,” Advances in Water Resources , 435–451 (2018).[125] P. Seil, S. Pirker, and T. Lichtenegger, “Onset of sediment transport in mono- and bidisperse beds under turbulent shearflow,” Computational Particle Mechanics , 203–212 (2018).[126] P. Gondret, M. Lance, and L. Petit, “Bouncing motion of spherical particles in fluids,” Physics of Fluids , 643 (2002).[127] F. L. Yang and M. L. Hunt, “Dynamics of particle-particle collisions in a viscous liquid,” Physics of Fluids , 121506(2006).[128] J. A. Simeonov and J. Calantoni, “Dense granular flow rheology in turbulent bedload transport,” International Journalof Multiphase Flow , 38–53 (2012).[129] N. Brodu, R. Delannay, A. Valance, and P. Richard, “New patterns in high-speed granular flows,” Journal of FluidMechanics , 218–228 (2015).[130] T. G. Drake and J. Calantoni, “Discrete particle model for sheet flow sediment transport in the nearshore,” Journal ofGeophysical Research , 19859–19868 (2001).[131] F. Chiodi, P. Claudin, and B. Andreotti, “A two-phase flow model of sediment transport: transition from bedload tosuspended load,” Journal of Fluid Mechanics755