Unified Model of Sediment Transport Threshold and Rate Across Weak and Intense Subaqueous Bedload, Windblown Sand, and Windblown Snow

Nonsuspended sediment transport (NST) refers to the sediment transport regime in which the flow turbulence is unable to support the weight of transported grains. It occurs in fluvial environments (i.e., driven by a stream of liquid) and in aeolian environments (i.e., wind‐blown) and plays a key role in shaping sedimentary landscapes of planetary bodies. NST is a highly fluctuating physical process because of turbulence, surface inhomogeneities, and variations of grain size and shape and packing geometry. Furthermore, the energy of transported grains varies strongly due to variations of their flow exposure duration since their entrainment from the bed. In spite of such variability, we here propose a deterministic model that represents the entire grain motion, including grains that roll and/or slide along the bed, by a periodic saltation motion with rebound laws that describe an average rebound of a grain after colliding with the bed. The model simultaneously captures laboratory and field measurements and discrete element method (DEM)‐based numerical simulations of the threshold and rate of equilibrium NST within a factor of about 2, unifying weak and intense transport conditions in oil, water, and air (oil only for threshold). The model parameters have not been adjusted to these measurements but determined from independent data sets. Recent DEM‐based numerical simulations (Comola, Gaume, et al., 2019; https://doi.org/10.1029/2019GL082195) suggest that equilibrium aeolian NST on Earth is insensitive to the strength of cohesive bonds between bed grains. Consistently, the model captures cohesive windblown sand and windblown snow conditions despite not explicitly accounting for cohesion.

atmospheric wind-driven (aeolian) transport of sand-sized minerals (windblown sand) and snow and ice (windblown snow).
Recently,  unified Q across most aeolian and fluvial environmental conditions (nonshallow flows), parametrized by the gravity constant g, bed slope angle α (slopes aligned with flow direction, for downward slopes, α > 0), kinematic fluid viscosity ν f , and ρ f and τ. These authors first defined the following dimensionless numbers: / , where     (1 1 / ) cos g s g is the vertical buoyancy-reduced value of g. The Shields number Θ is a measure for the ratio between tangential flow forces and normal resisting forces acting on bed surface grains , Ga controls the scaling of the dimensionless settling velocity   / s v sgd (Camenen, 2007), while s separates the settling velocity scale  sgd from the escape velocity scale  gd that grains need to exceed in order to escape the potential traps of the bed surface.  then separated Q into the transport load M (i.e., the total mass of transported grains per unit bed area) and the average streamwise sediment velocity v x via Q Mv x  and derived a parametrization for From the threshold model predictions, we also derive simple scaling laws for Θ t valid for different NST regimes. Special attention is paid to the predicted effect of α on Θ t for these regimes, which we compare with the classical bed slope correction of Θ t .
A further aspect that is addressed in this study is the effect of soil cohesiveness on Θ t and Q * . From DEMbased numerical simulations of equilibrium aeolian NST for Earth's atmospheric conditions, Comola, Gaume, et al. (2019) suggested that the strength of cohesive bonds between bed grains does neither significantly affect Θ t nor Q * even though it strongly affects the transient toward the equilibrium. However, this suggestion is very controversial. If it was indeed generally true, it would indicate a conceptual problem in most, if not all, existing aeolian transport threshold models. In fact, existing threshold models, regardless of whether they model transport initiation (e.g., Burr et al., 2015Burr et al., , 2020Iversen & White, 1982;Lu et al., 2005;Shao & Lu, 2000) or transport cessation (e.g., Andreotti et al., 2021;Berzi et al., 2016;Claudin & Andreotti, 2006;Kok, 2010;Pähtz et al., 2012), usually incorporate expressions that describe the entrainment of bed surface grains by the flow and/or grain-bed impacts, and both entrainment mechanisms are strongly hindered by cohesion. Furthermore, even those few existing models that do not consider bed sediment entrainment explicitly account for cohesive forces (Berzi et al., 2017;Pähtz & Durán, 2018a), which increase the calculated transport threshold.
Here, using numerical data provided by Comola, Gaume, et al. (2019), we show that, in contrast to previous threshold models, the conceptualization behind our threshold model supports the suggested insensitivity of aeolian NST to cohesion and propose a criterion for when cohesion can be expected to become important. Consistently, we validate the cohesionless coupled model with transport threshold and rate data not only for cohesionless aeolian and fluvial conditions but also for cohesive aeolian conditions, including aeolian NST of small mineral grains and of potentially very cohesive snow grains.
The reminder of the paper is organized as follows. Section 2 presents the transport threshold model, Section 3 the results, such as the evaluation of the coupled model with existing experimental and numerical data of transport threshold and rate, Section 4 the discussion of the results, and Section 5 conclusions drawn from it. Furthermore, there are several appendices, containing lengthy justifications of model assumptions and mathematical derivations, and a notations section at the end of the paper.
1. The mean motion of grains driven by a fluctuating turbulent flow along an inhomogeneous bed surface is represented by the motion of spherical, monodisperse grains driven by the mean (i.e., nonfluctuating) turbulent flow along the mean bed surface (i.e., grain-bed interactions are represented by their statistical mean effect). A partial justification of this idealization is presented in Appendix B 2. Given the idealization above, Θ t is defined as the smallest Shields number for which a nontrivial steady state grain trajectory exists. This definition of Θ t is conceptually equivalent to the one of Almeida et al. (2006Almeida et al. ( , 2008. It implies that, for Θ < Θ t , all transported grains lose more kinetic energy during their rebounds with the bed than they gain during their hops and thus eventually settle 3. The threshold steady state trajectory is modeled as a periodic saltation trajectory, where grain-bed interactions are described as the average grain-bed rebound, even for NST regimes in which a significant or predominant portion of grains roll and/or roll slide along the bed. However, for this trajectory to be consistent with a sustained rolling motion, we require that grains following this trajectory exhibit a sufficient kinetic energy to be able to roll (or hop) out of the most stable pockets of the bed surface assisted by the near-surface flow. A justification of this idealization is presented in Appendix C. It is partially based on a bed friction law associated with general equilibrium NST, which is introduced in Section 2.2 4. The quantities μ b and v x t  are calculated from the periodic saltation trajectory In the following subsections, we derive step-by-step the transport threshold model. First, we present basic assumptions that characterize flow, particles, and their interactions (Section 2.1). Second, we introduce the bed friction law that Equation 2b is based on and show that it leads to an expression linking the average difference between fluid and grain velocity to the bed friction coefficient μ b (Section 2.2). This friction law, when combined with insights from previous DEM-based numerical simulations of NST, supports representing the entire grain motion in equilibrium NST by grains saltating in identical periodic trajectories with rebound boundary conditions (Appendix C). Third, we present the mathematical description of this periodic saltation motion (Section 2.3). Fourth, we present the manner in which Θ t , μ b ,  x t v , and the equilibrium dimensionless sediment transport rate Q * are obtained from the family of identical periodic trajectory solutions (Section 2.4).

Flow Velocity Profile
We consider a mean inner turbulent boundary layer flow above the bed (justified in Appendix B) and assume that this flow is undisturbed by the presence of transported grains, since the nondimensionalized mass of transported sediment per unit bed area M * becomes arbitrarily small (i.e., M * → 0) when Θ approaches Θ t from above (i.e., Θ → Θ t ) because of Equation 2b. The inner turbulent boundary layer is defined by a nearly height-invariant total fluid shear stress in the absence of transported grains (i.e., dτ/dz ≃ 0, George, 2013). This definition implies that the boundary layer thickness or flow depth is much larger than the transport layer thickness (i.e., NST driven by water flows with relatively small flow depth, like in mountain streams, is excluded). The streamwise flow velocity profile u x (z) within the inner turbulent boundary layer (the law of the wall, Smits et al., 2011) is controlled by the fluid shear velocity u * and the shear Reynolds number

Fluid-Particle Interactions
Like recent numerical studies of the physics of aeolian and fluvial sediment transport (e.g., Durán et al., 2012;Schmeeckle, 2014), we consider the fluid drag and buoyancy forces as fluid-particle interactions but neglect other interaction forces because (i) they are usually much smaller than the drag force for grains in motion, (ii) there is no consensus about how these forces behave as a function of the distance from the bed surface, and (iii) we are only looking for the predominant effect and are content with an agreement between model and experimental data within a factor of 2. Details are explained in Appendix D.

Sedimentary Grains and Sediment Bed
We consider a random close packed bed made of nearly monodisperse, cohesionless, spherical sedimentary grains. Furthermore, we assume that flows have worked on this bed for a sufficiently long time such that it has reached a state of maximum resistance (Clark et al., 2017;. In this state, a quiescent bed surface is able to resist all flows but those whose largest value of the fluctuating fluid shear stress τ fluc , associated with the most energetic turbulent eddy that can possibly form (note that the size, and thus the energy content, of eddies is limited by the system dimensions, Smits et al., 2011), exceeds a certain critical resisting shear stress   Θ Y p gd (Clark et al., 2017). In particular, for laminar (i.e., nonfluctuating) fluvial conditions (τ fluc = τ), such a bed surface is able to resist all driving flows with Shields numbers Θ < Θ Y . The so-called yield stress Θ Y is therefore a statistical quantity encoding the resistance of the bed surface as a whole, though it may be interpreted as the nondimensionalized fluctuating fluid shear stress     fluc fluc Θ / ( ) p gd that is required to initiate rolling of grains resting in the most stable pockets of the bed surface (Clark et al., 2017). For a nonsloped bed of nearly monodisperse, cohesionless, frictional spheres, is expected to exhibit a universal value . Based on measurements for laminar fluvial driving flows (Charru et al., 2004;Houssais et al., 2015;Loiseleux et al., 2005;Ouriemi et al., 2007), we use the approximate value

Bed Friction Law
The transport threshold model derivation starts with describing general equilibrium NST by a bed friction law that goes back to Bagnold (1956); Bagnold (1966Bagnold ( , 1973, and which also led to the derivation of Equation 2b. In fact, the bed friction coefficient μ b in Equation 2b is rigorously linked to the streamwise (a x ) and vertical (a z ) components of the acceleration a of transported grains due to noncontact forces via (Pähtz & Durán, 2018a;) The bed surface is defined as the elevation at which p g d〈v x 〉/dz is maximal and assigned the vertical coordinate z = 0, where p g is the vertical particle pressure, v the grain velocity, and 〈⋅〉 denotes a local grain mass-weighted ensemble average (Pähtz & Durán, 2018b). The quantity describes the production rate of the cross-correlation fluctuation energy density − ρ〈v′ z v′ x 〉, where ρ is the particle concentration and v′ ≡ v−〈v〉 the fluctuation grain velocity (Pähtz & Durán, 2018b). Since the main source of production of −ρ〈v′ z v′ x 〉 are grain-bed rebounds, the bed surface elevation defined in this manner corresponds to the effective elevation of the center of mass of energetic grains when they rebound from the bed surface (Pähtz & Durán, 2018b). In DEM-based numerical simulations of NST, this elevation is about z Δ = 0.7 d above the virtual zero level of u x near threshold conditions (Pähtz & Durán, 2018a). This value is consistent with measurements of the distance Δz s ≈ 0.25 d between the crest level of static bed grains and the virtual zero level of u x in both laminar (Hong et al., 2015) and turbulent (Dey et al., 2012) flows. DEM, discrete element method; NST, nonsuspended sediment transport. v A similar link between u v x x  and v s as in Equation 7 was previously established by Bagnold (1973), while the expression for u v x x  has been validated with data from DEM-based numerical simulations of NST for a wide range of conditions (Pähtz & Durán, 2018a). Both facts support using a linearized drag law (Equation 6). Note that the terminal grain settling velocity in quiescent flow  s v , which is distinct from v s , obeys a modified version of the lower line of Equation 7 in which μ b − S is replaced by 1 (Camenen, 2007).
Equation 7, when combined with insights from previous DEM-based numerical simulations of NST into the physical nature of μ b , supports representing the entire grain motion in equilibrium NST by grains saltating in identical periodic trajectories with rebound boundary conditions (Appendix C). , in the context of periodic saltation. Below, we therefore derive expressions for both items in terms of three quantities that are intuitively associated with periodic saltation: the lift-off or rebound velocity v ↑ , impact velocity v ↓ , and hop time T.

Mathematical Description of Periodic Saltation
In periodic saltation, a grain crosses each elevation z < H, where H is the hop height, twice during a single saltation trajectory: once for times before (superscript ↑) and once for times after (superscript ↓) the instant t(H) at which it reaches its highest elevation z = H. In particular, due to mass conservation, the upward vertical mass flux of grains ϕ ↑ , which is equal to the negative downward vertical mass flux of grains, is constant for z < H, while ϕ ↑ = 0 for z > H (Berzi et al., 2016). Hence, for z < H, the concentrations of ascending and descending grains are given by   (Appendix C) and 〈v′ z v′ i 〉 = 〈v z v i 〉 because of 〈v z 〉 = 0 (mass conservation, Pähtz et al., 2015), it also follows that one can link μ b to the rebound and impact velocities: Now, we subdivide this subsection into further subsections. First, we present the deterministic laws governing the motion of a grain above the bed driven by the mean turbulent flow (Section 2.3.1). These laws directly map v ↑ to v ↓ . Second, we present the laws describing grain-bed rebounds (Section 2.3.2), mapping v ↓ back to v ↑ . For the grain trajectories to be identical and periodic, these laws must also be deterministic, which is achieved by representing them by their statistical mean effect. Third, we model the critical Shields number Θ roll needed for a grain of a given kinetic energy, associated with a given periodic saltation trajectory, to roll (or hop) out of the most stable pockets of the bed surface assisted by the near-surface flow (Section 2.3.3). For a periodic saltation trajectory to be consistent with a sustained rolling motion, we require that Θ > Θ roll .

Grain Motion Above the Bed
To make the analytical notation compact, we nondimensionalize location, velocity, acceleration, and time, indicated by a hat, using combinations of the terminal settling velocity and reduced gravity: , one then obtains the following system of differential equations describing the average trajectory from Equations 3, 5, and 6: The solution of Equations 10a-10c, with the initial condition   (0)v v , is straightforward and given in Appendix E. For the transport threshold model, the following expressions, which can be obtained from the solution (Appendix E), are crucial (written in a form that allows easy iterative evaluation, see Section 2.4): where W denotes the principal branch of the Lambert-W function.

Grain-Bed Rebounds
Grain collisions with a static sediment bed have been extensively studied experimentally (Ammi et al., 2009;Beladjine et al., 2007), numerically (Comola, Gaume, et al., 2019;Lämmel et al., 2017;Tanabe et al., 2017), and analytically (Comola & Lehning, 2017;Lämmel et al., 2017). In typical experiments, an incident grain is shot with a relatively high impact velocity (    gd v ) onto the bed and the outcome of this impact (i.e., the grain rebound and potentially ejected bed grains) statistically analyzed. We describe this process using a phenomenological description for the average rebound (vertical) restitution coefficient e ≡ |v ↑ |/|v ↓ | (e z ≡ −v ↑z /v ↓z ) as a function of sin θ ↓ = −v ↓z /|v ↓ |, with θ ↓ the average impact angle (Beladjine et al., 2007): where A = 0.87, B = 0.72, A z = 0.3, B z = 0.15, and C = 0. Equation 13c is our modification of Equation 13b, the original expression given by Beladjine et al. (2007). This modification accounts for the analytically derived asymptotic behavior of the rebound angle in the limit of small impact angle, (Lämmel et al., 2017), and for the requirement that θ ↑ → 90° when θ ↓ → 90°. Like the original expressions, the modified expressions are consistent with experimental data by Beladjine et al. (2007) for nearly monodisperse, cohesionless, spherical grains, as shown in Figures 2a-2c.
We assume that Equations 13a and 13c are roughly universal for monodisperse, cohesionless, spherical grains, independent of |v ↓ | and bed-related parameters, such as ρ p ,  g, and d, since the experimental data by Beladjine et al. (2007) have also been reproduced by a theoretical model that predicts the rebound parameters as a function of only θ ↓ (Lämmel et al., 2017). Furthermore, for conditions in which vertical drag is neg- Figures 2d and 2e show that, for the DEM numerical simulation data by Comola, Gaume, et al. (2019), e and e z are insensitive to the cohesiveness of the bed material. In fact, for a bed consisting of spherical grains with log-normally distributed size and a grain impacting with a relatively high-impact velocity (    gd v ), these authors varied the critical force F ϕ that is required to break cohesive bonds between bed grains over several orders of magnitude and found nearly no effect on the average rebound dynamics. This finding will play a crucial role in Section 4.2, where we discuss the importance of cohesion for NST.
Lastly, we note that viscous damping of binary collisions (e.g., Gondret et al., 2002), which can be important for fluvial NST, also does not seem to significantly affect the rebound laws (Pähtz et al., 2020, Section 4.1.1.4), in contrast to the assumptions in previous trajectory-based transport threshold models (Berzi et al., 2016(Berzi et al., , 2017. In particular, DEM-based simulations of aeolian and fluvial NST indicate that, even for nearly fully damped binary collisions (normal restitution coefficient ϵ = 0.01), grain-bed rebounds (and thus e and e z ) are not much affected when compared with nondamped (ϵ = 0.9) binary collisions (Pähtz & Durán, 2018a, see particularly their Movies S1-S3). This probably means that the tangential relative velocity component, which is not much affected by ϵ, dominates during the rebound process.

Critical Shields Number Required for Rolling
In Appendix C, we provide justifications for why one can represent the entire grain motion in NST, including grains that roll and/or slide along the bed, by a pure saltation motion. However, to be consistent with a sustained rolling motion of grains, we require that saltation trajectories are limited to Shields numbers Θ that are larger than the critical value Θ roll needed for a grain of a given kinetic energy, associated with a given periodic saltation trajectory, to roll (or hop) out of the most stable pockets of the bed surface assisted by the near-surface flow. In this section, we derive an expression for Θ roll using a highly simplified approach. First, since a grain located in such a pocket just changed its direction of motion from downward to upward, we assume that it exhibits the rebound kinet- where m is the grain mass, of the given periodic saltation trajectory. Second, we assume that this grain first rolls along its downstream neighbor until E ↑ has been fully converted into potential energy   * ( ) This rolling motion increases the pocket angle from the value ψ s corresponding to the most stable bed surface pocket to the value ψ * via ( Figure 3) We then model Θ roll as the critical Shields number required to push the grain out from this new pocket angle position, assuming that the mean grain motion driven by a turbulent flow is the same as the mean grain motion driven by the mean turbulent flow (the first idealization in Section 2).
Our definition of Θ roll implies that Θ roll depends on   / ( ) E mgd via Equation 14 and vanishes for sufficiently large values of   / ( ) E mgd , for which the grain can roll (or hop) out of the pocket all by itself (i.e., for ψ * ≥ 90°). Furthermore, in the limit since the yield stress Θ Y can be interpreted as the Shields number required to initiate rolling of grains resting in the most stable pockets of the bed surface in the absence of turbulent fluctuations (Section 2.1.3). This limit is relevant for Ga 2 s ≲ 1, typical for NST driven by laminar fluvial flows (Pähtz & Durán, 2018a), where the average energy of transported grains scales with  2 Ga smgd (Charru et al., 2004) and For a nonsloped (S = 0) bed of triangular or quadratic geometry and a laminar driving flow, Agudo et al. (2017) derived nearly exact analytical expressions for the critical Shields number Θ c required to push a grain out from an arbitrary pocket angle position ψ. As triangular arrangements are the most probable ones in disordered configurations (Agudo et al., 2017), we assume that these expressions approximately apply also to natural sediment beds. For the special case ψ = 25°, these expressions predict Θ c | ψ=25° ≃ 0.13. Since the nonsloped yield stress exhibits the same value, we identify this pocket angle as the one of the most stable bed surface pocket (i.e., ψ s = 25°). Agudo et al. (2017) further noted that their analytical expressions are reasonably well approximated by the nonsloped version of the model of Wiberg and Smith (1987): Θ c ∝ cot ψ for S = 0. Assuming that the general version of the model of Wiberg and Smith (1987), Θ c ∝ cot ψ − S, works reasonably well for arbitrarily sloped beds (which have not been studied by Agudo et al., 2017) where sgn denotes the sign function (note that sgn(0) = 0). Note that Θ roll = 0 does not imply that a flow with Shields number Θ = 0 is able to sustain a rolling motion of grains because Θ roll depends on the value of v ↑ , which is associated with a given periodic saltation trajectory, and for Θ = 0, such a trajectory does not exist in the first place.

Computation of Threshold and Rate of Equilibrium NST
From solving Equations 3, 7, 9, 11, 12, 13a, and 13c, we obtain a family of identical periodic trajectory solutions  Θ( , , , ) z Ga s S v . In detail, for given values of Ga, s, S, and  z v , Θ is obtained in the following manner: Ga s S v Ga s S v , we obtain the transport threshold from minimizing Θ as a function of  z v : Note that, for conditions corresponding to the bedload regime, associated with low-energy threshold trajectories (a precise definition is provided in Section 3.2), periodic trajectory solutions ). That is, the threshold trajectory calculated from Equation 16 tends to satisfy t for such conditions. From the threshold trajectory, we obtain the threshold bed friction coefficient μ bt and dimensionless average streamwise grain velocity v x t * using Equations 9, 10a, E11a, and E11b: Lastly, from μ bt , v x t * , Θ t , and c M = 1.7, we calculate the dimensionless rate Q * of equilibrium NST via Equations 2a and 2b using μ bt as the value of μ b in Equation 2b, since DEM-based numerical simulations of NST indicate that μ b does not significantly change with Θ (Pähtz & Durán, 2018b).

Model Evaluation With Experimental and Numerical Data
This section compares the model predictions with NST data from many experimental studies and with data from DEM-based numerical simulations of NST by Pähtz and Durán (2018a);  using the numerical model of Durán et al. (2012) (a snapshot and brief description are provided in Figure 4). We start with the comparison to the numerical data in order to explore the range of validity of the model. To this end, the form drag coefficient in Equation 7 is modified to the value   0.5 d C used by Pähtz and Durán (2018a); . Furthermore, since Pähtz and Durán (2018a); Pähtz and Durán (2020) simulated a quasi-two-dimensional system, parameters characterizing the bed surface also need to be modified. We did so manually and found that the following modified values lead to good overall agreement with the numerical data:   Durán et al. (2012) showing a quarter of the longitudinal simulation domain. This model couples a DEM for the particle motion under gravity, buoyancy, and fluid drag with a continuum Reynolds-averaged description of hydrodynamics. Spherical particles (∼10 4 ) with mild polydispersity are confined in a quasi-two-dimensional domain of length ∼10 3 d, with periodic boundary conditions in the flow direction, and interact via normal repulsion (restitution coefficient ϵ = 0.9) and tangential friction (contact friction coefficient μ c = 0.5). The bottom-most particle layer is glued on a bottom wall, while the top of the simulation domain is reflective but so high that it is never reached by transported particles. The Reynolds-averaged Navier-Stokes equations are combined with a semiempirical mixing length closure that ensures a smooth hydrodynamic transition from high to low particle concentration at the bed surface and quantitatively reproduces the mean turbulent flow velocity profile in the absence of transport. DEM, discrete element method; NST, nonsuspended sediment transport.  Galileo number and density ratio: Ga ∈ [0.1, 100] and s ∈ [2. 65,2000]. Furthermore, Figure 5b shows that the modified model also captures the simulated transport rate data for conditions that satisfy within a factor of 2 (closed symbols), whereas conditions that do not satisfy this constraint are not captured (open symbols). The reason behind this restriction in the model's validity range is that s 1/2 Ga is a Stokes-like number that encodes the importance of grain inertia relative to viscous drag forcing (Pähtz & Durán, 2017;Pähtz & Durán, 2018a). In fact, the validity of Equation 2a requires that, within the dense, highly collisional part of the transport layer, the motion of a transported grain is not much affected by fluid drag between two subsequent collisions with other transported grains . Note that collisions between transported grains are not relevant for the transport threshold model, which is why its validity is not affected by s 1/2 Ga (Figure 5a). Figures 5c-5g show that the nonmodified model simultaneously captures several experimental data sets corresponding to fluvial and aeolian NST within a factor of about 2. Note that all the threshold data chosen for the model evaluation were measured using methods that are known or suspected to be consistent with the definition of Θ t by Equation 2b. In the following subsection, the experimental data sets are described in detail.

Fluvial Transport Threshold Data Sets
For ). It has been argued that the resulting threshold is close to Θ t because such and similar critical transport conditions are associated with an abrupt sharp transition in the behavior of Q * as a function of Θ . Furthermore, we obtained Θ t from extrapolating the paired flume measurements of Q * and Θ by Meyer-Peter and Müller (1948, MPM48, those in Figure 5e, s ∈ (2.6, 2.7)) to vanishing Q * using the function (consistent with Equations 2a and 2b, , where we treated c 1 , c 2 , and Θ t as fit parameters.

Aeolian Transport Threshold Data Sets
For aeolian NST, we consider data that have been obtained from extrapolating equilibrium transport conditions to vanishing transport, consistent with the definition of Θ t in Equation 2b as the value at which the equilibrium value of M * vanishes. Following previous studies (Comola, Kok, et al., 2019;Martin & Kok, 2018;, we also assume that the cessation threshold of intermittent saltation is identical to Θ t . However, we do not consider measurements of the continuous transport threshold (discussed in Section 4.1) and of the threshold of saltation initiation.
Wind tunnel measurements for windblown sand from Ho (2012, H12, s ≃ 2060) and Zhu et al. (2019Zhu et al. ( , Z19, s ≃ 2210 are shown in Figure 5c, who carried out an indirect extrapolation to vanishing transport to obtain Θ t using a proxy of Q * : the surface roughness z o , which undergoes a regime shift when transport ceases. Furthermore, visual wind tunnel measurements for windblown sand and snow by Bagnold (1937, B37, s ≃ 2210), Chepil (1945, C45, s ∈ (1,400, 1800)), and Sugiura et al. (1998, S98, s ≃ 650) are shown, who obtained the cessation threshold of intermittent saltation (Θ t ) as the smallest value of Θ for which energetic mineral or snow grains fed at the tunnel entrance are able to saltate along the bed without stopping. Direct field measurements of this threshold based on the Time Frequency Equivalence Method (Wiggs et al., 2004) by Kok (2018, MK18, s ≃ 2210) are also shown. Moreover, wind tunnel measurements of Θ t from extrapolating to vanishing transport by Clifton (2006) did not feed snow at the tunnel entrance, we have chosen only their data points for freshly fallen snow. Unlike freshly fallen snow, old snow, used for the other measurements by these authors, is very cohesive (Pomeroy & Gray, 1990), and NST of old snow therefore requires a distance to reach equilibrium that is very likely much longer than the length of the wind tunnel of Clifton et al. (2006) in the absence of snow feeding (Comola, Gaume, et al., 2019). However, we have not excluded cohesive measurements if sediment feeding occurred, such as the two windblown sand data points at Ga ≈ 5 by Bagnold (1937) and Chepil (1945), corresponding to small and thus cohesive mineral grains with d ≈ 75 μm, and the windblown snow data point by Sugiura et al. (1998), corresponding to potentially very cohesive old snow.

Fluvial Transport Rate Data Sets
Since the model does not capture Q * for conditions for which s 1/2 Ga is too small (Figure 5b), we compare it only to subaqueous bedload measurements of Q * , for which s 1/2 Ga is sufficiently large. In Figures 5d and 5e, the flume measurements of Q * by Meyer-Peter and Müller (1948, MPM48, s ∈ (2.6, 2.7)), as corrected by Wong and Parker (2006), for relatively weak driving flows and small slope numbers (S ≃ 0) and by Smart and Jaeggi (1983, SJ83, s ∈ (2.6, 2.7)) and Capart and Fraccarollo (2011, CF11, s ≃ 1.51) for relatively intense driving flows and large S are shown. For all these data sets, the applied fluid shear stress is defined as τ = ρ f gh m sin α (Appendix A), where h m is the flow depth, including the sediment-fluid mixture above the bed surface z = 0, and we corrected h m for side wall drag using the method described in Section 2.3 of Guo (2014). For the latter two data sets, the model predictions can depend significantly on α and S for a given Θ. In order to make Q * only dependent on a single rather than three independent external control parameters, we have approximated Equations 2a and 2b in Figure 5e (but not in Figure 5d). We have used and f S predominantly matter when α (and thus M * ) is large. Combined, these approximations yield which are expressions independent of α and S for the rescaled transport rate  . Note that f S ≃ f α for aeolian NST and slope-driven NST in viscous liquids, implying that this slope correction of Θ and Q * is consistent with previous results for windblown sand (Iversen & Rasmussen, 1999;Wang et al., 2021), while f S < f α for slope-driven NST in turbulent liquids. These differences reflect the different physical origins of f S and f α . The physics behind f S involves the streamwise buoyancy force acting on transported grains, which is associated with the gradient of the viscous fluid shear stress (Maurin et al., 2018), while the physics behind f α involves the average streamwise fluid momentum balance (Appendix A), which is associated with the gradient of the sum of the viscous and Reynolds-averaged fluid shear stress.  Figures 5d and 5f. Note that the experiments by Sugiura et al. (1998) were carried out using potentially very cohesive old snow, while the data set by Ralaiarisoa et al. (2020) corresponds to the first controlled measurements of Q * for intense windblown sand, which are not captured by standard expressions for Q * from the literature. Furthermore, field measurements of Q * for windblown sand by Kok (2017, MK17, s ≃ 2210) are shown in Figures 5d and 5g. These authors estimated the intermittent (i.e., nonequilibrium) transport rate Q in and the fraction f Q of active windblown sand, from which we obtained the equilibrium transport rate via Q = Q in /f Q (Comola, Kok, et al., 2019).

NST Regimes
Pähtz and Durán (2018a) provided a criterion to distinguish bedload, defined as NST in which a significant portion of transported grains is moving in enduring contacts with the bed surface (e.g., via rolling and sliding), from saltation, defined as NST in which this portion is insignificant. This criterion states that saltation PÄHTZ ET AL.

10.1029/2020JF005859
14 of 36 occurs when more than 90% of the transport layer thickness are due to the contact-free motion of grains: Here, for threshold conditions, we distinguish bedload from saltation using the threshold trajectory value of the critical Shields number required for rolling roll Θ t , which is associated with the rebound energy of the threshold trajectory. When  roll Θ 0 t , grains are able to escape the pockets of the bed surface solely due to their saltation motion, that is, without the assistance of the near-surface flow. Hence, we identify this regime as saltation and distinguish it from bedload where  roll Θ 0 t . Figure 6a shows that this criterion is consistent with the one by Pähtz and Durán (2018a) for these authors' data obtained from their DEM-based simulations of NST.
Based on the hop height calculated from the transport threshold model, ), and the thickness of the viscous sublayer of the turbulent boundary layer, δ νt = 10 d/Re dt , we distinguish between viscous (H t ≲ δ νt ) and turbulent (H t ≳ δ νt ) conditions, giving rise to totally four transport regimes, which are indicated by text in Figure 6b: viscous bedload, turbulent bedload, viscous saltation, and turbulent saltation. It can be seen that viscous bedload occurs when the Stokes-like number s 1/2 Ga falls below about 75, implying that the validity criterion for the model's transport rate predictions (Equation 19) is only disobeyed for viscous bedload conditions. Note that the transition from viscous bedload to viscous PÄHTZ ET AL.  saltation coincides with a kink in the threshold curves (Figures 5a, 5c, and 6b) and that turbulent saltation occurs when s 1/4 Ga ≳ 200 (Section 3.3). Figure 6b shows that, for viscous saltation, the transport threshold approximately scales as   1/2 2 Θ ( ) t s Ga . To demonstrate the origin of this scaling, we approximate the flow velocity profile in Equation 3 as u x ≃ u * Re d z/d, since saltation trajectories are relatively large (H ≫ Z Δ d) and fully submerged within the viscous sublayer. Using this profile in Equation 12 and approximating the dimensionless terminal settling velocity by its Stokes drag limit, v s* ≃ Ga/18 (Equation 7 for small Ga), yields where ) and Θ t = 0. Hence, vertical drag is not negligible in viscous saltation. In fact, the spikes in Figure 6c indicate a substantial positive deviation from e z = 1 for viscous saltation, whereas e z is close to unity for the other regimes. Larger values of e z mean that vertical drag suppresses the vertical motion of grains more strongly, causing a slight positive deviation from μ bt = const (spikes in Figure 6d) and a substantial negative deviation from ( ) / v v z t x t 2 1 2  (spikes in Figure 6e). In Appendix C, we use these two scaling relations and their violation for viscous saltation to justify the rebound boundary conditions. We now perform a similar analysis for turbulent saltation. We approximate the flow velocity profile in Equation 3 as the log-layer profile u x ≃ κ −1 u * ln(z/z o ), where z o is the roughness height, since saltation trajectories are relatively large (H ≫ Z Δ d) and substantially exceed the viscous sublayer and buffer layer. Using this profile in Equation 12 and neglecting vertical drag (as e z ≃ 1, see Figure 6c) yields   Figure 7 shows how the slope number S affects the transport threshold predictions.

Bed Slope Dependency of Transport Threshold
For the different NST regimes, different scaling laws are found:  Ga s S f f Ga f s (25b) the validity of which is shown in Figures 7a-7c, respectively. Equation 25a follows from the fact that the term  is substantially larger than S for the threshold trajectory, since grain velocities become comparable to the terminal settling velocity v s (i.e.,   1 zt v ) because of vertical drag in viscous saltation, implying that the effect of S is small. In contrast, for turbulent NST, vertical drag can be neglected (i.e.,   1 zt v ), which leads to Equation 25b (Appendix F). Note that, for turbulent saltation, the validity of Equation 25b requires s 1/4 Ga ≳ 200, where s 1/4 Ga is a dimensionless number that roughly controls the occurrence of the transport threshold minimum for saltation (Pähtz & Durán, 2018a) Classically, the term 1 − S/tan α r , where α r is the angle of repose (i.e., the slope angle at which the entire granular bed, as opposed to only the bed surface, fails), is used to correct Θ and/or Θ t in NST (Iversen & Rasmussen, 1994;Maurin et al., 2018), in reasonable agreement with threshold measurements (Chiew & Parker, 1994;Iversen & Rasmussen, 1994). Both the slope correction term  in Equation 25b and the slope correction term 1 − S/cot ψ s in Equation 25c resemble the functional structure of 1 − S/ tan α r . However,  o b is a purely kinematic quantity and entirely unrelated to tan α r even though its value ) is very close to typical values of tan α r . Likewise, cot ψ s ≃ 2.14 is substantially larger than typical values of tan α r and the predicted slope correction for viscous bedload therefore much milder than the classical one. Consistently, Figure 7d shows that, for the viscous bedoad experiments by Loiseleux et al. (2005), the model (solid line) reproduces the measured behavior that Θ t changes only mildly with S for |S| ≲ 0.5, while the classical slope correction (dashed line) fails to capture these data. The deviations between model and measurements for |S|≳ 0.5 are likely due to the fact that |S| approaches tan α r , weakening the resistance of the bulk of the bed. In fact, once the bulk of the bed is close to yield, this will probably affect the resistance of bed surface grains via long-range correlations (since yielding is probably a critical phenomenon, , which the model does not account for. Lastly, we emphasize that the model predictions do not take into account that large bed slopes in nature (e.g., for mountain streams) are usually accompanied by very small flow depths of the order of 1d, which cause bed mobility to decrease rather than increase with S (Prancevic et al., 2014;Prancevic & Lamb, 2015).

Transport Threshold Interpretation
In Section 2, we first idealized equilibrium NST in a manner that eliminates all kinds of grain trajectory fluctuations. We then defined the transport threshold Θ t as the smallest Shields number for which a nontrivial steady state grain trajectory exists. This definition raises two important questions: • What values of the initial lift-off velocity  o v are sufficient for a grain to approach the nontrivial steady state threshold trajectory? • How sensitive is this threshold trajectory to trajectory fluctuations?
This section gives answers to these questions and discusses consequences for the physical interpretation of Θ t arising from these answers.
We reiterate that, for most bedload conditions, the model predicts  roll Θ Θ t t (Section 2.4). This implies that the average flow is barely able to sustain a rolling motion of grains out off the most stable pockets of the bed surface, meaning that grain motion would stop for an arbitrarily small negative fluctuation of the threshold rebound velocity v ↑t , since smaller grain energies make rolling out of the most stable bed surface pockets more difficult (Section 2.3.3).

10.1029/2020JF005859
18 of 36 For saltation, the model predicts a similar behavior as for bedload. To show this, we calculate the saltation trajectory evolution from equations E1a, E1b, E1d, and E1f using Equation 3 and the periodic saltation trajectory (i.e., steady) value of the dimensionless terminal settling velocity v s* , the calculation of which was described in Section 2.4. Based on this calculation, Figure 8a shows  Figures 8c and 8d for an exemplary trajectory) from those conditions that approach no motion (subcritical, southwest of the critical lines, see dashed line in Figure 8c for an exemplary trajectory). Furthermore, Figure 8b shows that the distance and  o c v obeys critical scaling behavior for sufficiently small Θ −Θ t and vanishes in the limit Θ → Θ t . Hence, arbitrarily small negative fluctuations of v ↑t will cause saltation to cease in time.
The model prediction that both bedload and saltation threshold trajectories are unstable against arbitrarily small negative fluctuations of v ↑t implies that, for realistic natural settings associated with substantial grain trajectory fluctuations, continuous NST cannot be sustained once Θ is too close to Θ t . In other words, Θ t must be strictly smaller than the continuous transport threshold Θ cont regardless of how continuous NST is defined (an unambiguous definition is currently an open problem, Pähtz et al., 2020, Section 4.1.3.2).
This subsection is subdivided into further subsections. Section 4.1.1 models Θ cont for aeolian NST and compares the model predictions with experimental data. Section 4.1.2 discusses intermittent NST dynamics occurring for Θ < Θ cont .

Continuous Transport Threshold
Understanding the behavior of Θ cont is crucial because equilibrium transport rate expressions, such as Equations 2a and 2b, are invalid for intermittent (i.e., nonequilibrium) NST conditions (Comola, Kok, et al., 2019;. To this end, let us consider NST conditions with Θ > Θ t and suppose that the system departs more and more from the equilibrium by depositing grains on the bed surface. The more grains are deposited, the more the flow will be undisturbed by the presence of transported grains. To drive this system back to equilibrium, it is required that bed surface grains are entrained and subsequently net accelerated by the undisturbed flow. We therefore propose that Θ cont corresponds to the minimal Shields number for which the average velocity of such an entrained grain  where  Θ( , , , ) z Ga s S v denotes, like before (Section 2.4), a periodic trajectory solution. Consistent with this definition, Figure 8a shows that, for saltation, the range of supercritical initial dimensionless lift-off ve- ) substantially increases with Θ/Θ t . Note that our proposed definition of Θ cont is similar to the one by Doorschot and Lehning (2002). The main and possibly only difference is that these authors' definition referred to the average grain lifting off from the bed surface, including entrained and rebounding grains, rather than only the average entrained grain.
For aeolian NST, we can model Θ cont using Equation 26 and the fact that, on average, bed surface grains entrained by impacts of transported grains are more energetic than grains entrained directly by the flow and therefore more likely to be supercritical. In fact, the former grains are literally ejected and their average ejection velocity  e v is weakly but significantly correlated with the impact velocity v ↓ (Beladjine et al., 2007).
An empirical scaling relation that is consistent with experimental data is where c e is a proportionality constant close to unity (Beladjine et al., 2007).  (Andreotti et al., 2021) or ρ f = 1.2 kg/m 3 and varying d (Carneiro et al., 2015;Ho, 2012;Martin & Kok, 2018). The symbol color indicates the value of the density ratio s. Lines in (e) correspond to model predictions of Θ cont and the transport threshold Θ t . Θ cont for aeolian NST (Andreotti et al., 2021;Carneiro et al., 2015;Ho, 2012;Martin & Kok, 2018). Note that Andreotti et al. (2021) carried out their experiments in a wind tunnel with varying air pressure and defined threshold conditions as the transition between saltation of groups of particles (bursts) to intermittent saltation of single particles (at high pressure) or no transport (at low pressure). These authors interpreted their measurements as measurements of the transport threshold Θ t , that is, they assumed Θ cont = Θ t . However, Figure 8e shows that, for many of their measured conditions, the model predictions of Θ t are nearly an order of magnitude below those of Θ cont , and only the latter are consistent with the measurements.

Intermittent NST
For Shields numbers Θ below Θ cont , NST may remain intermittent. There are two distinct kinds of transport intermittency. The first kind occurs when turbulence-driven bed sediment entrainment events associated with energetic turbulent eddies (Cameron et al., 2020;Paterna et al., 2016;Valyrakis et al., 2011;Zheng et al., 2020) generate intermittent rolling events of entrained grains . This kind of intermittency is negligible for saltation, since the transport rate of rolling grains is much smaller than that of saltating grains. However, it is important for bedload, where it is known to occur also below Θ t . However, since the average flow cannot sustain the average motion of grains, such rolling events end very quickly below Θ t . Second, turbulent fluctuation events that temporarily push Θ fluc above Θ cont cause a different kind of intermittency, which is usually associated with saltation (Comola, Kok, et al., 2019), though the exact mechanism of this intermittency depends on the physical processes behind Θ cont . In the context of our proposed definition of Θ cont , such events will cause grains entrained by grain-bed impacts to approach a periodic saltation trajectory, thus generating a saltation chain reaction that rapidly increases the transport rate Q. Once the turbulent fluctuation event is over, provided that Θ > Θ t , saltation may maintain a large value of Q for a relatively long time , though not indefinitely as trajectory fluctuations will cause grains to eventually settle (Section 4.1). Nonetheless, this leads to a substantial hysteresis of Q for Θ t < Θ < Θ cont (Carneiro et al., 2015).
Lastly, we emphasize that, although turbulence plays a crucial role for the complex intermittent behavior of NST for Θ t < Θ < Θ cont , both Θ t and Θ cont are statistical quantities referring to the grain motion averaged over long times. That is, the influence of turbulence on Θ t and Θ cont is probably relatively weak.

Importance of Cohesion
Despite not explicitly accounting for cohesion, the coupled model captures measurements of the transport threshold Θ t for small (d ≈ 75 μm), that is, cohesive windblown sand grains by Bagnold (1937) and Chepil (1945) and measurements of Θ t and the dimensionless transport rate Q * for potentially very cohesive windblown snow by Sugiura et al. (1998), as shown in Figures 5c, 5d, and 5f. In particular, the model suggests that the increase of Θ t with decreasing grain size d for sufficiently small d, which was previously attributed to cohesion (e.g., Berzi et al., 2017;Shao & Lu, 2000), is solely due to NST entering the viscous saltation regime. In this regime, which is a stronger decrease than the one (Θ t ∼ d −2 ) predicted by standard cohesion-based models (Berzi et al., 2017;Shao & Lu, 2000).
The agreement of the model with cohesive data supports the controversial suggestion by Comola, Gaume, et al. (2019) that, for equilibrium aeolian NST, Θ t and Q * are nearly unaffected by the strength of cohesive bonds between bed grains. Further support comes from the model conceptualization. In fact, for the saltation regime, to which aeolian NST belongs, the only manner in which bed grains affect the model conceptualization is via the rebound laws, since both Q *  and Θ t (Section 2) have been conceptually introduced as bed sediment entrainment-independent physical quantities. However, the rebound laws are insensitive to the strength of cohesive bonds between bed grains for the simulation data by Comola, Gaume, et al. (2019) (Figures 2d and 2e). Note that the DEM-based numerical simulations of equilibrium aeolian NST by Comola, Gaume, et al. (2019) were limited to Earth's atmospheric conditions, for which saltating grains perform relatively large hops on average (10 … 100d, see Ho et al., 2014, Table A4). However, our model conceptualization suggests that the insensitivity of equilibrium NST to cohesion is not limited to Earth's atmospheric conditions but applies to equilibrium NST in general as long as transported grains are unable to form cohesive bonds with bed grains during their contacts with the bed. This requirement is probably satisfied in the saltation regime (dashed lines in Figure 6b) as the duration of grain-bed rebounds PÄHTZ ET AL.

10.1029/2020JF005859
21 of 36 is very short. However, it may be violated in the bedload regime (solid lines in Figure 6b) as grains tend to move in enduring contact with the bed. We therefore propose that the effects of cohesion tend to become negligible once the model predicts  roll Θ 0 t , which is the criterion with which we identify the saltation regime ( Figure 6a). However, we recognize that this proposition is controversial.

Conclusions
In this study, we have proposed and evaluated a model of the two arguably most important statistical properties of equilibrium NST: the transport threshold Shields number Θ t and dimensionless transport rate Q * . The model captures, within a factor of about 2, experimental and numerical data of Θ t across the entire range of environmental conditions and experimental and numerical data of Q * across weak and intense conditions in which a small critical value of the Stokes-like number s 1/2 Ga (Equation 19) is exceeded ( Figure 5). Such conditions include subaqueous bedload, windblown sand, and windblown snow. The conditions that are not captured by the transport rate model correspond solely to viscous bedload ( Figure 6b). Note that the agreement between model and experimental data includes the first controlled measurements of Q * for intense windblown sand by Ralaiarisoa et al. (2020), which are not captured by standard expressions for Q * from the literature.
The agreement between the model and experimental data is not the result of parameter adjustment. In fact, all model parameters have been obtained from independent data sets: c M and Z Δ from DEM-based numerical simulations of NST (Pähtz & Durán, 2018a, A, B, and C from grain-bed rebound experiments (Beladjine et al., 2007), and Θ o Y and ψ s from experiments and numerical simulations of bed yielding and grain entrainment driven by laminar flows (Agudo et al., 2017;Charru et al., 2004;Houssais et al., 2015;Loiseleux et al., 2005;Ouriemi et al., 2007).
NST is a highly fluctuating physical process (Ancey, 2020b;Durán et al., 2011) because of turbulence, surface inhomogeneities, and variations of grain size and shape and packing geometry. Furthermore, the energy of transported grains varies strongly due to variations of their flow exposure duration since their entrainment from the bed (Durán et al., 2011). However, such internal variability is completely neglected in the model, since it represents equilibrium NST by grains moving in identical periodic saltation trajectories. The high predictive capability of the model therefore suggests that crucial statistical properties of NST are relatively insensitive to its internal variability.
Although the model represents threshold conditions by a continuous grain motion, we have shown that Θ t must be strictly smaller than the continuous transport threshold Θ cont for realistic natural settings (Section 4.1). In particular, a semiempirical extension of the model for aeolian NST predicts values of Θ cont that are consistent with measurements ( Figure 8e). For Θ < Θ cont , NST can exhibit complex intermittency characteristics.
The model straightforwardly provides a criterion, which we validated with numerical data from the literature (Figure 6a), that distinguishes bedload, defined as NST in which a significant portion of transported grains is moving in enduring contacts with the bed surface (e.g., via rolling and sliding), from saltation, defined as NST in which this portion is insignificant (Section 3.2). Based on the conceptualization of the model, we have proposed that, in the saltation regime, equilibrium NST is insensitive to the cohesiveness of bed grains (Section 4.2). Consistently, the cohesionless model captures measurements for aeolian saltation of cohesive sand and snow grains. In particular, the increase of Θ t with decreasing grain size d for windblown sand with d ≲ 100 μm, previously attributed to cohesion (e.g., Berzi et al., 2017;Shao & Lu, 2000), is predicted to be solely caused by NST entering the viscous saltation regime (Figure 6b), corresponding to saltation within the viscous sublayer of the turbulent boundary layer. In this regime, the model predicts     1/2 2 3 Θ ( ) t s Ga d , which is a stronger decrease than the one (Θ t ∼ d −2 ) predicted by standard cohesion-based models (Berzi et al., 2017;Shao & Lu, 2000). Hence, the model supports the controversial notion that the strength of cohesive bonds between bed grains does neither significantly affect Θ t nor Q * , which Comola, Gaume, et al. (2019)  Classically, the transport threshold has been corrected for a nonzero slope number S via Θ t = (1 − S/tan α r ) Θ t | S=0 , where α r is the angle of repose (Iversen & Rasmussen, 1994;Maurin et al., 2018). However, the model predicts that the predominant slope correction factor for turbulent NST is actually  is very close to typical values of tan α r , its physical meaning in the model is fundamentally different. It is a purely kinematic bed friction coefficient associated with the laws that describe a grainbed rebound. Furthermore, for viscous bedload, the model predicts a much milder bed slope dependency than the classical one (Equation 25c), in agreement with measurements ( Figure 7d). Note that the model predictions do not take into account that large bed slopes in nature (e.g., for mountain streams) are usually accompanied by very small flow depths of the order of 1d, which cause bed mobility to decrease rather than increase with S (Prancevic et al., 2014;Prancevic & Lamb, 2015).
Both the insensitivity of Θ t to cohesion and the insensitivity of its slope dependency to α r are features of the fact that Θ t in the model is not associated with the entrainment of bed surface grains, neither by the flow nor grain-bed impacts. Instead, Θ t in the model can be characterized as a rebound threshold, as defined and extensively discussed in a recent review on the matter .
In the future, the model may be used to reliably predict equilibrium NST in extraterrestrial environments, such as on Venus, Titan, Mars, and Pluto. However, while the model can probably be applied to Venus and Titan conditions, since they are well within the range of environmental conditions for which we evaluated the model, the application of the model to conditions with very large density ratios s ≳ 10 4 (e.g., Mars and Pluto) requires further model validation. For this reason, DEM-based numerical simulations of NST for conditions with large s are planned in the future.