The physics of sediment transport initiation, cessation, and entrainment across aeolian and fluvial environments

Predicting the morphodynamics of sedimentary landscapes due to fluvial and aeolian flows requires answering the following questions: Is the flow strong enough to initiate sediment transport, is the flow strong enough to sustain sediment transport once initiated, and how much sediment is transported by the flow in the saturated state (i.e., what is the transport capacity)? In the geomorphological and related literature, the widespread consensus has been that the initiation, cessation, and capacity of fluvial transport, and the initiation of aeolian transport, are controlled by fluid entrainment of bed sediment caused by flow forces overcoming local resisting forces, whereas aeolian transport cessation and capacity are controlled by impact entrainment caused by the impacts of transported particles with the bed. Here the physics of sediment transport initiation, cessation, and capacity is reviewed with emphasis on recent consensus-challenging developments in sediment transport experiments, two-phase flow modeling, and the incorporation of granular physics' concepts. Highlighted are the similarities between dense granular flows and sediment transport, such as a superslow granular motion known as creeping (which occurs for arbitrarily weak driving flows) and system-spanning force networks that resist bed sediment entrainment; the roles of the magnitude and duration of turbulent fluctuation events in fluid entrainment; the traditionally overlooked role of particle-bed impacts in triggering entrainment events in fluvial transport; and the common physical underpinning of transport thresholds across aeolian and fluvial environments. This sheds a new light on the well-known Shields diagram, where measurements of fluid entrainment thresholds could actually correspond to entrainment-independent cessation thresholds.


Introduction
When an erodible sediment bed is subjected to a shearing flow of a Newtonian fluid, such as air or water, bed particles may be entrained (i.e., set into motion) by the action of flow forces and then transported by the flow, initiating a process known as sediment transport. The critical conditions that are required for the initiation of sediment transport have been studied for more than two centuries [e.g., Brahms, 1757]. Dating back to the pioneering studies for water-driven transport by Shields [1936] and for wind-driven transport by Bagnold [1936Bagnold [ , 1937Bagnold [ , 1938 (summarized in his book [Bagnold, 1941]), the initiation of sediment transport in both cases has been commonly described by threshold values of the timeaveraged shear stress τ that the flow applies onto the bed [see reviews by Durán et al., 2011;Kok et al., 2012;Merrison, 2012;Dey andAli, 2018, 2019;Yang et al., 2019, and references therein]. The idea of a threshold value of τ is natural, since a necessary condition for flowdriven entrainment (or fluid entrainment) is that flow forces and/or flow-induced torques act-ing on bed surface particles must overcome resisting forces and/or torques. Consistently, for wall-bounded flows (to which sediment transport belongs) at a given shear Reynolds number Re * ≡ u * d/ν f , the shear velocity u * ≡ τ/ρ f controls the near-surface profile of the streamwise flow velocity when averaged over the entire spectrum of turbulent fluctuations [see review by Smits et al., 2011, and references therein], where ρ f is the fluid density, ν f the kinematic fluid viscosity, and d a particle diameter characteristic for bed particles. As forces resisting entrainment of a bed particle scale with the submerged gravity force (∝ (ρ s − ρ f )gd 3 ), where ρ s is the particle density and g the gravity constant, it has been common among geomorphologists to nondimenionalize τ via Θ ≡ τ/[(ρ p − ρ f )gd] [Shields, 1936], which is known as the Shields number or Shields parameter. In the aeolian research community, the threshold parameter √ Θ [Bagnold, 1941, p. 86] is also often used. Shields [1936] and numerous researchers after him have measured transport thresholds for water-driven transport [see reviews by Miller et al., 1977;Buffington and Montgomery, 1997;Paphitis, 2001;Dey and Papanicolaou, 2008;Dey and Ali, 2019;Yang et al., 2019, and references therein]. These measurements are usually summarized in a diagram showing the threshold Shields number Θ t as a function of Re * (the Shields curve Θ t (Re * )), which is known as the Shields diagram.
However, the concept of a threshold shear stress for incipient motion (i.e., for the initiation of sediment transport by fluid entrainment) has had several consistency problems. First, for wind-driven transport, the most widely used incipient motion models [Iversen and White, 1982;Shao and Lu, 2000], when applied to Martian atmospheric conditions, predict threshold shear stresses for fine sand particles that are so large that transport should occur only during rare strong Mars storms [Sullivan and Kok, 2017]. However, this prediction is contradicted by modern observations indicating widespread and persistent sediment activity [Bridges et al., 2012a,b;Silvestro et al., 2013;Chojnacki et al., 2015], even of very coarse sand [Baker et al., 2018].
A second inconsistency, which has long been known, concerns water-driven sediment transport and is tacitly acknowledged whenever the concept of an incipient motion shear stress is applied: the sediment transport rate Q (i.e., the average particle momentum per unit bed area) seems to never truly vanish for nearly any Θ > 0 in water flume experiments because of occasional strong turbulent fluctuation events causing entrainment by bursts of much-larger-than-average flow forces. That is why measurements of Θ t have relied either on indirect extrapolation methods or on vague criteria defining the value of Q (or a proxy of Q) at which transport is critical [Buffington and Montgomery, 1997]. Such criteria had been introduced even before Shields [Gilbert, 1914;Kramer, 1935]. In particular, the experiments by Paintal [1971] suggest a power law relationship between Q, appropriately nondimensionalized, and Θ for weak flows over gravel beds: Q * ≡ Q/[ρ p d (ρ p /ρ f − 1)gd] ∝ Θ 16 (it was necessary to measure Q over tens of hours for the weakest flows), which describes a dramatic but not infinitely rapid decrease of Q * with decreasing Θ. Qualitatively similar observations were reported by Helland- Hansen et al. [1974]. Largely because of Paintal's experiments, Lavelle and Mofjeld [1987] strongly argued in favor of stochastic sediment transport models that do not contain a threshold shear stress [e.g., Einstein, 1950] in a highly cited paper with the title, "Do Critical Stresses for Incipient Motion and Erosion Really Exist?" Despite the fact that many researchers have been well aware of this inconsistency, the concept of a threshold shear stress has remained alive and never been truly questioned by the majority of scientists working on water-driven sediment transport [Dey andAli, 2018, 2019;. There are two main reasons for the trust in this concept. First, above a value of Q * that roughly coincides with typical criteria defining critical transport (Q * ≈ 0.007), the relationship between Q * and Θ turns into a much milder power law [Paintal, 1971]: Q * ∝ Θ 2.5 , suggesting a clear physical meaning of the threshold Shields number associated with this transition (Θ t ≈ 0.05). Second, descriptions of water-driven sediment transport that are based on a threshold shear stress (i.e., expressions Q * (Θ) with Q * (Θ ≤ Θ t ) = 0) have been quite successful in reproducing transport rate measurements for well-controlled conditions when using very similar values of Θ t . For example, the scaling Q * ∝ (Θ − Θ t ) 3/2 by Meyer-Peter and Müller [1948] with Θ t ≈ 0.05 is one of the most widely used expressions in hydraulic engineering for gravel transport driven by water [Wong and Parker, 2006]. However, if this value of Θ t has a real physical meaning, what is it? Does it truly describe incipient motion, which has always been the predominant interpretation [see reviews by Miller et al., 1977;Buffington and Montgomery, 1997;Paphitis, 2001;Dey and Papanicolaou, 2008;Dey andAli, 2018, 2019;Yang et al., 2019, and references therein], despite the fact that Q * (Θ ≤ Θ t ) > 0 (in Paintal's experiments, Q * > 0 even for Θ ≈ 0.007 Θ t )?
A third inconsistency in the concept of an incipient motion shear stress, which also concerns water-driven sediment transport, is also old but much less well known, perhaps because one of the key papers [Graf and Pazis, 1977] is published in French language. Graf's and Pazis' measurements show that increasing the shear stress on the bed due to the water flow from zero up to a certain value τ (a transport initiation protocol) results in smaller transport rates Q than decreasing the shear stress from a larger value down to τ (a transport cessation protocol). This clearly indicates an important role of particle inertia in sustaining water-driven sediment transport. Hence, any measurement of Θ t is affected by particle inertia because, regardless of whether an initiation or cessation protocol is used, particles are already transported when Θ approaches Θ t (see the second inconsistency discussed above). Hence, Θ t is not, or at least not only, associated with fluid entrainment and thus incipient motion. The importance of particle inertia was proposed and indirectly shown even earlier, in a largely ignored study (only eight citations indexed by Web of Science today, half a century after publication) by Ward [1969]. In this study, Ward [1969] measured smaller values of Θ t for a larger particle-fluid-density ratio s ≡ ρ p /ρ f (which is a measure for particle inertia) at the same shear Reynolds number Re * . A slight downward trend of Θ t with s even existed in the pioneering experiments by Shields [1936]. Interestingly, a particle inertia effect in waterdriven sediment transport has actually been studied. It is well known, although often not considered to be crucial in the context of transport thresholds, that the flow strength at which a transported particle can come to rest at the bed surface is weaker than the one at which it can reenter transport [e.g., Francis, 1973;Reid et al., 1985;Drake et al., 1988;Ancey et al., 2002]. In contrast, another potentially important effect of particle inertia in water-driven sediment transport has not received the same attention: the interaction between particles that are already in transport and particles of the bed surface (e.g., particle-bed impacts) may support bed particle entrainment or even be predominantly responsible for it (impact entrainment).
Particle inertia and particularly impact entrainment have been widely recognized as crucial for sustaining wind-driven sediment transport since the pioneering studies by Bagnold [1941]. Yet, in contrast to water-driven transport, there seems to be a clear-cut shear stress threshold when applying an initiation protocol in wind tunnel experiments [e.g., Bagnold, 1941]. This rather curious difference between wind-driven and water-driven transport is usually not discussed in the context of incipient motion. Why is it necessary to define critical transport rates for measuring an incipient motion shear stress threshold in waterdriven transport but not in wind-driven transport? A complete description of incipient motion should be generally applicable and not limited to a subset of possible sediment transport conditions, since there is no reason to believe that the physical mechanisms involved in the entrainment of a bed particle by a turbulent flow depend much on the nature of the flow. In fact, frameworks unifying sediment transport across driving fluids (not only in regard to transport thresholds) are scarce in general (e.g., apart from modern studies, only Bagnold [1956Bagnold [ , 1973 seems to have attempted unifying water-driven and wind-driven transport conditions).
One of the most desired aspects of a general framework of sediment transport would be its ability to reliably predict the general dependency of Q * on Θ and other dimensionless environmental parameters, such as the density ratio s. However, there is an obvious problem: since measured transport rates may depend on the experimental protocol for a given condition, as was the case in the experiments by Graf and Pazis [1977] (see third inconsistency), does the concept of a general relationship even make sense? The consensus is, yes, it does make sense when referring to transport capacity (also known as transport saturation in aeolian geomorphology), which loosely defines the maximal amount of sediment a given flow can carry without causing net sediment deposition at the bed. However, a precise definition of transport capacity is very tricky and controversial [see review by Wainwright et al., 2015, and references therein]. For example, the fact that equilibrium transport rates may depend on the experimental protocol for a given condition implies that not every equilibrium transport condition is equivalent to transport capacity and that transport capacity is in some way linked to particle inertia. In fact, that the latter may be the case was recognized by no and García [1998], who numerically modeled water-driven sediment transport as a continuous motion of particles hopping along a flat wall. In particular, these authors mentioned that the capacity relation obtained from their numerical simulations contains a threshold Shields number that may not be associated with fluid entrainment, demonstrating the necessity for a good understanding of transport capacity and its relationship to particle inertia in the context of sediment transport thresholds.
While this introduction has focused on introducing issues in our understanding of fluid entrainment, shear stress thresholds, particle inertia, transport capacity, and their mutual relationships from a historical perspective, there have been major developments in these topics in the last two decades, largely because of the emergence of novel experimental designs and modeling techniques. The purpose of this review is to draw the attention of the involved research communities to these developments that, if put together, resolve the above issues and provide a largely improved conceptual understanding of sediment entrainment and transport thresholds.
A large portion of recent developments in the field can be attributed to numerical studies modeling the particle phase using the discrete element method (DEM). In comparison to other methods modeling the particle phase (e.g., continuum models), this method has the big advantage that it approximates the laws of physics at a very basic level, namely, at the level of intergrain contacts. In fact, the force laws commonly used to model intergrain contacts are known to produce system results that match experiments extremely well [e.g., Stewart et al., 2001;Lätzel et al., 2003;Clark et al., 2016]. Additionally, granular continuum models are formulated using DEM simulations [da Cruz et al., 2005] but reproduce complex experiments on granular flows often very accurately [Jop et al., 2006]. In the context of sediment transport, the main uncertainty of DEM-based models lies therefore in the modeling of the coupling between the particle phase and the Newtonian fluid driving transport. However, many of the simulations that are described in this review show that the results are often insensitive to the details of how this coupling is treated. The authors of this review thus argue that new physics uncovered by DEM-based numerical simulations are on a relatively solid footing.
To limit the scope of this review, it focuses on studies of mildly sloped beds of relatively uniform sediments unless mentioned otherwise. Also, because of the focus on physical processes involving the bed surface, this review largely concerns nonsuspended sediment transport (i.e., the fluid turbulence is unable to support the submerged particle weight), in which transported particles remain in regular contact with the bed surface (typical for particles of sand size and larger) and which is the relevant transport mode for the morphodynamics of planetary landscapes, riverscapes, and seascapes. In contrast, in suspended transport (typical for particles of silt or dust size and smaller), transported particles can remain out of contact with the bed surface for very long times (e.g., as atmospheric dust aerosols). In typical nonsuspended wind-driven (aeolian) sediment transport, many particles move in large ballistic hops and the transport layer thickness h is therefore much larger than the particle diameter d. In the aeolian geomorphology community, such hopping particles are said to move in saltation and explicitly distinguished from particles rolling and sliding along the surface. However, this terminology is not used in this review. Instead, the term saltation transport is used for general transport regimes with h d, that is, it refers to all rather than a subset of transported particles. In typical nonsuspended liquid-driven transport (henceforth referred to as fluvial transport for simplicity although this mode is not limited to fluvial environments), h is of the order of d because the largest particle hops are small. Following the fluvial geomorphology community, transport regimes with h ∼ d are termed bedload transport.
This manuscript is organized into sections that focus on specific topics (sections 2-4) followed by a summary and outlook section (section 5) and a Notation section describing the definitions of technical terms and mathematical symbols. It is noted that readers may find it useful to read section 5 first in order to organize the contents of the manuscript, and then consult sections 2-4 for more detailed information on a particular topic. Section 2 reviews recent insights into the mechanics of beginning sediment motion and fluid entrainment gained from studying sediment transport as a dense granular flow phenomenon. For example, it has become increasingly clear that granular material can flow even when a macroscopic motion does not occur, such as for a collapsed pile of sand, because of a process known as creeping, which describes an irreversible superslow granular motion associated with sporadic microscopic rearrangements. That is, it is crucial to clearly define what kind of motion one refers to when introducing sediment transport thresholds. Likewise, forces resisting the entrainment of a bed particle do not only depend on the local arrangement of bed particles but also on granular interactions with regions within the bed that are far away from the entrainment location (i.e., sediment entrainment is a nonlocal phenomenon). This is because of collective granular structures that particles can form. Section 3 reviews insights gained from recent experimental and theoretical studies showing that the fluid shear stress applied onto the bed surface alone only poorly characterizes the critical conditions required for fluid entrainment by turbulent flows. These studies have provided more suitable criteria for sediment entrainment that take into account turbulent fluctuation events and, in particular, their durations. However, section 3 also explains that a critical fluid shear stress for incipient motion does make sense when referring to the shear stress at which the fluid entrainment probability exceeds zero (which, for turbulent fluvial bedload transport, occurs much below the Shields curve [Paintal, 1971]). For example, in wind tunnel studies (but not necessarily in the field), aeolian saltation transport is initiated at about this threshold. Finally, section 4 reviews studies on the role of particle inertia in sediment transport, a topic that has very recently undergone a dramatic change. In fact, while it is well established that impact entrainment is crucial for aeolian saltation transport [see reviews by Durán et al., 2011;Kok et al., 2012;Valance et al., 2015, and references therein], very recent experimental and theoretical studies revealed that it is also crucial for sustaining fluvial bedload transport. Likewise, a very old argument by Bagnold [1941], which was forgotten or deemed unimportant, has recently been revived. Bagnold [1941] pointed out that, for aeolian saltation transport, a predominant role of impact entrainment requires that the flow is able to sustain the motion of transported particles. This is only possible if the energy loss of transported particles rebounding with the bed is compensated by their energy gain during their trajectories via fluid drag acceleration. Models that explicitly incorporate this requirement have been able to partially unify aeolian saltation and viscous and turbulent fluvial bedload transport. When combined, the insights from the studies reviewed in sections 2-4 provide a conceptual picture free of inconsistencies, which is described in section 5. For example, the shear stress threshold compiled in the Shields diagram seems to characterize the cessation of sediment bulk motion and an appropriately defined transport capacity rather than incipient motion. Section 5 also summarizes important open problems and provides a brief outlook into related problems that have not been discussed in this review, such as the effects of particle size heterogeneity on transport thresholds and bed sediment entrainment.

Yield and Flow of Dense Granular Media in the Context of Sediment Transport
In theoretical considerations of a problem as complex as the mechanics of beginning sediment motion, simplifying assumptions must be made. This often means that the granular phase is treated extremely coarsely, as a continuum with a Coulomb-like friction coefficient [Terzaghi, 1951;Drucker and Prager, 1952], or very finely, where the pocket geom-etry of individual grains sets the bed strength [Wiberg and Smith, 1987]. However, recent advances in granular mechanics have shown that Coulomb-like behavior of granular materials is inherently nonlocal, so it must be treated on intermediate length scales. This is due to the fact that the yielding condition, defined as the minimum shear stress required to achieve permanent granular flow, is set by emergent, collective networks of grains. These networks can couple different sections of the material together over large distances. The purpose of this section is to provide an overview of recent work on yield and flow of dense granular materials in the context of sediment transport, with a particular focus on the nonlocal nature of granular yielding. To simplify the discussion, it is assumed throughout this section that the granular bed is subjected to a constant bed shear stress (like for laminar flows), in which case the existence of a fluid entrainment threshold associated with bed failure does make sense. However, this is no longer true for turbulent flows, as reviewed in section 3. For more information on dense granular flow, readers might consult recent reviews [Forterre and Pouliquen, 2008;Jop, 2015;Kamrin, 2018] devoted exclusively to the topic of dense granular flow. For the connection between granular flow and sediment transport, the perspective and review by Church [2009, 2011] are also recommended.

Yielding of Granular Media
Surface grains sit in pockets on top of the bed, and the geometry of the pocket determines the entrainment conditions for that particular grain via its protrusion (i.e., the grain height above surrounding grains) and friction angle. When the downstream drag force from the fluid overcomes resistive forces from gravity and from contact forces with the pocket, the grain will begin to move. This conceptually simple scenario appears in many theoretical studies [e.g., Wiberg and Smith, 1987;Ling, 1995;Dey, 1999;Dey and Papanicolaou, 2008;Ali and Dey, 2016]. However, this picture has several conceptual problems. For example, there are many different pocket geometries [Kirchner et al., 1990;Buffington et al., 1992] implying a distribution of entrainment thresholds. Kirchner et al. [1990] made a similar argument, advocating for a statistical treatment of pocket geometries, where only the grains with the smallest entrainment thresholds would be relevant. Additionally, when transport thresholds are discussed, one typically does not include transient behavior, after the flow has pushed grains from less stable to more stable pockets. For example, an entrained grain that then restabilizes in a nearby pocket would not constitute sediment transport. After such a rearrangement, the resulting bed would have a different intergrain force and contact structure, which would be more suited to resisting the applied flow forces [Masteller and Finnegan, 2017]. Thus, determining the fluid entrainment threshold amounts to determining the strongest bed that can be formed by the grains, subject to the flow forces and dynamics. This process necessarily involves transient behavior, as grains search for stable configurations, and spatial correlations, since information about each grain's movement is transmitted through the intergrain force network.
While this represents a very challenging problem, it is exactly the picture that has emerged in recent years regarding the physical origin of frictional behavior in noncohesive soils or sediments. The yield criterion of granular materials is defined by the maximum internal shear stress that a granular material can achieve, but grains must rearrange to find this maximum stress, sometimes for a long time [Clark et al., 2018;Srivastava et al., 2019]. The yield criterion has the form of a friction coefficient, where flow occurs only when µ ≡ τ p /P > µ s , where τ p and P are the granular shear stress and pressure P, respectively, that arise from intergrain contacts, and µ s is the static friction coefficient of the material. At first glance, this is not surprising, since the grains themselves have a surface friction coefficient µ g . However, µ s is only weakly dependent on µ g [da Cruz et al., 2005], as shown in Figure 1. Even frictionless spheres have µ s ≈ 0.1 [Peyneau and Roux, 2008a,b], which arises from a preferred orientation for intergrain contacts that aligns with the compressive direction of the applied shear deformation. This effect is independent of whether the grains interact via linear spring forces [Thompson and Clark, 2019] or more realistic Hertzian interactions [Peyneau and Roux, 2008a]. Similar behavior is observed for grains and Kamrin and Koval [2014] showing a measurement of the bulk static friction coefficient µ s as a function of µ g , which is the static friction coefficient between the surfaces of two grains (simulated as two-dimensional disks).
with surface friction and irregular shape [Radjaï et al., 1998;Radjaï, 2010, 2014;Trulsson, 2018], but the maximum stress anisotropy is enhanced by these effects, since graingrain contacts can have both normal and tangential components. This raises the yield stress slightly: frictional disks have µ s ≈ 0.2 − 0.3 [da Cruz et al., 2005] and frictional spheres have µ s ≈ 0.3 − 0.4 [Jop et al., 2006], with only a weak dependence on µ g for µ g > 0.1. Additionally, µ s is nearly independent of polydispersity [Voivret et al., 2009]. This picture assumes grains are slowly moving with persistent intergrain contacts, but µ s can be lowered significantly for more energetic kinds of driving, like vibration [Gaudel and De Richter, 2019] or in aeolian saltation transport , probably because the tendency of the contact orientation to align with the compressive direction is somewhat suppressed . Thus, frictional behavior in granular media arises primarily from the anisotropic structure of force and contact networks, and grain-grain friction, shape, and polydispersity play secondary roles.
Here, µ is used to denote the local nondimensional shear stress in the granular material itself, while the Shields number Θ is the dimensionless shear stress applied to the granular bed surface, so the two quantities are not equivalent but are closely related. At the surface of the bed, µ ≈ Θ if lift forces are neglected. The existence of a maximum shear stress that can be supported by a granular material (which is independent of grain size) suggests that, for noncohesive sediments, there should be a theoretical upper limit to the threshold Shields number Θ t , Θ max t ≈ µ s . This implies that the Shields curve must plateau at low values of the shear Reynolds number Re * for laminar flows. This fact has been a subject of debate for many years, with some authors [Shields, 1936;Mantz, 1977;Miller et al., 1977;Yalin and Karahan, 1979;Govers, 1987;Buffington and Montgomery, 1997;Dey, 1999;Hong et al., 2015] showing a trend where Θ t continues to grow as Re * gets smaller, while other studies [Wiberg and Smith, 1987;Paphitis, 2001;Pilotti and Menduni, 2001;Ouriemi et al., 2007] show a plateau at low Re * . Recent work by the present authors [Clark et al., 2015aPähtz and Durán, 2018a] has investigated sediment transport thresholds over a wide range of Re * and density ratio s using simulations based on the DEM to model noncohesive grains that are coupled to fluid-driven shear forces. These studies all suggest that Θ t is a constant at low Re * and s, corresponding to the strongest possible state of the bed. It is noted that cohesive effects become important for very small grains, which can cause Θ t to continue to grow for smaller Re * .

Open Problem: Value of Viscous Yield Stress Θ max
t Measured values of the viscous yield stress Θ max t vary substantially. For nearly monodisperse beds of spherical particles, most studies reported Θ max t ≈ 0.12 [Charru et al., 2004;Loiseleux et al., 2005;Ouriemi et al., 2007;Seizilles et al., 2014;Houssais et al., 2015], but larger values of up to about 0.37 have also been reported [Lobkovsky et al., 2008;Hong et al., 2015]. Also, some measurements suggest that Θ max t depends on the median grain size [Hong et al., 2015], in contradiction to the grain size independence of µ s , while other studies find no such dependence [Ouriemi et al., 2007]. To the authors' knowledge, there is currently no convincing explanation for these contradicting observations. However, the scatter in the reported values for Θ max t (between 0.12 and 0.37) is within the range reported for the yield stress of granular materials, ranging from low-friction spheres to rougher, more frictional particles. Thus, the yield stress of the bulk granular material may at least play some role in setting the scatter in Θ max t . In this context, it is worth noting that, for the entrainment of particles resting on an idealized substrate by a laminar flow, threshold Shields numbers range from zero to very large values depending on the packing arrangement [Agudo et al., 2017;Deskos and Diplas, 2018;Topic et al., 2019;Shih and Diplas, 2019].

Rheological Descriptions
The existence of a yield stress is one piece of a rheological description, which is a constitutive law that mathematically connects the strain rate to the local stress at each point in a material. For granular materials, dissipation implies that more force is required for faster strain rates, so µ will increase with strain rate γ. For the case of sediment transport, formulation of a constitutive law has obvious practical benefits, namely that it would allow an analytical prediction of transport rates Q at varying Shields number Θ for transport conditions dominated by granular interactions. However, note that a bulk constitutive law may not be able to capture certain cases, particularly very near to the onset or cessation of fluvial bedload or aeolian saltation transport, where the transport layer is dominated by the isolated motion of a single grain along the bed (which is the typical situation in gravel-bed rivers [Parker, 1978;Phillips and Jerolmack, 2016]). Despite the fact that the force and contact networks discussed above are spatially extended, some progress has been made by considering so-called local rheologies. Based on dimensional analysis, da Cruz et al. [2005] showed that µ for dry, uniform granular flows must depend on γ via a single dimensionless number, I ≡ γd/ P/ρ p , where I is called the inertial number, similar to the Savage [Savage, 1984] or Coulomb [Ancey et al., 1999] numbers. A functional form for µ(I) can then be measured from experiments or DEM simulations (a crude approximation is given by µ = µ s +c I I, where c I is a constant parameter). If one then assumes that a three-dimensional, tensorial generalization of this law is locally satisfied at each point in space in arbitrary geometries, then the equations of motion are closed and one can predict (at least numerically) flow in any arbitrary geometry where the forces and boundary conditions are known. Experimental measurements of rapid, dense flow in several geometries show good agreement with the local rheology [MiDi, 2004;Jop et al., 2005Jop et al., , 2006].

Open Problem: Rheology of Nonsuspended Sediment Transport
There are many physical mechanisms that are relevant to nonsuspended sediment transport that are not included in the inertial number description, but recent work has suggested that appropriate dimensional analysis can be used to find a general rheological description that is relevant in all contexts. For example, viscous effects from the fluid can be included [Boyer et al., 2011;Trulsson et al., 2012;Sun, 2015, 2016;Houssais et al., 2016;Amarsid et al., 2017;Houssais and Jerolmack, 2017;Guazzelli and Pouliquen, 2018] by replacing the inertial number I with the viscous number J ≡ ρ f ν f γ/P. This description is valid when the Stokes-like number I 2 /J is small, and the standard µ(I) rheology again takes over for large I 2 /J. This crossover can be heuristically written in terms of a viscoinertial number K ≡ J + c K I 2 , where c K is an order-unity fit parameter [Trulsson et al., 2012;Sun, 2015, 2016;Amarsid et al., 2017], and the rheology takes the form µ(K).
The previous paragraph describes a unification of dry and wet, viscous granular flows, but some situations, like turbulent bedload or aeolian saltation transport, do not fit neatly into this description. Maurin et al. [2016] showed that, for intense turbulent bedload transport, the inertial number I (used for dry flows) collapses the data best, but with a different µ(I) relation compared to dry flows. Additionally, the presence of more severe velocity fluctuations and grain-grain collisions can weaken the material, giving a µ that is smaller than would be predicted by a µ(I) or µ(K) rheology at a given shear rate [Pähtz and Durán, 2018b]. Another option is to build a rheological description that explicitly accounts for these fluctuations and collisions via the Péclet number Pe ≡ γd/ √ T , where the granular temperature T equals the mean square of kinetic particle velocity fluctuations. The advantage of Pe is that it is applicable to a wide range of different granular flows (e.g., it unifies intense fluvial bedload and aeolian saltation transport), whereas K is limited to relatively homogeneous flows. The disadvantage is that Pe involves another granular property (T) that requires modeling.

Creep and Nonlocal Rheologies
As discussed in section 1, some water flume experiments suggest that fluvial bedload transport never truly ceases for nearly any Θ > 0, which is usually attributed to turbulent fluctuations. However, as discussed in this section, the granular material itself may be partially responsible. In fact, it is well known that granular creep can be observed in a variety of observational geophysical contexts [Boulton and Hindmarsh, 1987;Pierson et al., 1987;Ferdowsi et al., 2018] as well as more idealized granular flows in a laboratory setting [Roering et al., 2001;Komatsu et al., 2001;Nichol et al., 2010;Moosavi et al., 2013;Amon et al., 2013], including sediment transport explicitly [Houssais et al., 2015;Allen and Kudrolli, 2018], as depicted in Figure 2. Generally, creeping refers to slow, typically intermittent flow (not limited to the bed surface) that occurs below a macroscopic yield criterion.
One class of creeping flow involves systems where regions with µ > µ s and µ < µ s exist nearby each other, which often occurs in systems with stress gradients (e.g., due to gravity or curvature). In this case, creeping flow is observed in regions with µ < µ s [Fenistein and van Hecke, 2003;MiDi, 2004;Crassous et al., 2008;Koval et al., 2009]. This creeping flow is not steady or continuous, but occurs in a series of intermittent, avalanche-like slips, which are triggered by the nearby steadily flowing region with µ > µ s . The time-averaged shear rate profiles decay quasi-exponentially with spatial distance to the steadily flowing region. Various nonlocal theories have been proposed [Baran et al., 2006;Pouliquen and Forterre, 2009] that include a spatial length scale ξ over which flow can be triggered in this way. The most successful theories [Kamrin and Koval, 2012;Henann and Kamrin, 2013;Kamrin and Henann, 2015;Bouzid et al., 2013Bouzid et al., , 2015 suggest that the cooperative length scale ξ diverges at the yield stress (i.e., ξ ∝ | µ − µ s | −ν , where ν ≈ 0.5). This means that, near the yield stress, flow events can be triggered over arbitrarily large distances; this point is revisited below. The grain-scale physical origin of the nonlocal models and associated spatial correlations [Zhang and Kamrin, 2017] as well as how exactly to best mathematically formulate a nonlocal rheology [Bouzid et al., 2017;Li and Henann, 2019] is still a subject of debate in the literature.
The creeping flow captured by these nonlocal models is also apparent in laboratory flumes used to model fluvial sediment transport. Houssais et al. [2015Houssais et al. [ , 2016 showed that sediment transport involves the coexistence of three regimes: a dilute suspension above the bed surface, the bedload layer at the bed surface, and creeping behavior below the surface. These regions are depicted in Figure 2a. The shear rate profile in the creeping regime follows an exponential decay, which is consistent with the predictions of nonlocal models. Similar behavior was also observed by Allen and Kudrolli [2017], shown in Figure 2b, who also stressed that the apparent agreement with nonlocal models formulated for dry granular ma- terials implies that the fluid stress is not playing a major role in the observed creeping behavior. In the creeping regime, µ < µ s , but flow events are triggered via the bedload transport regime at the top of the bed via spatial correlations in the force network. These creeping events, although slow and intermittent, can lead to segregation effects over long times (∼10-100 hr), where large particles are sorted to the top [Ferdowsi et al., 2017]. Thus, creep and nonlocal rheology may play a crucial role in armoring of gravel-bedded rivers, as opposed to size sorting in the transported layer. Additionally, recent computational work [Pähtz and Durán, 2018b] has shown that sediment transport rheology is nonlocal even relatively far from the sediment transport threshold.
There is a second class of creeping flow, which is currently not explained by any rheological model. In the above discussion, creeping granular flow at µ < µ s was always induced by nearby regions with µ > µ s . In some cases, creeping flow can be observed at µ < µ s without any apparent granular flow nearby at µ > µ s [Amon et al., 2013]. This class of creep is often accompanied by compaction of the bed. Slow shear and compaction interact in a complex way that is not fully understood but can be crucial in regulating slow (e.g., millimeters to meters per day) geophysical flows [Moore and Iverson, 2002]. Similar behavior was also observed in laboratory sediment transport experiments by Houssais et al. [2015] and further studied by Allen and Kudrolli [2018], as shown in Figure 3. The latter authors observed a granular bed with an overlying laminar shear flow and showed that slow (less than 0.1 grain diameters in 90 min) creeping flow persisted even for Θ Θ t (meaning that µ < µ s everywhere in the granular bed). The grain motion in the direction of fluid flow followed an exponential decay with depth, similar to the creep described by nonlocal models. However, it was not induced by granular flow but somehow by the laminar fluid flow. Streamwise creep was also accompanied by compaction of the bed, which can strengthen the material and thus reduce creep. This second class of creep is therefore similar to com- where the color goes from dark to light with increasing particle movement. We see that particles move even at no shear, but there is greater movement at higher shear. Looking more closely in a short segment where we can track all the particles t = 30-90 seconds are the displacement of particles as a function of distance from the surface in (c) the flow x and (d) gravity z directions. We see an exponential behavior in the flow movement while the bed shifts down linearly with depth compacting uniformly. The inset of (d) is the measure of the strain γ z from fitting the slopes of (d).  Fig. 1(d) in order to have well-defined shear history conditions. Here grains which remain within y = ±0.16d are tracked and analyzed. The magnitude of displacement of the individual grains s in the plane over this time interval is denoted using the colormap to capture the bed evolution. One observes that both examples rearrange, including the one with no applied shear, with greater motion occurring for τ * /τ * c = 0.8. We examine the displacements inside the bed by observing the motion of the grains over a time t = 30-90 s. The grain displacements in the same flow and gravity directions are plotted as a function of depth z in Figs. 4(c) and 4(d), respectively. One observes that the bed creeps forward faster and settles further with increasing shear stress. Moreover, the creep along the flow direction appears to decay exponentially with depth as shown by the fits in Fig. 4(c). The decay length from the exponential fit in the case of the higher shear rates, where a meaningful variation occurs, is found to be 2.5d ± 0.1d. This decay is similar to the length scale over which grain speeds exponentially decay into the bed for τ * > τ * c [10] and was observed to be common to dry granular beds in gravity which are sheared horizontally at the top [27,28].

A. Rearrangements with depth
At the same time, the linear compaction with depth at all shear rates implies that the bed settles uniformly as grains rearrange in gravity. Such a linear increase would imply that the volume fraction of the bed increases uniformly into the bed, an issue we will examine more closely later in the discussion. The strain gradient γ z = − z/z obtained from the linear fit is shown in the inset of 074305-7 paction [Knight et al., 1995;Ribière et al., 2005] and creep [Divoux et al., 2008;Candelier and Dauchot, 2009] that is induced by tapping or vibrations, despite the fact that no explicit vibrations were applied. The existence of this class of creep implies that sediment is likely always transported (albeit slowly) for arbitrarily small values of Θ, even in the absence of turbulence. Another recent experimental flume study [Masteller and Finnegan, 2017] showed a similar result, where conditioning a bed by applying weak fluid flow led to zero net transport but a smoother bed profile with fewer protruding grains. Then, when the fluid flow rate was increased to a value associated with significant transport for a conditioned bed, sediment transport rates were smaller when compared with an unconditioned bed.

Open Problem: Physical Origin of Creeping Below Macroscopic Yield
The physical mechanisms that lead to the second class of creep, where µ < µ s everywhere in the system, are not known. One possible mechanism is contact aging [Jia et al., 2011], where the microscopic contact structure between two solid objects (i.e., grains) can evolve and weaken with time for reasons that are not fully understood [Liu and Szlufarska, 2012]. Additionally, Pons et al. [2016] showed that this second class of creep could be induced in dry granular flow by applying small pressure fluctuations to the interstitial air, with resulting shear rates of the order of 10 −7 . Similar fluctuations likely always exist in natural systems. These two hypotheses are supported by the fact that, to the authors' knowledge, this class of creep does not occur in DEM simulations, which use a Cundall-Strack model [Cundall and Strack, 1979] or similar Coulomb-like yield criterion for the frictional forces between grains, and fluctuating forces or slow variations in grain-grain friction are not included. Some DEM studies have observed creeping below a macroscopic yield criterion like the angle of response [Ferdowsi et al., 2018], but the results from these studies seem to always include some region of µ > µ s .

Critical Behavior and Weak Links
Many experimental and computational studies [Carneiro et al., 2011;Heyman et al., 2013;Houssais et al., 2015] have observed that, near sediment transport thresholds (including the impact entrainment threshold, reviewed in section 4.1.3), the time t conv required for some system measurement (e.g., the sediment transport rate Q) to converge to its steady state value appears to grow very large. A common form [Clark et al., 2015a] to capture these long time scales is t conv ∝ |Θ − Θ t | −β , where β is some positive exponent. A diverging time scale can arise in many ways, but one possibility is a critical phase transition. The study of phase transitions, where a material abruptly changes as a control parameter is smoothly varied, originated in thermal physics (e.g., liquid-gas or ferromagnetic transition), but it has also been successful in describing many other kinds of systems where thermal physics is not applicable. The key feature of a critical phase transition is a diverging correlation length, such that small changes near the critical point can have system-spanning effects that last for arbitrarily long times. The system is thus said to be scale-free at the critical point, since there is no largest length or time scale that is affected by a perturbation.

Open Question: Is Flow-Induced Bed Failure a Critical Phenomenon?
Bed failure at the yield stress describes by definition a phase transition, but whether this transition is critical and how it arises from grain-grain and grain-fluid interactions remain open questions. However, there is a growing body of work [Clark et al., 2018;Srivastava et al., 2019;Thompson and Clark, 2019] suggesting that the yielding transition for granular media is a critical transition. This is also suggested by the diverging correlation length ξ ∝ | µ − µ s | −ν that is present in the nonlocal models discussed above [Kamrin and Koval, 2012;Bouzid et al., 2013]. In addition to describing creeping flow for µ < µ s , nonlocal theories are also able to correctly predict other size-dependent effects, like strengthening of thin layers [MiDi, 2004;Kamrin and Henann, 2015]. The idea that yielding of granular media is a critical transition helps to explain certain experimentally observed behaviors in laboratory and computational models of sediment transport. For example, using a laboratory flume near the viscous limit, Houssais et al. [2015Houssais et al. [ , 2016 found a diverging time scale near the critical Shields number that is "associated with the slowing down, and increasing variability, of the particle dynamics; it is unrelated to hydrodynamics." Evidence of scale-free channeling patterns [Aussillous et al., 2016] was also observed during erosion of granular beds, which was attributed to the fact that the onset of erosion was behaving like a critical phase transition.
When the physics controlling the onset of grain motion is no longer just the yield strength of the granular material itself, then the picture changes somewhat. For example, once particle inertia becomes important in sustaining nonsuspended sediment transport (see section 4), the granular phase may not have a frictional state µ that is close to µ s , and thus it may be far from the critical point. For viscous bedload transport (small Re * ), when particle inertia is not important, computational studies typically show that t conv obeys system size dependence that is consistent with a critical phase transition Clark et al., 2018]. However, under steady driving conditions, when grain inertia starts to play a role (e.g., for larger Re * ), then t conv still diverges, t conv ∝ |Θ − Θ t | −β , but systems of different sizes will have the same t conv [Clark et al., 2015a. Thus, Θ t for inertial particles appears to be more similar to a dynamical instability rather than a true critical point.
However, nonlocal effects still likely play a role in the initiation of permanent bed failure. For example, if particle inertia plays a crucial role in sustaining sediment transport, as argued below in section 4, then a bed could be above the threshold needed to sustain motion but not have any way to get started. Returning to the argument from Kirchner et al. [1990] discussed above, if only the grains with the lowest entrainment thresholds are susceptible to being moved by the fluid, then these grains might be thought of as weak links in the bed. Motion that is initiated by these weak links could trigger flow elsewhere in the system, via the redistribution of forces or by collision. Clark et al. [2015aClark et al. [ , 2017 showed that the initiation of motion did indeed obey statistics consistent with a Weibullian weakest link scenario.

Summary
This section has described recent advances in the physics of sheared granular flows, with a focus on application to sediment transport. The main ideas are as follows. First, the yield condition for granular materials (e.g., a sediment bed) has the form of a static friction coefficient µ s , but it is not set directly by grain-grain friction. Instead, µ s is an emergent property that arises from the maximum structural anisotropy that the grain-grain contact network can support. Friction plays a minor role in determining this maximum anisotropy, and grain shape and polydispersity also play minor roles. Second, although these contact networks are extended in space (and thus inherently nonlocal), local rheological descriptions (i.e., constitutive laws) can be very successful in many contexts. Recent advances suggest that a unified, local rheological description might be within reach. This rule could be used to model any context of wet or dry granular flow with appropriate boundary conditions. Such a description could be used to predict sediment transport rates and thresholds if the grain properties (i.e., size distribution, friction coefficient, grain shape, etc.) were known, even approximately. Third, the inherently nonlocal nature of yielding is dominant when the material is near its yield condition. This causes creeping behavior in regions where a local rheology would predict no flow, which complicates the search for a unified rheological description. However, the results described in Figures 2 and 3 showed that creeping is similar in wet and dry flows, since it very slow and thus dominated by grain rearrangements (not fluid). This suggests that the nonlocal descriptions for wet and dry flows might also be unified in a relatively simple way. The underlying physics behind this nonlocal behavior is not fully understood, but there is mounting evidence that yielding of granular materials represents a kind of critical transition, where different parts of the system can be correlated over arbitrary distances. Remarkably, for sediment transport, creep seems to occur even much below the yield transition, that is, for seemingly arbitrarily small Shields numbers Θ.
This section has considered only sediment beds sheared by nonfluctuating flows and usually neglected the effects of particle inertia in sustaining sediment transport. That is, except for the occurrence of creep, many of the results of this section do not apply to turbulent flows nor flows with significant particle inertia effects that are near the threshold for grain motion (occurring for sufficiently large Re * and/or s, see section 4). In particular, the average fluid shear stress at which turbulent flows are able to entrain bed particles is usually much below the yield stress of the granular phase. Nonetheless, both creep and the viscous yields stress Θ max t will play crucial roles in the new conceptual picture of sediment transport thresholds and sediment entrainment that is presented in section 5.

Fluid Entrainment by Turbulent Flows
This section reviews the state of the art on the entrainment of bed particles by a turbulent flow of Newtonian fluid. This process is not equivalent to the initiation of overall sediment motion, which occurs even in the absence of bed sediment entrainment because of creeping (see section 2.3). It is also not equivalent to the comparably simple physics of fluid entrainment by a nonfluctuating flow. For example, when a laminar flow of a Newtonian fluid shears a target particle resting on the sediment bed, there are critical values of the fluid shear stress τ, which depend on the local bed arrangement, above which this particle begins to roll and slide, respectively [Agudo et al., 2017;Deskos and Diplas, 2018]. Once motion begins, resisting forces weaken and, since the flow does not fluctuate, the particle will inevitably leave its bed pocket (i.e., become entrained). The entrained particle will travel along the bed until it comes to rest in another pocket in which it can resist the flow, provided such a pocket exists and is accessible (when the sediment bed has yielded, particles can no longer find stable resting place, see section 2.1). In contrast, in turbulent flows, even though resist-ing forces weaken when a bed particle becomes mobilized, such a mobilized particle may not find its way out of its initial bed pocket (i.e., incomplete entrainment). The prototype for this situation is a turbulent fluctuation of the flow that exerts a large force on the particle, but the fluctuation is too short-lived for it to become entrained. Hence, there are two important ingredients that need to be considered to accurately describe sediment entrainment by turbulent flows for a given pocket geometry: the magnitude and duration of turbulent fluctuations (evidence for this statement is briefly reviewed in section 3.1). Only entrainment criteria that account for both aspects are able to accurately describe fluid entrainment experiments (section 3.2). Shear stress-based criteria, in general, do not belong to this category. Yet one can still define the critical shear stress τ In t above which the probability of fluid entrainment exceeds zero. This and related thresholds have received a lot of attention in studies on aeolian and planetary transport (section 3.3).

The Role of Turbulent Fluctuations in Fluid Entrainment
Turbulent fluctuations have been known to play a crucial role in fluid entrainment for a long time. For example, Einstein and El-Samni [1949], and later Mollinger and Nieuwstadt [1996], measured large fluctuating lift forces on a fixed rough surface induced by pressure gradient fluctuations of the order of the mean pressure gradient. These authors concluded that such pressure gradient fluctuations must be important also for the mobilization of bed sediment. In fact, numerous laboratory, field, and theoretical studies have advocated the viewpoint that the magnitude of peaks of the instantaneous flow force acting on a bed particle, consisting of both lift and drag forces, is a key aspect of fluid entrainment [e.g., Kalinske, 1947;Sutherland, 1967;Paintal, 1971;Heathershaw and Thorne, 1985;Apperley and Raudkivi, 1989;Kirchner et al., 1990;Nelson et al., 1995;Papanicolaou et al., 2001;Sumer et al., 2003;Zanke, 2003;Hofland et al., 2005;Schmeeckle et al., 2007;Vollmer and Kleinhans, 2007;Giménez-Curto and Corniero, 2009;Dwivedi et al., 2010a,b;Cameron et al., 2019Cameron et al., , 2020. However, while such force peaks explain certain observations, such as the episodic character of very weak turbulent bedload transport [Paintal, 1971;Helland-Hansen et al., 1974;Hofland, 2005] or the strong increase of weak turbulent bedload transport in the presence of vegetation [Yager and Schmeeckle, 2013;Yang and Nepf , 2018, 2019], they do not explain all observations. In fact, experiments in which a target particle was placed on an idealized rough substrate and exposed to an electrodynamic force revealed that very high force pulses do not lead to entrainment if their duration is too short [Diplas et al., 2008]. Likewise, moderate force pulses that only barely exceed resisting forces lead to entrainment if their duration is sufficiently long. That the duration of force peaks is as important as their magnitude has also been experimentally confirmed both for particles resting on idealized, fixed beds [Diplas et al., 2008;Celik et al., 2010Celik et al., , 2013Celik et al., , 2014Valyrakis et al., 2010Valyrakis et al., , 2011Valyrakis et al., , 2013Valyrakis, 2013] and natural erodible sediment beds [Salim et al., 2017[Salim et al., , 2018. However, note that, for sediment transport along erodible beds (with the exception of viscous bedload transport), the vast majority of entrainment events are triggered by particle-bed impacts, except for very weak transport conditions (see sections 4.1.2 and 4.1.3). In the following, criteria are reviewed that account for both the magnitude and duration of turbulent fluctuation events.

Impulse Criterion
The initiation of movement of a target particle resting in a pocket of the bed surface necessarily requires that the instantaneous flow forces (or torques) F(t 0 ) acting on it at the instant t 0 of initial motion overcome resisting forces (or torques) F c : However, this criterion is not sufficient for entrainment to occur as the target particle may merely move back to its initial resting place if F(t) becomes subcritical for times t too soon after t o so that its gained kinetic energy is insufficient to overcome the potential barrier of its bed pocket. For this reason, Diplas et al. [2008] proposed that the fluid impulse I f associated with larger-than-critical flow forces must exceed a critical value: where T is the duration of the impulse event (i.e., the duration of the particle acceleration phase of a turbulent fluctuation event). Note that T can be much smaller than the time needed to leave the bed pocket as the latter also includes the particle deceleration phase. Diplas et al. [2008] confirmed their hypothesis with idealized experiments in which they subjected an isolated target particle with a constant electrodynamic, horizontal force F D for a given time T D , for which I f = F D T D . In fact, their measured data of the force that is required for entrainment roughly obey the relationF D ≡ F D /F min where F min D is the minimal force required for measurable particle motion (but not necessarily entrainment) and T max D the associated time that is needed for F min D to cause entrainment ( Figure 4). In order to use equation (2) for predicting particle entrainment, one needs to know the impulse threshold I f c . For entrainment into a rolling motion, Valyrakis et al. [2010] derived an expression for the critical impulse I f c = F t T t (F t is defined below) assuming a constant pulse of a hydrodynamic force, separated into a horizontal drag and vertical lift component (F = (F D , F L )), of short duration T t (so that the angular displacement ∆ψ of the particle remains small for t ∈ (t 0 , t 0 + T t )): where F t = F D sin ψ + F L cos ψ and F n = −F D cos ψ + F L sin ψ are the tangential and normal components, respectively, of the driving flow force at the rest position, m p = 1 6 ρ p πd 3 is the particle mass, F tc = m p g cos(ψ + α)/sin ψ − (m p g/s) cot ψ the resisting force, L arm the lever arm length, C m = 1/2 the added mass coefficient, and f (ψ, α, s) = cos(ψ + α) sin α + [1 − sin(ψ + α)](cos α − 1/s), with α the bed slope angle and ψ the pivoting angle ( Figure 5). For many conditions, this expression can be well approximated by Lee et al. [2012] derived an alternate expression for short turbulent fluctuation events. Instead of a pure rolling motion, they considered entrainment into a combined rolling and sliding motion (however, note that rolling is usually the preferred mode of entrainment) without bed slope (α = 0), assuming that the associated tangential motion is described by a Coulomb friction law with friction coefficient µ C . Furthermore, instead of the pivoting angle, they described the pocket geometry by the horizontal (∆X) and vertical (∆Z) particle displacement (in units of d) that is needed for the particle to escape (equivalent to ψ +α = π/2 in Figure 5). The expression by Lee et al. [2012] reads where F e = F D (sin ψ − µ cos ψ) + F L (cos ψ + µ C sin ψ) is an effective hydrodynamic force and F ec = m p g(1 − 1/s)(sin ψ + µ C cos ψ) its critical value. For entrainment into a hopping motion, defined as a lift force-induced particle uplift by a vertical distance ≥ 1d, Valyrakis et al. [2010] derived where the resistance force is given by F Lc = m p g(1 − 1/s) cos α. Note that equation (6) with α = 0 is equivalent to equation (5) if the critical dimensionless displacement ∆Z + µ C ∆X = 1 and F L(c) replaced by F e(c) .

Equations (3)-(6) reveal that the impulse threshold I f c is constant only if the driving flow force is very strong (F(t)
F c ). However, for near-critical fluctuation events (F(t) → F c ), I f c diverges. This motivates the introduction of an energy-based entrainment criterion.

Energy Criterion
The impulse criterion (equation (2)) accounts for the available momentum of the turbulent fluctuation event in comparison to the momentum required for entrainment. However, close observation of near-bed turbulence reveals that fluctuation events are scarcely ever square pulses or even single-peaked [Valyrakis, 2013]. Instead, turbulent flows in nature exhibit a wide range of flow patterns and structures, some of which may be more efficient for particle entrainment than others. For example, the transfer of energy from flow to particles in turbulent fluctuation events with large driving flow forces (F(t) F c ) is expected to be much more efficient than in fluctuation events with near-critical flow forces (F(t) ∼ F c , see section 3.2.1). This motivates the characterization of entrainment using the energy of the fluctuation event that is effectively transferred to the particle : where W c is the minimal amount of work required for complete particle entrainment and ] the instantaneous flow power, parameterized by the cube of the local flow velocity, and C eff is the coefficient of energy transfer efficiency of the turbulent fluctuation event. The energy transfer coefficient C eff is expected to increase with F /F c (see section 3.2.1), where · denotes the time average over the event. Water flume experiments on the entrainment of a particle resting on an idealized substrate have confirmed that C eff tends to increase with F /F c ( Figure 6). However, one has to keep in mind that C eff incorporates also other effects such as grain orientation and shape.
In order to use equation (7) for predicting particle entrainment, one needs to know the energy threshold W c . Valyrakis et al. [2013] derived Hopping: For typical sediment beds, the ratio between both energy thresholds ([1 − sin(ψ + α)]L arm /d) is of the order of 0.1, demonstrating that a rolling motion is much more easily initiated upon entrainment than a hopping motion. Note that, in contrast to the expressions for the critical impulse for rolling (equations (3) and (4)), equation (8) does neither require the assumption of a small angular particle displacement ∆ψ during the acceleration phase of a turbulent fluctuation event nor the assumption of a short duration of this phase.

Shear Stress Threshold of Incipient Motion and Initiation of Aeolian Saltation Transport
The entrainment criteria reviewed in section 3.2 are able to predict whether a certain turbulent fluctuation event is capable of entraining a target particle, whereas a criterion based on a critical shear stress would not suffice for this purpose. However, one can still define a shear stress threshold τ In t (the initiation threshold) at which the fluid entrainment probability exceeds zero (i.e., below which entrainment never occurs). Such a threshold must exist because the size of turbulent flow eddies is limited by the system dimensions, such as the boundary layer thickness δ. In fact, a limited size of turbulent flow eddies implies that also Figure 6. (a-c) Flow power P f (t) versus time t for three different turbulent fluctuation events that lead to entrainment of a target particle resting on a prearranged substrate. The solid lines corresponds to experimental data . The dashed lines indicate the start of the respective fluctuation event. The dotted lines indicate the critical flow power that must be exceeded in order to overcome the resisting forces (i.e., u > u c (t)), which depend on time because resisting forces weaken once the target particle starts to move. the magnitude of peaks of the flow force is limited. That is, one can always find a nonzero shear stress below which even the largest fluctuation peaks do not exceed the resisting forces acting on bed particles (however, note that the existence of sufficiently large flow force peaks does not guarantee a nonzero entrainment probability because their durations may always be too short). Like for Θ max t , transient behavior associated with the flow temporarily pushing particles from less stable to more stable pockets is excluded in the definition of τ In t , which implies Θ In t Θ max t for laminar flows at sufficiently low shear Reynolds number Re * . Furthermore, surface inhomogeneities that can generate a lot of turbulence, such as vegetation [Yager and Schmeeckle, 2013;Nepf , 2018, 2019], are also not considered in the definition of τ In t . While τ In t is usually not measured for turbulent fluvial bedload transport (it is much below the Shields curve [Paintal, 1971]), it has often been measured in wind tunnel experiments (briefly reviewed in section 3.3.1), including those that sought to determine the initiation threshold of aeolian saltation transport. The reason is that as soon as the first particles of the initially quiescent bed surface are entrained (i.e., begin to roll as rolling requires the smallest flow forces), the flow is usually nearly sufficient to net accelerate them during their downstream motion, resulting in larger and larger particle hops (i.e., the initiation threshold of aeolian saltation transport is only slightly larger than τ In t ) [Bagnold, 1941;Iversen et al., 1987;Burr et al., 2015]. This occurs because, for typical wind tunnels, τ In t is significantly above the cessation threshold of saltation transport (see section 4.3). However, it will become clear that this statement may not apply to aeolian field conditions. Section 3.3.2 briefly reviews models of τ In t derived from wind tunnel experiments, while section 3.3.3 reviews recent evidence that indicates that such models, in general, are unreliable, particularly when applied to field conditions.

Wind Tunnel Experiments of the Initiation of Aeolian Rolling and Saltation Transport
Two distinct experimental setups have been used to measure τ In t . In the first setup, small isolated patches of particles are placed at the bottom of a wind tunnel and then the fluid shear stress τ is increased until particles in such patches start to roll or detach [Williams et al., 1994;Merrison et al., 2007;de Vet et al., 2014]. In the second setup, a complete bed of particles is prepared at the tunnel bottom and then the fluid shear stress τ is increased until saltation transport begins [e.g., Bagnold, 1937;Chepil, 1945;Lyles and Krauss, 1971;Iversen et al., 1976;Greeley et al., 1976Greeley et al., , 1980Greeley et al., , 1984Gillette et al., 1980;Greeley and Marshall, 1985;Nickling, 1988;Iversen and Rasmussen, 1994;Dong et al., 2003;Cornelis and Gabriels, 2004;Burr et al., 2015;Carneiro et al., 2015;Swann et al., 2020] (see also Raffaele et al. [2016, and references therein]). It is worth noting that, according to the definition of τ In t , beginning saltation transport refers to the mere occurrence of saltation transport, even if very sporadic, which is also the definition used by Bagnold [1937]. However, many experimental studies defined beginning saltation transport through a critical loosely defined saltation transport activity (similar to the definition of the fluvial transport thresholds compiled in the Shields diagram), which yields slightly larger threshold values [Nickling, 1988].

Open Problem: Qualitative Discrepancy Between Threshold Measurements
For cohesionless particles (d 100 µm), existing threshold measurements based on the second setup show that τ In t increases relatively strongly with the particle diameter d [Raffaele et al., 2016]. In contrast, for the first setup, measurements indicate that τ In t remains constant with d for d 100 µm [Merrison et al., 2007;de Vet et al., 2014]. The reason for this qualitative inconsistency is not understood. Merrison et al. [2007] suggested that the initiation of rolling (measured in their experiments) may be different to that of saltation transport. However, this suggestion is inconsistent with the observation that saltation transport in wind tunnels is preceded by rolling further upwind [Bagnold, 1941;Iversen et al., 1987;Burr et al., 2015]. Furthermore, in contrast to standard wind tunnel experiments, for experiments in pressurized wind tunnels with Venusian air pressure, both an equilibrium rolling (lower initiation threshold) and an equilibrium saltation transport regime (higher initiation threshold) exist, and both initiation thresholds strongly increase with d [Greeley and Marshall, 1985].

Models of the Initiation of Aeolian Rolling and Saltation Transport
Nearly all existing models of the initiation of aeolian rolling and saltation transport (including sand transport [Bagnold, 1941;Iversen et al., 1976Iversen et al., , 1987Iversen and White, 1982;Shao and Lu, 2000;Cornelis and Gabriels, 2004;Lu et al., 2005;Claudin and Andreotti, 2006;Kok and Renno, 2006;Merrison et al., 2007;Durán et al., 2011;Duan et al., 2013a;de Vet et al., 2014;Burr et al., 2015;Edwards and Namikas, 2015], drifting snow [Schmidt, 1980;Lehning et al., 2000;He and Ohara, 2017], and the transport of regolith dust by outgassed ice on the comet 67P/Churyumov-Gerasimenko [Jia et al., 2017]) predict τ In t from the balance between aerodynamic forces and/or torques and resisting forces and/or torques acting on a bed particle. Even though many of these models do not consider peaks of the aerodynamic force, and some of them do not treat τ In t as what it is (i.e., the threshold at which the fluid entrainment probability exceeds zero, see above), they are conceptually very similar and mainly differ in the empirical equations that they use for the aerodynamic and cohesive interparticle forces. For this reason, only one of the most popular and simple models, the model by Shao and Lu [2000], is discussed here. It reads where A N = 0.0123 is an empirical scaling factor and γ C = 3 × 10 −4 kg/s 2 an empirical constant that accounts for cohesive interparticle forces. More complex models [e.g., Iversen and White, 1982;Claudin and Andreotti, 2006;Durán et al., 2011] involve additional dependencies of Θ In t on the shear Reynolds number Re * or, equivalently, on the Galileo number Ga ≡ (s − 1)gd 3 /ν f ≡ Re * / √ Θ (also called Yalin parameter [Yalin, 1977]).

Effects of the Boundary Layer Thickness on the Initiation of Aeolian Rolling and Saltation Transport
The size of turbulent flow eddies, and thus the duration of turbulent fluctuation events, is limited by the system dimensions, more specifically, the boundary layer thickness δ [see review by Smits et al., 2011, and references therein]. However, in most wind tunnel experiments and the field, the produced turbulent boundary layer should be so thick that any turbulent fluctuation has a nonzero probability to last sufficiently long for entrainment to occur . That is, the mere existence of aerodynamic force peaks that exceed resisting forces is sufficient for τ In t to be exceeded. However, this is no longer true when δ becomes too small, at which point turbulent fluctuation events may cause particles to rock (i.e., vibrate or wobble or oscillate) within their bed pockets but no fluctuation lasts long enough for the particles to completely leave them.  physically modeled such situations and derived an expression for the ratio between τ In t and the shear stress threshold τ In t of incipient rocking (equivalent to the Shields number ratio Θ In t /Θ In t ). These authors' derivation uses the impulse criterion of section 3.2.1 (even though  start with the energy criterion, their analysis is effectively equivalent to assuming a constant impulse threshold) and the fact that the maximal duration T max of turbulent fluctuation events is controlled by δ and the local mean flow velocity u via T max ∝ δ/u [Alhamdi and Bailey, 2017]. The derived expression reads where α f ≡ u m /u ≥ 1 is the ratio between the characteristic flow velocity u m associated with the largest positive fluctuations and u, and f (G) is a factor that encodes information about particle shape, orientation, and the pocket geometry. Equation (11) encompasses three different regimes. In one extreme, if there is a nonzero probability that turbulent fluctuation events associated with the largest positive fluctuations last sufficiently long for particle entrainment, then there will be a nonzero probability that incipient rocking evolves into incipient rolling (i.e., Θ In t /Θ In t 1). In the other extreme, if all positive fluctuation events always last too short, the mean flow must exceed the torque balance for entrainment to occur (i.e., Θ In t /Θ In t α 2 f ). In the intermediate regime between these two extremes, Θ In t /Θ In t is proportional to the square of the inverse dimensionless boundary layer thickness (d/δ) 2 . Although weak logarithmic dependencies on δ/d are also incorporated in α f and Θ In t [Lu et al., 2005], they are dominated by this proportionality. In fact, Figure 7 shows that the prediction for the intermediate regime is roughly consistent with the experimental data by Williams et al. [1994] if one uses that the Shields number for incipient rocking (Θ In t ) is approximately constant, neglecting the logarithmic dependency of Θ In t on δ (and further minor dependencies on Ga). Williams et al. [1994] set up their wind tunnel in a manner that produces a relatively thin developing turbulent boundary layer (i.e., δ increases with downstream distance). However, once the intermediate regime is exceeded (i.e., Θ In t Θ In t ) because δ becomes too (only for illustration purposes as the actual values of Θ In t in the experiments by Williams et al. [1994] are unknown). It is suspected that the one extreme outlier for d = 165 µm may either have been a faulty measurement or be associated with the observation that the boundary layer for this particular sand sample was not always fully turbulent [Williams et al., 1994]. large, as for most wind tunnel experiments with fully developed boundary layers, the logarithmic dependency of Θ In t on δ/d via Θ In t may become significant ( Figure 8). For example, for the same Galileo number Ga, the threshold values measured by Burr et al. [2015] in Figure 8, which were carried out in a pressurized wind tunnel with δ ≈ 1.9 cm, are significantly larger than those measured by Iversen et al. [1976], which were carried out in a wind tunnel with δ ≈ 1.2 m.

Open Problem: Unexpected Behavior of Saltation Transport Initiation Threshold for Large Density Ratio
The very recent measurements by Swann et al. [2020], who used a very-low pressure wind tunnel and three different beds of cohesionless particles (d = [310, 730, 1310] µm) to mimic Martian conditions, indicate that Θ In t unexpectedly increases substantially with Ga and thus d (Figure 8). A possible explanation could be that, because of the very large density ratio s, some of the experimental conditions may have been in the intermediate regime (i.e., 1 ≤ C ≤ α f in equation (11)), in which Θ In t scales with d ( Figure 7). In fact,  Figure 7).  Burr et al., 2015;Swann et al., 2020] versus the Galileo number Ga. The color indicates the thickness of the boundary layer δ relative to the particle diameter d, which controls the relative amplitude of turbulent fluid velocity fluctuations for a constant Ga. Circles correspond to threshold values obtained from the raw data by Swann et al. [2020]. The threshold values for the experiments by Iversen et al. [1976] are found in Iversen and White [1982].

Controversy: Dependency of Saltation Transport Initiation Threshold on Density Ratio
Based on comparisons between experiments in pressurized wind tunnels with comparably very thin boundary layers but larger-than-normal air density [Greeley et al., 1984;Burr et al., 2015] and nonpressurized wind tunnels with comparably very large boundary layers  (and normal air density), Iversen et al. [1987] and Burr et al. [2015] argued that there is an underlying decrease of the saltation transport initiation threshold (which is slightly larger than Θ In t for aeolian transport in typical wind tunnels, see above) with the density ratio s for a constant shear Reynolds number Re * (equivalent to a constant Ga). However, this dependency on s may be an artifact of huge differences in the dimensionless boundary layer thickness δ/d . In fact, even though the dependency of Θ In t on δ/d is logarithmic once the intermediate regime is exceeded (like for the measurements in question), such weak dependencies can still have significant effects once differences in δ/d become very large. This point of view is supported by Figure 8, in which δ/d is color-coded. It can be seen that the yellow, open diamond (a measurement from a nonpressurized wind tunnel) exhibits a similar value of s as the blue symbols (measurements from a pressurized wind tunnels), which was achieved by using a very light particle material (ρ p = 210 kg/m 3 ). Nonetheless, the threshold Θ In t of the former is significantly smaller than those of the latter. Also, the former measurement relatively smoothly connects to the other measurements carried out in the same nonpressurized wind tunnel, which exhibit much larger values of s. On the other hand, the measurements by Swann et al. [2020], for which s is comparably very large and δ/d of a similar size as for the measurements by Iversen et al. [1976], support the density ratio hypothesis because of comparably small values of Θ In t . Note that, for the discussion of threshold values, one has to keep in mind that threshold measurements are highly prone to measurement errors of various sources [Raffaele et al., 2016].
Such errors are likely much larger than often reported because measurements of Θ In t can vary by more than a factor of 2 for a given condition, even for cohesionless particles [Raffaele et al., 2016].

Open Problem: Aeolian Bedload Transport in the Field
In wind tunnel experiments, rolling is being initiated at threshold values that are significantly above the cessation threshold of saltation transport (see section 4.3). This is why rolling seems to always evolve into saltation transport (i.e., equilibrium rolling and thus aeolian bedload transport does not seem to exist) [Bagnold, 1941;Iversen et al., 1987;Burr et al., 2015]. However, atmospheric boundary layers are several orders of magnitude thicker than those of wind tunnels [Lorenz et al., 2010;Petrosyan et al., 2011;Kok et al., 2012;Lebonnois et al., 2018] and may therefore exhibit a significantly smaller rolling threshold. In contrast, the cessation threshold of saltation transport is predominantly a property of the mean turbulent flow (see section 4.3) and therefore rather insensitive to the boundary layer thickness δ. Hence, for atmospheric boundary layers, it is possible that equilibrium rolling transport exists. Note that equilibrium rolling transport has been observed in pressurized wind tunnels with Venusian air pressure for a narrow range of Shields numbers Θ [e.g., Greeley and Marshall, 1985].

Open Problem: Reliable Models of the Initiation Threshold of Planetary Saltation Transport
The most widely used models for the initiation of aeolian saltation transport (see section 3.3.2), which have been adjusted to wind tunnel measurements, do not take into account the dependency of the relative magnitude of turbulent fluctuations on the dimensionless boundary layer thickness δ/d. This may be the reason why these models, when applied to Martian atmospheric conditions, predict threshold shear stresses for fine sand particles that are so large that transport should occur only during rare strong Mars storms [Sullivan and Kok, 2017], in contradiction to modern observations indicating widespread and persistent sediment activity [Bridges et al., 2012a,b;Silvestro et al., 2013;Chojnacki et al., 2015], even of very coarse sand [Baker et al., 2018]. For example, for the Martian conditions reported by Baker et al. [2018] (ρ p = 2900 kg/m 3 , ρ f = 0.02 kg/m 3 , g = 3.71 m/s 2 , d = 1.5 mm), equation (10) predicts for the threshold shear velocity: u In * t ≡ Θ In t (ρ p /ρ f − 1)gd 3.7 m/s, which corresponds to winds that are more than twice as fast as the strongest Mars storms. Note that Lu et al. [2005] proposed a model for the initiation of rolling that includes the effect of δ/d. The authors of this review therefore recommend to use the model by Lu et al. [2005] in combination with models of the cessation threshold of saltation transport (see section 4.3) for the estimation of the occurrence of saltation transport in real atmospheric boundary layers. However, it remains to be demonstrated that this approach yields reliable predictions. In fact, in the field, atmospheric instability, topography gradients, and surface inhomogeneities, such as obstacles and vegetation, can dramatically enhance local turbulence and thus fluid entrainment. Likewise, sublimation of subsurface ice in cold environments (the so-called solid-state greenhouse effect [Kaufmann et al., 2006]) can generate airborne particles of carbon dioxide, methane, and nitrogen ice [Hansen et al., 1990;Thomas et al., 2015;Jia et al., 2017;Telfer et al., 2018]. Given that even a few entrained particle can result in fully developed saltation transport provided that the fetch is sufficiently long [Sullivan and Kok, 2017], it may well be that saltation transport in the field can almost always be initiated close to the cessation threshold [Sullivan and Kok, 2017;Telfer et al., 2018]. Evidence for this hypothesis is seen on Pluto, where aeolian dunes and wind streaks have been observed even though saltation transport initiation had been thought to be virtually impossible because of Pluto's very thin atmosphere (pressure P = 1 Pa) and relatively weak 10 m winds (u max 10m ≈ 10 m/s) [Telfer et al., 2018].

Open Problem: Lack of Direct Aeolian Sediment Transport Initiation Measurements in the Field
The overarching problem associated with the rather poor current knowledge of aeolian sediment transport initiation in the field (see open problems above) is that, to the authors' knowledge, there are no direct field measurements of the transport initiation threshold Θ In t . In fact, existing field experiments have focused on detecting aeolian saltation transport [Barchyn and Hugenholtz, 2011, and references therein] rather than on how the fluid entrainment of individual bed particles, which usually starts out as a rolling motion, leads to saltation transport. Hence, we currently do neither know the wind speeds that are required in the field to initiate rolling transport of individual bed particles nor whether such rolling transport, like in wind tunnels, always evolves into saltation transport (see open problems above). What adds to the problem is that existing field studies either obtain saltation transport threshold estimates using methods that do not seek to distinguish saltation transport initiation and cessation [Barchyn and Hugenholtz, 2011, and references therein] or assume that Θ In t coincides with the continuous saltation transport threshold  (which is a controversial assumption, see section 4.1.3).

The Role of Particle Inertia in Nonsuspended Sediment Transport
As discussed in section 1, old experimental studies [e.g., Ward, 1969;Graf and Pazis, 1977] strongly indicated that the fluvial transport threshold measurements that are compiled in the Shields diagram are to a nonnegligible degree affected by particle inertia. As the Shields diagram shows a rough data collapse of the threshold Shields number Θ t as a function of the shear Reynolds number Re * , this raises the question of whether Re * is in some way associated with particle inertia. Indeed, while Re * has usually been interpreted as the ratio between the particle size and the size of the viscous sublayer of the turbulent boundary layer, Clark et al. [2017] showed that it can also be interpreted as a number that compares the viscous damping time scale to the ballistic time scale between bed collisions. Importantly, these authors showed that the shape of the Shields curve can be partly explained by the fact that inertial particles at high Re * are harder to stop.
In general, the role of particle inertia in nonsuspended sediment transport can be twofold. On the one hand, entrainment by or supported by particle-bed impacts may be able to supply the transport layer with bed particles and thus compensate captures of transported particles by the bed (section 4.1). This mechanism gives rise to a shear stress threshold associated with impact entrainment. On the other hand, although the mean turbulent flow is usually too weak to initiate transport (which instead usually requires turbulent fluctuation events, see section 3), it may be able to sustain the motion of particles that are already in transport. This mechanism gives rise to a physical process-based definition of transport capacity and a shear stress threshold, which has often been misidentified as an entrainment threshold by Shields [1936] and others (section 4.2). Various models for both shear stress thresholds that have been proposed in the literature are compared with one another in section 4.3. Bagnold [1941] was the first to recognize that impact entrainment is crucial for sustaining aeolian saltation transport. Based on his wind tunnel and field observations, he explained [Bagnold, 1941, p. 102], "In air, the grains, when once set in motion along the surface, strike other stationary grains, and either themselves bounce high (a distance measured in hundreds if not thousands of grain diameters) into the relatively tenuous fluid, or eject other grains upwards to a similar height." Largely because of Bagnold's observations, the statistics of particle impacts onto a static granular packing have been subject of many experimental and theoretical investigations (section 4.1.1). Bagnold [1941, p. 102] also believed that impact entrainment is negligible for fluvial bedload transport: "If the physics of this impact-ejection mechanism is applied to sand in water, it is found that the impact momen-tum of the descending grains is insufficient to raise surface grains to a height greater than a small fraction of one grain diameter." However, Bagnold, and numerous researchers after him, did not consider that even a marginal uplift of a bed particle can make it much easier for a turbulent fluctuation event to entrain it (section 4.1.2) and that, once bedload transport becomes sufficiently strong, multiple particle-bed impacts occur in so short sequence that the bed can no longer be considered as static. In fact, for continuous transport, recent studies revealed that impact entrainment alone can sustain bedload transport (section 4.1.3).

Impact of an Incident Particle Onto a Static Granular Packing
The collision process between an incident particle and a static granular packing has been investigated in many experimental [Mitha et al., 1986;Werner, 1990;Rioual et al., 2000Rioual et al., , 2003Tanaka et al., 2002;Nishida et al., 2004;Beladjine et al., 2007;Oger et al., 2008;Ammi et al., 2009;Clark et al., 2012Clark et al., , 2015bClark et al., , 2016Bachelet et al., 2018;Chen et al., 2019] and theoretical [Werner and Haff , 1988;Haff , 1988, 1991;Haff and Anderson, 1993;McElwaine et al., 2004;Oger et al., 2005Oger et al., , 2008Zheng et al., 2005Zheng et al., , 2008Namikas, 2006;Crassous et al., 2007;Bourrier et al., 2008;Kok and Renno, 2009;Valance and Crassous, 2009;Ho et al., 2012;Duan et al., 2013b;Xing and He, 2013;Comola and Lehning, 2017;Huang et al., 2017;Tanabe et al., 2017; studies in order to better understand aeolian saltation transport and other geophysical phenomena (e.g., rockfall [Bourrier et al., 2008;Bachelet et al., 2018]); see also [White and Schulz, 1977;Rice, 1986, 1989;McEwan et al., 1992;Nalpanis et al., 1993;Rice et al., 1995Rice et al., , 1996Dong et al., 2002;McElwaine et al., 2004;McKenna Neuman, 2009, 2011] for collision statistics during ongoing aeolian saltation transport. In typical experiments, a spherical incident particle of diameter d and mass m is shot (e.g., by an airgun) at a given speed v i and angle θ i onto a static packing of spheres of the same size. As shown in Figure 9 and sketched in Figure 10, as a result of its impact on the packing, the incident particle may rebound (velocity v r , angles θ r , φ r ) and/or eject bed particles into motion (number N e , velocity v e , angles θ e , φ e ), where a particle is typically counted as ejected if its center is lifted by more than d above the top of the bed surface. The statistics of this process has been the subject of several recent experimental and numerical studies [e.g., Beladjine et al., 2007;Ammi et al., 2009;Tanabe et al., 2017] (note that experimental studies that used only one camera measured quantities projected into the incident plane: v 2D r(e) ≡ v 2 r(e)x + v 2 r(e)z and tan θ 2D r(e) ≡ tan θ r(e) /cos φ r(e) ). These studies have yielded the following insights: (i) The incident particle loses much more energy in head-on than in grazing collisions. In fact, the average restitution coefficient and its two-dimensional projection obey the following empirical relationships for 10 • ≤ θ i ≤ 90 • : where the overbar denotes an ensemble average over collision experiments, and the A and B coefficients are empirical constants that vary slightly between the studies (e.g., A ≈ A 2D ≈ 0.87, B ≈ 0.62 [Ammi et al., 2009], and B 2D ≈ 0.72 ).
(ii) The average vertical restitution coefficient exceeds unity at small impact angles and obeys the following empirical relationship for 10 • ≤ θ i ≤ 90 • : where A z ≈ 0.3 and B z ≈ 0.15 for the experiments by Beladjine et al. [2007]. Pähtz et al. [2020] suggested the following modification of equation (13): This modification, which is also consistent with the experimental data, ensures the correct asymptotic behavior of the average rebound angle, θ r ∼ √ θ i , in the limit θ i → 0.
(iii) The average rebound angle and its two-dimensional projection are independent of the incident speed, increase with the impact angle, and obey the following empirical relationships for 10 • ≤ θ i ≤ 90 • : where θ 0 ≈ 20 • and χ ≈ 0.19 for the experiments by Ammi et al. [2009].
(iv) The average energy that the incident particle transfers to the bed is spent for the ejection of bed particles. That is, it is proportional to the average of the sum of the kinetic energy of ejected particles (E e = 1 2 mv e 2 and E 2D e = 1 2 mv 2D2 e ). In fact, the following empirical relationships are obeyed for 10 • ≤ θ i ≤ 90 • : where r ≈ 0.04 and r 2D ≈ 0.038 for the experiments by Ammi et al. [2009]. Note that r and r 2D decrease with the coordination number of the particle packing [Rioual et al., 2003].
(v) The average number of ejected particles is a linear function of the incident speed for 10 • ≤ θ i ≤ 90 • : where n 0 ≈ 13 and ζ ≈ 40 for the experiments by Ammi et al. [2009]. Note that n 0 decreases with the coordination number of the particle packing [Rioual et al., 2003].
(vi) The average horizontal and lateral velocities of ejected particles are nearly independent of the incident velocity, but the average vertical velocity increases slightly with the incident velocity and is independent of the impact angle for 10 • ≤ θ i ≤ 90 • [Ammi et al., 2009]: (vii) The average ejection angle θ e is constant for 10 • ≤ θ i ≤ 90 • [Ammi et al., 2009]. However, its projection into the incident plane increases with the impact angle : Open Problem: Behavior of the Rebound Probability Mitha et al. [1986] measured that about 94% of all impacting particles are not captured by the bed (i.e., they successfully rebound). However, the range of impact velocities in their experiments was very narrow (|v i | ∈ (106, 125) √ gd). More systematic measurements of the rebound probability P r are needed.
Studies have attempted to physically describe both the rebound [Zheng et al., 2005[Zheng et al., , 2008Namikas, 2006; and ejection dynamics [McElwaine et al., 2004;Crassous et al., 2007;Kok and Renno, 2009;Valance and Crassous, 2009;Ho et al., 2012;Comola and Lehning, 2017;. For example, the rebound dynamics can be analytically calculated for an idealized packing geometry and a given rebound location assuming a binary collision between the incident particle and hit bed particle. From averaging over all possible rebound locations, one can then determine the rebound angle and restitution coefficient distributions. Using this procedure,  derived the following expressions for e 2D , e z , θ 2D r , and P r in the limit of shallow impact angles (θ i 20 • ): where α r and β r are the normal and tangential rebound restitution coefficients, respectively, in the impact plane, which depend on the binary normal and tangential restitution coefficient (i.e., the ratio between the postcollisional and precollisional relative particle velocity component normal and tangential, respectively, to the contact plane). agreement with the data with θ i 20 • is acceptable considering that the theory has been derived mostly from first physical principles. Equation (23), which is the modified version of equations (41) and (42) of  that these authors describe in their text, cannot be tested because of the lack of systematic measurements of the rebound probability P r . A widely used alternative expression for P r was given by Anderson and Haff [1991]: P r ≈ 0.95[1 − exp(−γ r |v i |)]. However, this expression is empirical and contains the dimensional parameter γ r (note that Andreotti [2004] assumed γ r ∝ 1/ √ gd). Because e z 1, which is a precondition for sustained aeolian saltation transport (from energy conservation), requires shallow impact angles, equations (20)-(23) can be used for the theoretical modeling of aeolian saltation transport.
For the description of the ejection dynamics, there have been two distinct approaches: solving an underdetermined momentum and/or energy balance of the particles involved in the collision process [Kok and Renno, 2009;Comola and Lehning, 2017] and treating the collision process as a sequence of binary collisions, in which the energy is split between the collisional partners (i.e., incident and bed particle or two bed particles) [McElwaine et al., 2004;Crassous et al., 2007;Valance and Crassous, 2009;Ho et al., 2012;. A minimal numerical model that is based on the latter approach has been able to reproduce experimental data of both the rebound and ejection dynamics, including the measured log-normal distribution of the vertical ejection velocity [Crassous et al., 2007]. Furthermore, based on this approach and the derivation by Ho et al. [2012],  derived the following analytical expression for the distribution of the ejection energy E e : from which they further obtained expressions for N e , E e , and |v e |: where r = 0.06. Figure 12 shows that these expressions are roughly consistent with experimental data considering that they have been derived mostly from first physical principles. Note that equations (20) Solid lines correspond to equations (25) and (27) combined with the approximation |v e | |v 2D e |.
uations in which the size of the impacting particle differs from the size of the particles of the granular packing . Further note that equation (38) of , which is the equivalent of equation (27), contains a typo (a σ is missing in the denominator).

Open Problem: Impacts Onto Mobile Beds
The findings from collision experiments with static beds are often applied to model fluvial bedload [Berzi et al., 2016;Pähtz et al., 2020] and aeolian saltation transport [Andreotti, 2004;Claudin and Andreotti, 2006;Creyssels et al., 2009;Kok and Renno, 2009;Kok, 2010a;Jenkins et al., 2010;Lämmel et al., 2012;Pähtz et al., 2012;Huang et al., 2014;Valance, 2014, 2018;Zheng, 2014, 2015;Berzi et al., 2016Berzi et al., , 2017Bo et al., 2017;Pähtz et al., 2020]. However, if the time between successive particle-bed impacts is too short for a bed particle to fully recover from each impact, it can accumulate more and more kinetic energy with each impact. Hence, for a sufficiently large impact frequency and impact energy (both increase with the sediment transport rate Q), the bed can no longer be treated as static and the findings from such collision experiments may no longer apply. For example, the simultaneous impact of two particles onto the bed leads to a significantly different outcome compared with the situation in which each particle impacts separately [Duan et al., 2013b]. For these reasons, future studies should try to systematically investigate the effects of disturbances from the static bed on the outcome of a particle-bed impact.

Open Problem: Effects of Particle Shape and Size Distribution
Chen et al.
[2019] investigated the particle-bed collision process for natural sand particles, which exhibit nonspherical shapes and nonuniform particle size distributions. They found significant quantitative and qualitative deviations from the laws describing spherical, uniform particles. More systematic experimental studies are needed to pinpoint the exact manner in which particle shape and size distribution affect the collision process.

Controversy: Effects of Viscous Damping
Binary collisions that occur within an ambient fluid can be significantly damped depending on the Stokes number St ≡ s|v r |d/(9ν f ) [Gondret et al., 2002;Yang and Hunt, 2006;Schmeeckle, 2014;Maurin et al., 2015], where v r is the relative particle velocity just before a collision. For example, experiments suggest that the effective normal restitution coefficient of a damped binary collision vanishes for St 10 [Gondret et al., 2002]. The question that then arises is how does viscous damping affect the rebound and ejection dynamics of a particle-bed impact. Berzi et al. [2016Berzi et al. [ , 2017 assumed that the rebound restitution coefficients e 2D and e z , like , also vanish when St falls below a critical value. In contrast, DEMbased simulations indicate that the dynamics of saltation [Pähtz and Durán, 2018a] and particularly bedload transport [Drake and Calantoni, 2001;Maurin et al., 2015;Elghannay and Tafti, 2017;Durán, 2017, 2018a,b] are not much affected by the value of , which suggests that the rebound and ejection dynamics of a particle-bed impact may not be much affected by viscous damping. A possible explanation for this unexpected behavior could be that a nearly elastic particle-bed impact may be roughly equivalent to a sequence of binary collisions between particles in contact at the instant of impact. In fact, a theoretical model based on this hypothesis reproduced experiments of the collision process [Crassous et al., 2007;Valance and Crassous, 2009]. For the perfectly elastic case ( = 1), the impactor would then transfer all of its momentum in the direction normal to the contact plane to the particle it hit (which is the expected result of an elastic binary collision between a mobile and a resting particle) and, therefore, rebound with zero normal momentum. A complete loss of normal momentum is also expected for the completely inelastic case ( = 0). This suggests that the rebound process is not much affected by , which would imply that the momentum in the direction tangential to the contact plane is what mainly matters. Collision experiments in an ambient viscous liquid could resolve this controversy.

Open Problem: Effects of Cohesion
Cohesive interparticle forces, including van der Waals [Castellanos, 2005], water adsorption [Herminghaus, 2005], and electrostatic forces [Lacks and Sankaran, 2011], become significant in the collision process for sufficiently small particles (on Earth, for d 100 µm) because they scale with a lower power p in the particle diameter (F coh ∼ d p ) than the gravity force (F g ∼ d 3 ). However, collision experiments with so small particles have not been carried out because it is very difficult to detected their dynamics with cameras. Numerical studies are also very scarce. To the authors' knowledge, only the very recent study by Comola et al. [2019a] studied cohesive forces, by implementing them in a numerical DEM-based model of aeolian saltation transport. These authors investigated the impact of a particle onto the bed for a large range of the strength of cohesive forces and found that cohesion decreases N e via solidifying the bed, while e slightly and |v e |/|v i | considerably increase. However, more systematic studies are needed to confirm these results and determine scaling laws describing the effects of cohesion on the outcome of a particle-bed impact.

Collision-Enhanced Turbulent Entrainment in Fluvial Bedload Transport
To the authors' knowledge, only a single study has resolved the effects of particle-bed impacts on entrainment by turbulent fluctuation events in bedload transport [Vowinckel et al., 2016]. However, this study provided one of the largest, if not the largest, data sets of entrainment events associated with fluvial bedload transport with a very high resolution in space and time. Vowinckel et al. [2016] coupled direct numerical simulations (DNS) for the fluid phase (i.e., the Navier-Stokes equations are directly solved without using turbulent closure assumptions) with DEM simulations for the particle phase (i.e., particles interact with each other according to a contact model) using the immersed boundary method, which fully resolves the geometry of particles (and thus the hydrodynamic forces acting on them) without remeshing the grid during their motion [Vowinckel et al., 2014]. Because of the sophistication of this numerical method (i.e., resolving all relevant physical processes at very small scale), the produced data can be considered to be very reliable. The simulated setup consisted of two layers of grains resting on the simulation bottom wall, the lower of which was fixed, arranged in a hexagonal packing, and exposed to a unidirectional open channel flow of thickness H = 9d (Reynolds number Re ≡ U b H/ν f = 2941, where U b is the bulk flow velocity). The Shields number was at Θ = 0.0255, which is about 25% below the Shields curve for the simulated condition. That is, the nondimensionalized transport rate Q * was likely below the value associated with critical transport conditions (see section 1), which is consistent with Vowinckel et al. [2016] reporting that only 3% of all particles were in motion on average. For these conditions, it was found that, in the vast majority of cases (overall 96.5%), a particle-bed impact and a subsequent turbulent fluctuation event are responsible for entrainment, even when one or more of the six pockets surrounding the target particle were not occupied by other particles (in which case the target particle experiences a larger exposure to the flow). For an entrainment event following this pattern, Figure 13 shows the time evolution of (a) the vertical displacement (y p ) and (b) velocity of a bed surface particle (u p ), while Figure 14 shows the simulation domain and contour plots of the instantaneous flow field. It a b can be seen that, at the instant of entrainment, the instantaneous streamwise flow velocity (u) exhibits larger-than-average values (Figure 14c). In fact, Vowinckel et al. [2016] reported that 82% of the entrainment events were caused by sweep, characterized by positive fluctuations The results by Vowinckel et al. [2016] were obtained for an idealized hexagonal packing and may not necessarily apply in their full extent to realistic sediment beds found in nature. While for a hexagonal packing, the vast majority of entrainment events are initiated by particle-bed impacts, it remains unclear whether this holds true also for natural sediment beds, in which bed surface particles tend to protrude much more strongly into the flow. On the one hand, a larger protrusion makes it easier for a turbulent fluctuation event to entrain a bed surface particle without a preceding particle-bed impact. On the other hand, particlebed impacts can result in entrainment without the need of a turbulent fluctuation event (see section 4.1.3). Pähtz and Durán [2017] numerically studied the role of particle-bed impacts in sustaining continuous nonsuspended sediment transport for transport conditions characterized by a large range of the Shields number Θ, density ratio s, and Galileo number Ga. These authors coupled quasi-two-dimensional DEM simulations for the particle phase with a Reynoldsaveraged description of the fluid hydrodynamics that neglects turbulent fluctuations around the mean turbulent flow. While such simulations cannot resolve entrainment by turbulent fluctuation events, they are able to elucidate the importance of entrainment by particle-bed impacts relative to entrainment by the mean turbulent flow. Also, the absence of turbulent fluctuations eliminates transport intermittency in the sense that transport in the simulation domain is either continuous (i.e., periods of rest are absent) or it completely stops after a finite time (except for potential creeping, see section 2.3). From their simulations, Pähtz and Durán [2017] determined an effective value of the local particle velocity averaged over elevations near the bed surface (V b ) relative to the critical velocity that is needed to escape the potential wells set by the pockets of the bed surface (∝ √ĝ d, whereĝ = [1+1/(s +C m )] is the value of the gravity constant reduced by the buoyancy and added mass force, with C m = 1/2 the added mass coefficient). They found that V b / √ĝ d exhibits a universal approximately constant value of order unity for continuous nonsuspended sediment transport if the following constraint is obeyed: Im ≡ Ga s + C m 20 or Θ 5/Im.

The Role of Particle-Bed Impacts in Sustaining Continuous Sediment Transport
The interpretation of V b / √ĝ d ≈ const is that particles located near the bed surface (which includes both particles of the bed and transported particles) are on average at the verge of leaving it or being captured by its potential wells, consistent with a dynamic equilibrium that is solely controlled by particle inertia. This implies that entrainment occurs solely due to the action of particle-bed impacts. Consistently, Pähtz and Durán [2017] observed from visually inspecting simulations that obey equation (28) that every entrainment event is initiated by a particle-bed impact, usually with a small time delay between the instant of impact and beginning visible motion. In contrast, for transport conditions that do not obey equation (28), V b / √ĝ d exhibits a smaller value, which means that the mean turbulent flow must assist particles located near the bed surface in escaping the potential wells. For bedload transport, the findings by Pähtz and Durán [2017] were independent of the effective normal restitution coefficient for a damped binary collision, which indicates that viscous damping does not suppress impact entrainment (see also the discussion of viscous damping in section 4.1.1).
The constraint set by equation (28) is obeyed by the vast majority of sediment transport regimes, including turbulent fluvial bedload transport. That is, for the absence of turbulent fluctuation events, only viscous fluvial bedload transport is significantly affected by the entrainment of bed sediment by the mean turbulent flow. The numerical prediction that impact entrainment dominates entrainment by the mean turbulent flow in turbulent fluvial bedload transport is supported by experiments [Heyman et al., 2016;Lee and Jerolmack, 2018]. Lee and Jerolmack [2018] studied bedload transport driven by a water flow in a quasi-twodimensional flume (i.e., its lateral dimension was only slightly larger than the particle diameter d). Because the size of turbulent structures, and thus turbulent fluctuation events, is strongly suppressed when the system dimensions are so strongly narrowed down, their experiments are somewhat comparable to the numerical simulations by Pähtz and Durán [2017] described above. Lee and Jerolmack [2018] fixed the water discharge and fed particles at the flume entrance with varying frequency f in (the tested range of f in was likely associated with a transport rate below capacity). In contrast to similar older experiments [Böhm et al., 2004;Ancey et al., 2008;Heyman et al., 2013], the bed was relatively deep, which ensured the complete dissipation of shock waves associated with particle-bed impacts [Rioual et al., 2003]. Lee and Jerolmack [2018] reported that, for all tested conditions, every entrainment event is initiated by a particle-bed impact, exactly as numerically predicted, and that the number of transported particles roughly scales with the energy transferred to the bed by rebounding particles. The latter finding is remarkably similar to the scaling of the average ejected particle number N e in static bed experiments (e.g., see equation (25)). Lee and Jerolmack [2018] also measured the frequency of particles passing an illuminated window near the flume exit ( f out ). They found that f out < f in for sufficiently small f in and that f out ≈ f in once f in exceeds a critical value. Heyman et al. [2016], who used a water flume with a narrow but larger width (W = 5d) than Lee and Jerolmack [2018] and who also used a relatively deep bed. Heyman et al. [2016] measured that the entrainment rate was proportional to the number of transported particles per unit bed area, which is indirect evidence supporting that the majority of entrainment events is caused by particle-bed impacts. These authors also reported for all their tested feeding frequencies f in that the entrainment and deposition rate are equal to one another, in resemblance of the measurement f out ≈ f in for sufficiently large f in by Lee and Jerolmack [2018]. Note that one expects the approximate equality f out ≈ f in to break down for large f in (when the influx exceeds transport capacity) because increasing momentum transfer from fluid to particles slows down the flow, which at some point can no longer sustain the particle motion.

Similar observations were made by
The results by Heyman et al. [2016] and Lee and Jerolmack [2018] suggest that mainly (but not solely) particles that were previously in motion are being entrained by particle-bed impacts. Otherwise, there would be no reason to expect that the entrainment and deposition rate are relatively equal to one another for a large range of f in (instead, one would expect that only for transport capacity). This can be explained when assuming that particle-bed impacts are effective in mobilizing a bed particle almost only when the bed particle exceeds a critical energy level just before the impact. On the one hand, this assumption would explain why bed particles that have never been transported only rarely become mobilized by particle-bed impacts. On the other hand, this assumption is consistent with the fact that a transported particle that has just been captured by a bed pocket exhibits a residual kinetic energy that takes some time to be completely dissipated, during which it can be remobilized by an impact from a particle coming from behind. It seems that, once f in exceeds a critical value, there is usually a particle coming from behind in time and transported particles can only rarely settle completely even though they may temporarily stop. Temporary particle stops and reentrainment make transported particles tend to move in clusters near the flume exit even though they are apart from one another at the flume entrance, which is exactly what Lee and Jerolmack [2018] reported and what can be observed in the numerical simulations by Pähtz and Durán [2018a, Movie S2].
There is evidence that the presumed impact entrainment mechanism described above may play an important role in nonsuspended sediment transport in general. In fact, in simulations of steady, homogenous sediment transport using DEM-based numerical models that neglect turbulent fluctuations around the mean turbulent flow, the steady state transport rate Q exhibits a discontinuous jump at a fluid shear stress τ ImE t [Carneiro et al., 2011[Carneiro et al., , 2013Clark et al., 2015aClark et al., , 2017Pähtz and Durán, 2018a]. That is, for τ ≥ τ ImE t , transport is significantly larger than zero (Q > 0) and continuous, whereas Q 0 when τ < τ ImE t . Assuming that only impact entrainment took place in all these simulations (as the mean turbulent flow is too weak for entrainment, see above), τ ImE t can be identified as the impact entrainment threshold. The discontinuous jump of Q thus means that, in order for impact entrainment to sustain transport, a critical transport rate must be exceeded. Like the critical feeding frequency in the experiments by Lee and Jerolmack [2018], this critical transport rate may be interpreted as the value above which most transported particles can be captured only temporarily by bed pockets as they are usually hit in time and thus reentrained by an impact from a particle coming from behind before dissipating too much of their kinetic energy. However, it is crucial to point out that impact entrainment of bed particles that have never been transported occasionally occurs in DEM-based sediment transport simulations as well, which is why a further interpretation of the physical origin of the discontinuous jump of Q has been proposed [Pähtz and Durán, 2018a]. It states that, at a critical transport rate, bed surface particles do no longer sufficiently recover between successive particle-bed impacts. They thus accumulate energy between successive impacts until they are eventually entrained. In contrast, for subcritical transport rates, particles sufficiently recover between impacts so that impact entrainment is inefficient, causing transport to eventually stop. The two interpretations above are based only on the energy of bed particles or temporarily captured transported particles. In contrast, in the context of an idealized continuous rebound modeling framework (see section 4.2), an alternative mechanism based on the critical amount of energy E c that bed particles need to acquire for entrainment (more precisely, for entering a quasi-continuous motion) can explain the discontinuous jump of Q without further assumptions (see section 4.2.1).

Open Problem: Precise Mechanism of Impact Entrainment in Continuous Transport
The proposed impact entrainment mechanisms described above and in section 4.2.1 are mostly speculative and based on indirect experimental or theoretical evidence, or idealized models. More direct investigations are therefore needed to uncover the precise nature of impact entrainment and the degree to which each of these mechanisms contributes. Such investigations may also help to better understand fluctuations of nonsuspended sediment transport. For example, the longer the average time t conv it takes for transport to stop (in the absence of turbulent fluctuations around the mean turbulent flow) when τ < τ ImE t (t conv obeys a critical scaling behavior at τ ImE t , see section 2.4), the larger are the transport autocorrelations, which can be quite substantial in fluvial bedload transport [Heathershaw and Thorne, 1985;Drake et al., 1988;Dinehart, 1999;Ancey et al., 2006Ancey et al., , 2008Ancey et al., , 2015Martin et al., 2012].

Open Problem: Precise Definitions of Intermittent and Continuous Transport
As explained above, in simulations of steady, homogenous sediment transport using DEM-based numerical models that neglect turbulent fluctuations, the steady transport rate Q (in a time-averaged sense) exhibits a discontinuous jump at the impact entrainment threshold τ ImE t . In contrast, for most natural conditions, fluid entrainment by turbulent events can reinitiate transport whenever it temporarily stops, meaning that Q remains significant below τ ImE t [Carneiro et al., 2011]. Hence, since turbulent events capable of fluid entrainment occur only at an intermittent basis (see section 3), Pähtz and Durán [2018a] suggested that τ ImE t is equivalent to the continuous transport threshold for most natural conditions and that transport becomes intermittent below τ ImE t . However, provided that fluid entrainment does occur, it is certain to find particles being in transport below τ ImE t at any given instant in time in the large-system limit, which renders the distinction between intermittent and continuous transport somewhat ambiguous. For this reason, Pähtz and Durán [2018a] referred to intermittent conditions as those that deviate significantly from transport capacity (defined as in section 4.2.2). Consistently, Martin and Kok [2018] and Comola et al. [2019b] found from aeolian field experiments that the long-term-averaged transport remains at capacity when the fraction f Q of active saltation transport is close to unity, that is, when transport quantified over a short but somewhat arbitrary time interval (2 s  or 0.04 s [Comola et al., 2019b]) almost never stops. Interestingly, Comola et al. [2019b] showed that the value of f Q can be indirectly estimated from the lowpass-filtered wind speed associated with large and very large scale turbulent structures (cutoff frequency Ω ≈ 0.04 Hz). Alternatively, for their coupled DNS/DEM simulations of fluvial bedload transport, González et al. [2017] fitted continuous functions to the distributions of the discrete transported particle number (defined as the number of particles faster than a somewhat arbitrary velocity threshold) at different τ and identified the onset of continuous transport as the value of τ at which these fitting functions predict a zero probability for a vanishing particle number. Future studies should investigate the compatibility of these and other definitions of continuous transport.

Controversy: Threshold of Continuous Aeolian Saltation Transport
In the opinion of the authors, the evidence reviewed above for the hypothesis that continuous transport occurs once impact entrainment alone is sufficient in compensating random captures of transported particles is quite strong. (In other words, significant fluid entrainment may occur in continuous transport-and does so quite likely in aeolian saltation transport given that the turbulent intensity within the saltation transport layer increases with the sediment transport rate [Li and McKenna Neuman, 2012]-but it is not needed to sustain continuous transport.) However, it is worth pointing out that most aeolian researchers prefer a different narrative for aeolian saltation transport. For example, Martin and Kok [2018] assumed that continuous aeolian saltation transport in the field occurs once the saltation transport initiation threshold (≈ τ In t ) is exceeded, whereas the impact entrainment threshold describes the cessation of intermittent saltation transport. This assumption is based on the idea that fluid entrainment continuously provides the transport layer with bed particles. However, this idea is problematic because turbulent events capable of fluid entrainment occur only at an intermittent basis (see section 3). The interested reader is also referred to the commentary by , in which this controversy is extensively discussed.

Continuous Particle Rebounds and Transport Capacity
In order for the mean turbulent flow to sustain the motion of particles that are already in transport, it needs to compensate, on average, the energy dissipated in particle-bed rebounds via drag acceleration during the particle trajectories. This mechanism, which is illus-trated in detail by means of a thought experiment in section 4.2.1, gives rise to a shear stress threshold of sediment transport (henceforth termed rebound threshold), as was already noted by Bagnold [1941, p. 94] for aeolian saltation transport: "Physically [the rebound threshold] marks the critical stage at which the energy supplied to the saltating grains by the wind begins to balance the energy losses due to friction when the grains strike the ground [and rebound]." It also suggests a clear-cut definition of transport capacity, which is otherwise difficult to define [see review by Wainwright et al., 2015, and references therein], that leads to an experimentally and numerically validated universal scaling of the transport load M (i.e., the mass of transported sediment per unit bed area) with the fluid shear stress τ (section 4.2.2). From the appearance of the rebound threshold in this scaling of M, one can conclude that at a significant if not predominant portion of the threshold measurements by Shields [1936] and others have been misidentified as measurements of the entrainment threshold (section 4.2.3).

Particle Rebounds Along a Flat Wall
To illustrate the concept of continuous particle rebounds, the motion of a particle along a flat wall driven by a constant flow (e.g., the mean turbulent flow) is considered. This particle shall never be captured and instead, for illustration purposes, always rebound with a constant angle and lose a constant fraction of its impact energy (the core of the argument will not significantly change if more sophisticated rebound laws, such as equations (20)-(22), are considered). For this idealized scenario, there are two extremes of possible particle trajectories depending on the initial particle velocity v ↑ , which are sketched in Figure 15. First, if the Figure 15. Sketch of continuous rebound mechanism. Depending on its initial kinetic energy E ↑ relative to a critical energy level E c that depends on the properties of the flow, a particle (yellow lines) either (a) gains sufficient energy in its hops along a flat wall (black lines) to approach a steady, periodic hopping motion or (b) net loses energy until it stops. corresponding initial kinetic energy E ↑ exceeds a critical value E c , the particle will spend sufficiently long within the flow so that it gains sufficient energy via fluid drag during its hops to approach a steady, periodic hopping motion, in which its energy gain via fluid drag is exactly balanced by its energy loss during its rebounds (Figure 15a). Henceforth, such particles are termed continuous rebounders. Second, if E ↑ < E c , the particle loses net energy in its initial and all subsequent hops until it stops (Figure 15b). The critical energy E c depends on properties of the flow. Crucially, if the flow is too weak, all possible trajectories fall into the second category (i.e., E c = ∞).
There are a few takeaways from the this simple thought experiment for realistic systems. First, as the mean turbulent flow is controlled by the fluid shear stress τ, it suggests the existence of a rebound threshold τ Rb t below which the energy losses in particle-bed rebounds cannot be compensated by the flow on average [Jenkins and Valance, 2014;Berzi et al., 2016Berzi et al., , 2017Pähtz and Durán, 2018a;Pähtz et al., 2020]. Second, the randomness introduced by inhomogeneities of the bed and turbulent fluctuations of the flow introduce trajectory fluctuations that can lead to random losses of continuous rebounders, particularly when the lift-off energy accidentally falls below E c [Pähtz and Durán, 2018a]. Such losses must be compensated by the entrainment of bed particles into the continuous rebound layer. Hence, the mere mobilization of bed particles is not sufficient because the lift-off energy of mobilized particles must also exceed E c . In particular, for rebound threshold models (see section 4.3), it has been shown that E c becomes equal to the average rebound energy of continuous rebounders in the limit τ → τ Rb t [Pähtz et al., 2020]. This implies that the impact entrainment threshold τ ImE t must be strictly larger than τ Rb t , since the energy of an entrained particle is much smaller than the energy of the particle that caused its entrainment (i.e., a continuous rebounder) because of energy conservation. In particular, τ ImE t > τ Rb t automatically explains the discontinuous jump of the sediment transport rate Q at τ ImE t that has been observed in the absence of fluid entrainment by turbulent fluctuation events (see section 4.1.3) because Q(τ ImE

Transport Capacity Interpretation Based on Continuous Rebounds
A third takeaway for realistic systems of the thought experiment described in section 4.2.1 involves the fact that, because of momentum transfer from flow to particles, the flow slows down with increasing transport load M. Hence, for a given τ > τ Rb t , provided that there is an abundance of impact and/or fluid entrainment, the system tends to entrain bed material until the mean turbulent flow becomes so weak that it can barely sustain the average motion of continuous rebounders [Pähtz and Durán, 2018b]. Any further slowdown of the flow would then spike the deposition rate, leading to a decrease of M and subsequent increase of the flow speed. That is, the system is at a dynamic equilibrium that may be interpreted as transport capacity. Pähtz and Durán [2018b] analytically showed that this interpretation of transport capacity leads to the capacity scaling whereg = (1 − 1/s)g is the buoyancy-reduced value of the gravitational constant g and µ b = τ pb /P b an approximately constant bed friction coefficient (i.e., the ratio between the particle shear stress τ pb and normal-bed particle pressure P b Mg evaluated at the bed surface). Note that the definitions of τ pb and P b (and thus µ b ), in contrast to the definitions of τ p and P (and thus the yield stress ratio µ s , see section 2.1), include contributions from stresses associated with the particle fluctuation motion in addition to contributions from intergranular contacts. The derivation of equation (29) by Pähtz and Durán [2018b] is based on two main steps: showing the approximate constancy of µ b starting from a geometric constraint on particle-bed rebounds in the steady state and assuming τ gb τ − τ Rb t , which expresses the aforementioned dynamic equilibrium condition associated with the continuous rebound motion. Interestingly, τ gb describes the momentum that is transferred from flow to transported particles per unit bed area per unit time, which implies that high-buoyant fluids (smallg), such as water, require a larger transport load M for a given rate of momentum transfer (i.e., for a given Mg ∝ τ gb ) than low-buoyant fluids (largeg), such as air. Pähtz and Durán [2018b] tested these derivation steps with numerical data from DEM-based simulations of nonsuspended sediment transport (the same as those by Pähtz and Durán [2017], see section 4.1.3). It turned out that these steps, and thus equation (29), are obeyed across nonsuspended sediment transport conditions with Ga √ s 10 (all but relatively viscous bedload transport) provided that the bed surface is defined as the effective elevation of energetic particle-bed rebounds.
The functional form of equation (29) is the foundation of the majority of theoretical and experimental shear stress threshold-based expressions for the capacity transport rate, Q Mv x (where v x is the average streamwise velocity of particles moving above the bed surface), and goes back to the pioneering theoretical descriptions of nonsuspended sediment transport by Bagnold [1956Bagnold [ , 1966Bagnold [ , 1973. However, Bagnold's physical interpretation of the assumptions leading to this scaling was inaccurate: µ b is not equal to µ s and τ Rb t is not an entrainment threshold, as Bagnold assumed [Pähtz and Durán, 2018b]. In fact, equation (29) has no association with sediment entrainment whatsoever, except for the fact that sediment entrainment is a necessary requirement to keep transport at capacity [Pähtz and Durán, 2018b].
As explained in section 4.1.3, Q, and thus M, is significantly larger than zero at the impact entrainment threshold τ ImE t . In particular, transport becomes intermittent for τ < τ ImE t or even stops in the absence of entrainment by turbulent fluctuation events (i.e., transport capacity cannot be sustained). Hence, equation (29) is, in general, valid only for τ ≥ τ ImE t and also consistent with the rebound threshold model prediction τ ImE t > τ Rb t (see section 4.2.1). Note that aeolian saltation transport experiments [Carneiro et al., 2015;Martin and Kok, 2018] and coupled DNS/DEM fluvial bedload transport simulation [González et al., 2017], indeed, very roughly suggest τ ImE t ≈ 1.5τ Rb t and τ ImE t ≈ 2τ Rb t , respectively. In order to extend the validity of equation (29), and thus of standard sediment transport rate relationships, to shear stresses τ with τ Rb t < τ < τ ImE t , one must abandon long-term averaging sediment transport data. Instead, it is necessary to conditionally average M (or Q) only over periods of nearcapacity transport (on short-term average), but ignore periods with transport significantly below capacity or even at rest [Bunte and Abt, 2005;Singh et al., 2009;Shih and Diplas, 2018;Comola et al., 2019b]. Likewise, for realistic fluvial bedload transport, it is necessary to exclude the turbulence-driven fluctuation motion (including turbulent entrainment events) when measuring M for equation (29) to remain valid; otherwise, transport does not vanish for τ → τ Rb t . Salevan et al. [2017] demonstrated that implementing such constraints in the analysis of experimental data is, in principle, possible. By separating the velocity distribution of all measurable particles (including those that are visually perceived as resting) into a Student's t-distribution associated with the turbulence-driven fluctuation motion and an exponential distribution associated with the bulk transport of particles (which automatically implies conditional averaging as periods of rest do not affect this distribution), they obtained a measure for the number of transported particles relative to the total number of bed surface particles (n tr /n tot ). This measure, indeed, vanishes within experimental precision below a Shields number threshold (Figure 16a), which can be interpreted as Θ Rb t , whereas the number of particles n v t that are faster than a certain velocity threshold v t remains nonzero for the entire range of Θ because of the turbulence-driven fluctuation motion (Figure 16b).

Does the Shields Diagram Truly Show Incipient Motion Thresholds?
The Shields diagram is a compilation of measurements of the threshold Shields number Θ t as a function of the shear Reynolds number Re * , which have been labeled as measurements of incipient sediment motion by numerous studies and reviews [e.g., Shields, 1936;Miller et al., 1977;Yalin and Karahan, 1979;Parker and Klingeman, 1982;van Rijn, 1984;Wiberg and Smith, 1987;Ling, 1995;Buffington and Montgomery, 1997;Dey, 1999;Paphitis, 2001;Cao et al., 2006;Dey and Papanicolaou, 2008;Ali and Dey, 2016;Dey andAli, 2018, 2019;Yang et al., 2019, and references therein]. However, incipient motion of turbulent fluvial bedload transport is much better characterized by impulse and energy-based criteria (section 3.2), unless one refers to the Shields number Θ In t at which the fluid entrainment probability exceeds zero (section 3.3), which is much below the Shields curve [Paintal, 1971]. Furthermore, in steady, homogenous turbulent fluvial bedload transport in which turbulence is suppressed (e.g., in narrow water flumes), the vast majority of entrainment events is caused by particle-bed impacts (see section 4.1.3). It is therefore here argued, based on the results reviewed in section 4.2.2, that many of the threshold data compiled in the Shields diagram are actually measurements of the rebound threshold Θ Rb t .
The Shields diagram shows two kinds of threshold measurements obtained using two different methods. The first method is the reference method, where one takes paired measure- ments of Θ and the nondimensionalized transport rate Q * (or transport load M * ≡ M/(ρ p d)) and extrapolates them to the Shields number at which Q * (or M * ) either vanishes [e.g., Shields, 1936] (it is slightly controversial whether Shields really used this method [Buffington, 1999]) or equals a small reference value [e.g., Parker and Klingeman, 1982]. This method yields approximately the rebound threshold Θ Rb t if an expression for Q * (or M * ) based on equation (29) is used for the extrapolation and provided that the data used for the extrapolation are at capacity (i.e., Θ ≥ Θ ImE t ). For example, Lajeunesse et al. [2010] extrapolated their measurements (many data points obeyed Θ ≥ 2Θ Rb t ≈ Θ ImE t ) to M * = 0 using exactly equation (29), yielding exactly Θ Rb t . That the reference method yields the rebound threshold Θ Rb t is further supported by the fact that the values of Θ Rb t obtained from the DEM-based fluvial bedload transport simulations by Pähtz and Durán [2018a] are consistent with the compilation of reference method-based threshold measurements by Buffington and Montgomery [1997].
The second method is the visual method, where one increases Θ until criteria defining what is considered critical transport are obeyed [e.g., Kramer, 1935] (see section 1). The threshold values obtained from this method depend significantly on the chosen criterion and are, on average, close to those obtained from the reference method [Buffington and Montgomery, 1997]. For example, the transition point (Θ, Q * ) ≈ (0.05, 0.007) at which the function Q * (Θ) measured in the gravel-bed experiments by Paintal [1971] changed from Q * ∝ Θ 16 to Q * ∝ Θ 2.5 (see section 1) is indistinguishable from the reference threshold for the same conditions within measurement uncertainty. In particular, a close examination of Paintal's and other gravel bed data has revealed that Paintal's power-16 region can actually be subdivided into two regions [Dey and Ali, 2019, Figure 5] (see also [Shih and Diplas, 2019, Figure 8b]): one region (Θ 0.04) with a milder power law and one with a stronger power law (0.04 Θ 0.05), which includes a jump of Q * by an order of magnitude at Θ 0.04. Such a jump is consistent with exceeding the rebound threshold Θ Rb t because transported particles suddenly become able to move along the surface for comparably long times before being captured by the bed. Hence, it seems that also the visual method, at least for typical critical transport criteria, approximately yields the rebound threshold Θ Rb t rather than an entrainment threshold.
The hypothesis that the Shields diagram shows measurements of the rebound threshold is further supported by the fact that certain rebound threshold models [Pähtz and Durán, 2018a;Pähtz et al., 2020] reproduce the Shields curve without fitting to the experimental data compiled in the Shields diagram (see section 4.3), even when limited to only visually measured data [Pähtz et al., 2020].

Sediment Transport Cessation Models
This section reviews theoretical models for both the rebound threshold Θ Rb t and impact entrainment threshold Θ ImE t . One of the early motivations for developing such models was to better understand the hysteresis between the initiation and cessation of aeolian saltation transport observed in wind tunnel experiments [e.g., Bagnold, 1941;Chepil, 1945;Iversen and Rasmussen, 1994;Carneiro et al., 2015]. While the difference between transport initiation and cessation is relatively small on Earth, wind tunnel experiments and observations suggested a substantial difference on Mars, which needed to be explained [Almeida et al., 2008;Kok, 2010b]. (However, note that extrapolating wind tunnel measurements of the initiation threshold Θ In t to field conditions using standard initiation threshold models is actually inappropriate because Θ In t depends on the boundary layer thickness δ, as discussed in section 3.3.) Later on, cessation threshold models were developed with the purpose to unify fluvial bedload and aeolian saltation transport in a single theoretical framework [Berzi et al., 2016;Pähtz and Durán, 2018a;Pähtz et al., 2020].
As cessation threshold models are associated with a sustained motion of transported particles, they require a physical description of the particle motion within the transport layer that is coupled with boundary conditions that describe the interaction between transported particles and the bed surface. In general, there have been two approaches to describe the transport layer and bed interactions. The first approach consists of representing the entire particle motion by particles moving in identical periodic trajectories along a flat wall that mimics the bed surface (section 4.3.1). The second approach consists of deriving general correlations between transport layer-averaged physical quantities and obtain the correlation coefficients from numerical simulations (section 4.3.2). It will be shown that the latter approach is probably a rough approximation of a variant of the former. Correlation-based model equations elucidate the role that the density ratio s plays for the rebound threshold Θ Rb t in a simple manner and therefore provide a simple conceptual explanation for why Θ Rb t is smaller in aeolian saltation than in fluvial bedload transport (section 4.3.3).

Open Problem: Effect of Cohesion on Transport Cessation Thresholds
Most of the sediment transport cessation threshold models reviewed here account for cohesive interparticle forces and do so in a similar manner as transport initiation threshold models. However, Comola et al. [2019a] recently revealed that the effects of cohesion on transport cessation and initiation thresholds are actually fundamentally different from one another, which is why this section only considers versions of existing cessation threshold models for cohesionless particles. The effect of cohesion on transport cessation thresholds remains a major open problem.

Identical Periodic Trajectory Models (IPTMs)
Most studies proposing cessation threshold models start with the assumption that the motion of transported particles can be represented by a system in which all particles hop in the same periodic trajectory, referred to as the average trajectory, driven by the mean turbulent flow along a flat wall, with which they interact according to certain boundary conditions [Claudin and Andreotti, 2006;Kok, 2010a;Berzi et al., 2016Berzi et al., , 2017Pähtz et al., 2020]. (Note that, although Kok [2010a] does not explicitly refer to identical periodic trajectories, his mathematical treatment of the problem is equivalent to IPTMs.) However, the assumption of identical periodic particle trajectories introduces a variety of potentially major weak-nesses, which has cast doubt on the reliability of IPTMs [Andreotti, 2004;Durán, 2017, 2018a]: 1. In IPTMs, the particle concentration increases with elevation z and jumps to zero when z exceeds the hop height [Anderson and Hallet, 1986]. In contrast, in real nonsuspended sediment transport, it monotonously decreases with z, often exponentially [e.g., Durán et al., 2012]. IPTMs that refer only to the motion of a well-defined species of particles (e.g., continuous rebounders) do not necessarily suffer from this weakness because the concentration profile associated with this species may behave differently from that of the entire ensemble of transported particles. 2. In IPTMs, the mean square of the vertical particle velocity ( v 2 z ) decreases with z. In contrast, in real nonsuspended sediment transport, it increases with z, except far from the bed surface [Pähtz and Durán, 2017]. This behavior is a signature of the fact that the transport layer, in general, consists of different species of particles with different characteristic velocities [e.g., Durán et al., 2011, Figure 21]. That is, IPTMs that refer only to the motion of a well-defined species of particles (e.g., continuous rebounders) do not necessarily suffer from this weakness. 3. Only particles that take off from the wall with an energy E ↑ that is larger than a critical value E c can continue their motion after the initial few hops ( Figure 15). That is, IPTMs that take into account the motion of entrained particles [Claudin and Andreotti, 2006;Kok, 2010a] effectively assume that all entrained particles obey E ↑ ≥ E c even though most of them do not [Pähtz and Durán, 2018a]. 4. IPTMs neglect particle motion via rolling and sliding, which is significant in bedload transport.
Depending on the boundary conditions, three conceptually different kinds of IPTMs can be distinguished: 1. Models of the rebound threshold Θ Rb t consider only the dynamics of continuous rebounders. Their rebounds are described, for example, by equations (12b) and (13), which link the streamwise (x) and normal-wall (z) components of the impact velocity v i to the streamwise and normal-wall components of the rebound velocity v r . Such models then look for the smallest Shields number that results in a periodic trajectory under the constraint that the hop height of particles exceeds one particle diameter (z h ≥ d). This constraint ensures consistency with the underlying model assumption that continuous rebounders are never captured by the bed surface. The threshold resulting from this constraint is denoted as Θ Rb * * t . Pähtz et al. [2020] modified this constraint to take into account that the near-surface flow can assist particles in escaping the bed surface and is even predominantly responsible for the escape in the viscous bedload transport regime. These authors' escape criterion reads Θ/Θ max t ≥ cot ψ/cot ψ Y , where Θ max t = 0.12 is the viscous fluid entrainment threshold (see section 2.1), ψ Y = 30 • the pocket angle for particles resting within the deepest pockets of the bed surface, and sin ψ = sin ψ Y + v r 2 /(2gd). This criterion means that the rebound kinetic energy only needs to uplift a particle rebounding within the deepest pocket to a point at which the near-surface flow is able to push it out of the pocket. The threshold resulting from this modified constraint is denoted as Θ Rb * t . 2. Models of the impact entrainment threshold Θ ImE t [Claudin and Andreotti, 2006;Kok, 2010a] do not neglect captures of continuous rebounders and therefore take into account the entrainment of bed particles. One possible way to do this is by combining rebound boundary conditions with an additional constraint that describes that one particle leaves the surface per impact on average (e.g., |v i | ∝ √g d [Claudin and Andreotti, 2006]). However, the incorporation of entrained particles as part of the average trajectory leads to consistency problems (see third point in the list above). 3. Hybrids between continuous rebound and impact entrainment models [Berzi et al., 2016[Berzi et al., , 2017 look for the smallest Shields number (denoted as Θ Rb |ImE t ) that results in a periodic trajectory under the constraint z h ≥ d (like before) and the additional constraint that particle-bed impacts do not lead to entrainment. Berzi et al. [2016Berzi et al. [ , 2017 modeled the latter constraint via |v i |/ √g d ≤ ζ/2 ≈ 20 (cf. equation (17)), which assumes that the fastest particles represented by the average trajectory of continuous rebounders do not exceed the value ζ of the nondimensionalized impact velocity that is associated with the onset of entrainment (which can be roughly justified by assuming an even impact velocity distribution between 0 and ζ). However, Pähtz and Durán [2018a] pointed out that this additional constraint is inconsistent with the experimental and numerical evidence that impact entrainment to be effective requires that the transport rate is significantly larger than zero (see section 4.1.3), which is never the case at the rebound threshold Θ Rb t in the absence of entrainment by turbulent fluctuation events (see section 4.2.2 and equation (29)). Consistently, Pähtz et al. [2020] showed that, in the limit Θ → Θ Rb t , identical periodic trajectories of continuous rebounders are unstable against trajectory fluctuations. That is, the energy that a particle must acquire upon entrainment to become a continuous rebounder is equal to the rebound energy of the continuous rebounder that has entrained it in this limit. This requirement contradicts the fact that the entrainment energy is much smaller than the rebound energy because of energy conservation, which implies that impact entrainment is impossible in this limit (see also discussion in section 4.2.1).
Apart from these conceptual differences, existing IPTMs differ in several details (partly summarized in Table 1): the form of the fluid drag law, the consideration or neglect of vertical drag forces on the particle motion, the form of the mean flow velocity profile (including the question of whether the viscous sublayer of the turbulent boundary layer is considered; for more details, see Appendix), and the bed boundary conditions (including the incorporation of viscous damping in the rebound laws). In this regard, it is reiterated that the effects of vis- cous damping on the dynamics of particle-bed rebounds are probably negligible for bedload transport (for which viscous damping is deemed as potentially significant), even for conditions with strongly damped binary particle collisions (see section 4.1.1).
In order to facilitate a comparison between the different model types that does not depend on modeling details but focuses only on conceptual differences, the same mean flow velocity profile (equation (A1), which includes the viscous sublayer), boundary conditions (equations (20) and (21)), and fluid drag law (the drag law by Camenen [2007]) are used for all model types. Following the trajectory calculation by Pähtz et al. [2020], the impact velocity v i as a function of the rebound velocity v r approximates aŝ v iz =v r z −t h , witht h = 1 +v r z + W − (1 +v rz ) e −(1+v r z ) , where t h is the hop time, W the principal branch of the Lambert-W function, V s ≡ v s / √ sgd the dimensionless value of the settling velocity v s (defined in equation (31)), Z ∆ d = 0.7d the average elevation of the particles' center during particle-bed rebounds (obtained from experiments [Dey et al., 2012;Hong et al., 2015]), and √ Θ f expresses the nonfluctuating wall-bounded flow after Guo and Julien [2007], with f the function given in equation (A1). Furthermore, the hat denotes nondimensionalized quantities using combinations ofg and v s , which is given by where C ∞ d = 1 and m = 1.5 are parameter values associated with the drag law for naturally shaped particles. Equations (20), (21), (30a), and (30b) can be iteratively solved for Θ(Ga, s,v r z ). Then the thresholds are obtained from where the hop height is given by In equations (32a)-(32d), the rebound threshold Θ Rb * t is the only modeled cessation threshold that is linked to the viscous yield stress Θ max t and thus to dense granular flow rheology (see section 2.1). In a complete model covering all transport regimes, such a connection must exist because Θ max t represents an upper limit to any kind of cohesionless sediment transport threshold. Also, a complete model of any kind of cohesionless transport threshold must reach this maximum value in the limit of vanishing particle inertia (i.e., when typical particle velocities during a trajectory become much smaller than √g d). The characteristic particle velocity scale in IPTMs is given by the settling velocity v s , which scales as v s ∝ Ga √ sgd in the viscous regime (Eq. (31) for small Ga). That is, a complete model of any kind of cohesionless transport threshold must approach Θ max t in the limit v s / √g d ∝ Ga √ s → 0, where Ga √ s can be interpreted as a Stokes-like number [Berzi et al., 2016[Berzi et al., , 2017Clark et al., 2017;Pähtz and Durán, 2018a].  (2.65, 40, 190, 2200, and 250000) corresponding to five different fluvial or aeolian conditions (Water, Venus, Titan, Earth, and Mars). These figures also show cessation threshold measurements obtained for nearly cohesionless conditions using different experimental methods. For turbulent bedload transport driven by water, the compilation of reference method-based measurements (measurement mean and its 95% confidence interval) by Buffington and Montgomery [1997], which make up a large portion of the Shields diagram, is shown. As explained in section 4.2.3, this method yields approximately the rebound threshold Θ Rb t . For viscous bedload transport driven by water-oil mixtures, the visual incipient motion measurements by Yalin and Karahan [1979] and Loiseleux et al. [2005] and cessation threshold measurements by Ouriemi et al. [2007] are shown (for viscous bedload transport, the differences between transport initiation, rebound, and impact entrainment threshold are very small [Pähtz and Durán, 2018a]). For aeolian saltation transport, a few studies [e.g., Ho, 2012;Zhu et al., 2019] carried out an indirect extrapolation to vanishing transport to obtain Θ Rb t using a proxy of Q: the surface roughness z o (see Appendix for its definition in the absence of transport), which undergoes a regime shift when saltation transport ceases. Furthermore, visual measurements of Θ Rb t by Bagnold [1937] and Chepil [1945] are shown, obtained from successively decrementing Θ until intermittent saltation transport stops. Direct measurements of the intermittent saltation transport threshold (and thus Θ Rb t ), based on the so-called Time Frequency Equivalence  (2.65, 40, 190, 2200, and 250000) corresponding to five different fluvial or aeolian conditions (Water, Venus, Titan, Earth, and Mars). Symbols correspond to threshold measurements (or measurement compilations) from various studies [Bagnold, 1937;Chepil, 1945;Buffington and Montgomery, 1997;Loiseleux et al., 2005;Ouriemi et al., 2007;Ho, 2012;Martin and Kok, 2018;Zhu et al., 2019] and methods (see text). Ouriemi et al. [2007] did not report single measurement values but a constant threshold 0.12 ± 0.03 for a large range of viscous conditions, indicated by the dotted square. Error bars correspond to 95%-confidence intervals of the compilation of reference method-based measurements by Buffington and Montgomery [1997], which make up a large portion of the Shields diagram.
Method (TFEM) [Wiggs et al., 2004], by Martin and Kok [2018] are also shown. Note that, although the evidence that the thresholds obtained from extrapolation to vanishing transport and from direct measurements of the cessation of intermittent saltation transport correspond to the rebound threshold Θ Rb t is quite strong (see section 4.1.3 and 4.2.2), many aeolian researchers believe that they correspond to the impact entrainment threshold Θ ImE t [e.g., Martin and Kok, 2018]. One of the reasons for this belief can be seen in Figure 17b: the prediction of Θ ImE t from equation (32c) is consistent with aeolian saltation transport data on Earth despite not containing fit parameters. In fact, for the range of conditions corresponding to these data, the predictions of Θ Rb * t and Θ Rb * * t by equations (32a) and (32b) are equivalent and, coincidentally, very close to the predictions of Θ ImE t and Θ Rb |ImE t by equations (32c) and (32d), which are also equivalent to each other. At this point, it is worth reiterating that differences between the models caused by differences in the modeling details (e.g., those in Table 1) have been excluded here. Such detail differences cause the predictions of existing models to differ more strongly from one another than shown here. Figures 17a, 17b, and 18a show that the predictions of Θ Rb * * t , Θ ImE t , and Θ Rb |ImE t from equations (32b)-(32d) overestimate threshold measurements for fluvial bedload transport by at least an order of magnitude. For Θ Rb * * t and Θ Rb|ImE t , this overestimation is caused by the constraint in the minimization of Θ that the particle hop height z h must exceed one particle diameter d to escape the bed surface (equations (32b) and (32d)), preventing solutions with small particle velocities that would have a smaller threshold. However, the prediction of Θ Rb * Venus, Titan, Earth, and Mars). Symbols correspond to threshold measurements (or measurement compilations) from various studies [Bagnold, 1937;Chepil, 1945;Buffington and Montgomery, 1997;Loiseleux et al., 2005;Ouriemi et al., 2007;Ho, 2012;Martin and Kok, 2018;Zhu et al., 2019] and methods (see text). aeolian and fluvial transport regimes strongly supports modeling nonsuspended sediment transport within the continuous rebound framework.

Models Based on Correlations Between Transport Layer-Averaged Physical Quantities
Existing correlation-based cessation threshold models start with the assumption of a constant bed friction coefficient µ b [Pähtz et al., 2012;Pähtz and Durán, 2018a] (µ b is the inverse of the parameter α in the model by Pähtz et al. [2012]). As discussed in section 4.2.2, the approximate constancy of µ b has been analytically linked to continuous rebounds [Pähtz and Durán, 2018b]. However, in contrast to the purely kinematic meaning of µ b in IPTMs (equation (31)), for realistic nonsuspended sediment transport, µ b conveys information about both the particle kinematics and interparticle contacts. Note that µ b const is also predicted by IPTMs when vertical drag forces are small (i.e., the buoyancy-reduced gravity force dominates the vertical motion) because this fixes e z 1 and thus e 2D , θ 2D r , and µ b via the rebound laws [Pähtz et al., 2020].
A constant µ b links the average horizontal fluid drag acceleration a dx to the buoyancyreduced gravityg via µ b a dx /g, where the overbar denotes a particle concentration-weighted height average [Pähtz and Durán, 2018a] (which is equal to the average over the fluvial environments. The predictions from these two models differ for turbulent bedload transport (large U x / Θ Rb * t and small transport layer: Z ∼ Z ∆ ), mainly because the scaling V x 2 Θ Rb * t /κ that equation (33e) predicts for large U x / Θ Rb * t does not capture the outcome of the constrained minimization in equation (32a) for Z ∼ Z ∆ . The predictions from these two models also differ for viscous saltation transport because the model by Pähtz and Durán [2018a] neglects vertical drag forces at various instances. However, note that this model does not completely neglect vertical drag forces because the scaling V x ∝ U x that equation (33e) predicts for viscous saltation transport is associated with vertical drag [Pähtz and Durán, 2018a], which is why deviations between this model and the analogous IPTM are only moderate in this regime.

Open Problem: Reliable Models of the Impact Entrainment Threshold and Planetary Saltation Transport
Existing models of the impact entrainment threshold [Claudin and Andreotti, 2006;Kok, 2010a;Pähtz et al., 2012], which is arguably also the continuous transport threshold, do not take into account that the transport rate Q is significantly larger than zero at Θ ImE t , even in the absence of entrainment by turbulent fluctuation events (see section 4.1.3). Instead, Q vanishes at the rebound threshold Θ Rb t , which is smaller than Θ ImE t (see section 4.2.2). Likewise, as mentioned before, existing models of Θ ImE t effectively assume that all entrained particles exhibit a kinetic energy that allows them to participate in the continuous rebound motion even though most of them do not [Pähtz and Durán, 2018a]. For these reasons, existing impact entrainment threshold models seem to be missing important physics and need to be improved. This is problematic for modeling and predicting extraterrestrial sediment transport and associated bedform evolution [e.g., Almeida et al., 2008;Bourke et al., 2010;Kok, 2010b;Ayoub et al., 2014;Lorenz, 2014;Rasmussen et al., 2015;Jia et al., 2017;Telfer et al., 2018;Durán Vinent et al., 2019] because most predictions of the aeolian saltation transport rate require that transport is continuous (i.e., at capacity).

Main Difference Between Aeolian and Fluvial Rebound Threshold
The most important difference between aeolian saltation and fluvial bedload transport is the largely different density ratio s, which ranges from close to unity for oil and water to the order of 10 5 for air on Mars. Equations (33a)-(33e) elucidate that s affects the modeled rebound threshold Θ Rb * t in a relatively simple manner. In fact, it can be seen that s explicitly appears only in Eq. (33d), which describes a monotonous increase of the dimensionless transport layer thickness Z with s. Subsequently, Z monotonously increases the dimensionless transport layer-averaged flow velocity U x via equation (33b). That is, given a certain solution Θ Rb * t (Ga, s) of equations (33a)-(33e), an increase of s leads to an increase of U x , which must be compensated by a decrease of U x via a decrease of Θ Rb * t to achieve a new steady solution for the same Galileo number Ga. This mathematical fact expresses the physical fact that particles that stay longer in the flow can feel a given effective flow forcing at a lower fluid shear stress, which is the ultimate reason for why Θ Rb * t decreases with s for a given Ga.

Summary and Outlook
Section 1 outlined five old, yet very significant, inconsistencies related to the concept of a threshold shear stress for incipient motion. For the concept of a threshold shear stress to be physically meaningful, these inconsistencies must be addressed and resolved. They can be briefly summarized as follows: 1. By design, existing models of incipient motion capture the conditions to which they have been adjusted (aeolian or fluvial transport on Earth). However, the predictions from standard models adjusted to aeolian transport on Earth [Iversen and White, 1982;Shao and Lu, 2000] are in stark disagreement with recent observations of aeolian transport on Mars [e.g., Sullivan and Kok, 2017;Baker et al., 2018]. 2. Because of turbulent fluctuation events, fluid entrainment gives rise to fluvial bedload transport even for Shields numbers much below the Shields curve [Paintal, 1971], which is the curve that is thought to describe incipient motion. 3. Below the respective shear stress threshold associated with incipient motion, turbulent fluctuation events are able to initiate fluvial bedload transport but not aeolian saltation transport. However, there is no reason to believe that the physics of incipient motion are different in aeolian and fluvial environments. 4. Old experiments indicate a nonnegligible role of particle inertia in fluvial bedload transport [Ward, 1969;Graf and Pazis, 1977], which is problematic because critical conditions are defined via nonzero transport rates (i.e., particles are in motion at threshold conditions). 5. In old numerical simulations of turbulent fluvial bedload transport [no and García, 1998], it was recognized that the threshold shear stress obtained from extrapolating the simulated capacity transport rate to vanishing transport may not be associated with fluid entrainment. This is problematic because many of the threshold data compiled in the Shields diagram have been obtained from such or similar extrapolation methods [e.g., Shields, 1936].
As a result of the latest research reviewed in sections 2-4, a new conceptual picture has emerged (section 5.1) that resolves these problems. However, it must be emphasized that this conceptual picture represents the authors' synthesis of the current state of the art, and many aeolian and fluvial geomorphologists may disagree. This is because, in some places, it stands in stark contrast to what has been a century old consensus. Likewise, there are still many open problems and controversies, summarized in section 5.2 (and highlighted in sections 2-4), as well as a number of issues that have not been discussed in this review (e.g., the effects of particle size heterogeneity on transport thresholds and sediment entrainment), into which section 5.3 presents a brief outlook. Figure 19 summarizes the various shear stress thresholds of nonsuspended sediment transport and their relations to and effects on the transport characteristics. Details, with references to the research reviewed in sections 2-4, are described below.

Creeping (Θ > 0)
Creeping (see section 2.3) refers to a superslow granular motion, usually in the form of intermittent local particle rearrangements within the sediment bed (not limited to the bed surface), that occurs below a macroscopic yield criterion (see section 2.1). One form of creeping is triggered by nearby regions above yield, while another form (the origin of which is not fully understood) occurs even in the absence of such regions. The existence of the latter form implies that sediment likely is always transported (albeit slowly) for arbitrarily small values of the Shields number Θ, even in the absence of turbulence. Creeping of both kinds is very important in determining the particle motion near transport initiation. It is fundamentally related to the granular material, not a purely fluid-driven effect.

Viscous Yield Stress Θ max
t Apart from creeping, which affects the entire granular bed, bed surface particles can be entrained directly by flow forces. When a sediment bed is subjected to a laminar flow at a sufficiently low shear Reynolds number Re * , there is a critical Shields number, the yield   stress Θ max t , above which motion of bed surface particles is initiated and then never stops, whereas potential transient motion below Θ max t will inevitably come to an end (see section 2.1). The viscous yield stress Θ max t constitutes the upper limit for any kind of cohesionless sediment transport threshold, including the Shields curve. The values of Θ max t reported in the literature are somewhat scattered (between 0.1 and 0.4), but these numbers are within the range of the bulk friction coefficients for granular materials (ranging from low friction spheres to more frictional, rough particles), suggesting that the granular material's yield condition (see section 2.1) is very important in determining the viscous yield stress.

Initiation Threshold Θ In t
While for laminar flows, the entrainment of individual bed surface particles is controlled by a critical Shields number, the entrainment of individual bed surface particles by turbulent flows is better described by an impulse (section 3.2.1) or energy (section 3.2.2) criterion. Nonetheless, one can still define an initiation threshold Shields number Θ In t (see section 3.3) at which the probability of fluid entrainment exceeds zero (i.e., Θ In t Θ max t for laminar flows at sufficiently low Re * ). Because fluid entrainment is predominantly caused by turbulent fluctuation events, Θ In t depends not only on Re * but also on properties that control the size of the largest turbulent flow eddies, such as the turbulent boundary layer thickness. This may be one of the reasons why aeolian incipient motion models adjusted to wind tunnel measurements fail when applied to atmospheric boundary layers (see first problem outlined at the beginning of section 5). Further possible reasons include atmospheric instability, topography gradients, surface inhomogeneities, such as obstacles and vegetation, and sublimation of subsurface ice in natural atmospheres (see section 3.3.3).

Rebound Threshold Θ Rb t (Generalized Shields Curve)
The rebound threshold Θ Rb t (see section 4.2) is largely unrelated to the entrainment of bed sediment (except for viscous bedload transport, for which Θ Rb t Θ max t ) but describes the minimal dimensionless fluid shear stress that is needed for the mean turbulent flow to compensate the average energy loss of rebounding particles by fluid drag acceleration during their trajectories. Hence, for Θ ≥ Θ Rb t , transported particles rebound for comparably longer periods before they deposit, whereas they deposit very quickly for Θ < Θ Rb the field, the magnitude of Θ In t relative to Θ Rb t is unclear as Θ In t is smaller than in wind tunnels because of a much larger boundary layer thickness δ, since δ controls the size of the largest turbulent eddies and thus entrainment by turbulent fluctuation events (see section 3.3).

Implications for Field Phenomenona
The new conceptual picture described above has been derived nearly entirely from theoretical and laboratory investigations. One may therefore wonder to what degree does the notion of various transport thresholds have implications for natural field conditions, such as bedload transport in rivers and saltation transport driven by planetary winds. There are three major aspects in which the field differs from most laboratory experiments: much broader particle size distributions, much larger and more unstable boundary layers (mainly for aeolian transport), and various kinds of surface inhomogeneities, such as bedforms, obstacles, and vegetation. The effects of particle size heterogeneity have been excluded from this review (they are briefly discussed in the outlook, section 5.3.1). The remaining two aspects are both associated with increasing turbulence and thus fluid entrainment (see section 3). In contrast, the rebound threshold Θ Rb 4. Is aeolian transport in the field on Earth and other planetary bodies, in contrast to wind tunnels with similar atmospheric pressure conditions, always being initiated close to the rebound threshold Θ Rb t because of thick boundary layers, atmospheric instability, topography inhomogeneities, and subsurface ice sublimation? The answer to this question is probably the most important one, since a positive answer would imply that a reliable model for Θ In t (i.e., answers to the previous three questions) is not required for predicting aeolian processes on such bodies. 5. Does equilibrium aeolian bedload transport (i.e., h ∼ d) exist in the field because of thick boundary layers? 6. Direct measurements of aeolian sediment transport initiation, which are currently missing, can help answering the questions above.
Section 4: 1. For a particle collision with a static sediment bed: how does the rebound probability P r depend on impact velocity and angle? 2. How do particle shape and size distribution affect particle-bed collisions? 3. Does viscous damping truly not much affect particle-bed collisions, as suggested by the insensitivity of DEM-based sediment transport simulations to the normal restitution coefficient of binary collisions? And if so, what is the physical reason? 4. How do cohesive interparticle forces affect the collision process and thus sediment transport cessation? 5. How do the laws describing a particle collision with a static bed change for a particle collision with a mobile bed? 6. It is straightforward to define intermittent and continuous sediment transport for the absence of fluid entrainment because the sediment transport rate exhibits a discontinuous jump from nearly zero to a finite value at the continuous transport threshold. However, how does one universally define intermittent and continuous transport if fluid entrainment does occur? 7. Is the transition from intermittent to continuous aeolian saltation transport associated with fluid entrainment (the current consensus) or with impact entrainment (the authors' opinion, based on recent developments in the field)? 8. What controls the impact entrainment threshold Θ ImE t and how does one model it?

Outlook
To limit the scope of this review, several important topics have been excluded. Two of them are briefly discussed below.

Effects of Particle Size Heterogeneity on Sediment Transport Initiation, Cessation, and Entrainment
Perhaps the most important topic that has been excluded from this review is the effects of the heterogeneity of the size of bed surface particles on sediment transport initiation, cessation, and entrainment. Naturally, sediment transport initiation and entrainment are sizeselective. However, it is unclear whether this is also true for sediment transport cessation. While the continuous rebound mechanism (see section 4.2) is clearly a size-selective process (coarser particles are less accelerated during their trajectories), impact entrainment may not be [Martin and Kok, 2019;Zhu et al., 2019]. Furthermore, in heterogeneous sediment beds, relatively fine particles tend to be surrounded by coarser ones and their protrusion (i.e., the particle height above surrounding sediment) is thus smaller than on average, whereas relatively coarse particles tend to have a larger-than-average protrusion. Because driving forces decrease and resisting forces increase with decreasing protrusion [Yager et al., 2018], relatively fine particles are more difficult to be entrained when compared with a bed made only of such fine particles. The ability of fine particles to continuously rebound is also suppressed by the presence of coarse particles [Zhu et al., 2019]. All these effects can make heterogeneous sediment beds much less mobile than homogeneous ones of the same median particle size. For example, for both fluvial bedload [MacKenzie and Eaton, 2017;MacKenzie et al., 2018] and aeolian saltation transport [Zhu et al., 2019], it was found for certain heterogeneous beds that the largest particles of the particle size distribution (larger than the 90th percentile) have a very strong control on overall mobility. However, the manner and degree of the heterogeneousness seem to play an important role as not all kinds of heterogeneous beds are so strongly affected by the presence of large particles [e.g., Wilcock, 1993;Martin and Kok, 2019]. In particular, in the early stages of bed armoring, the sediment transport rate can increase because collisions between transported fine particles and coarse bed particles are more elastic than collisions between particles of the same size [Bagnold, 1973].

Effects of Steep Bed Slope on Sediment Transport Initiation, Cessation, and Entrainment
Another important topic that has been excluded from this review is the effects of steep bed slope angles on sediment transport initiation, cessation, and entrainment. For example, horizontal downslopes should, if everything else stays the same, increase bed mobility because of the additional horizontal gravity force acting on particles [Maurin et al., 2018]. However, in fluvial environments, steep slopes are usually accompanied by a very small water depth of the order of one particle diameter (or even lower), which strongly suppresses the magnitude of hydrodynamic forces acting on particles, thus decreasing rather than increasing bed mobility [Prancevic and Lamb, 2015]. Then again, an increasing downslope angle α increases the bulk friction coefficient µ within the sediment bed (for turbulent flows, [Maurin et al., 2018], where φ b is the bed volume fraction). Once µ exceeds the static friction coefficient µ s associated with the yielding transition (see section 2.1), the entire bed fails and a debris flow forms [Takahashi, 1978;Prancevic et al., 2014;Cheng et al., 2018]. Kudrolli [2017, 2018], Diplas et al. [2008], Valyrakis [2013], Williams et al. [1994], Iversen and White [1982], Burr et al. [2015], Swann et al. [2020], Beladjine et al. [2007], Vowinckel et al. [2016], Salevan et al. [2017], Buffington and Montgomery [1997], Yalin and Karahan [1979], Loiseleux et al. [2005], Ouriemi et al. [2007], Bagnold [1937], Chepil [1945], Ho [2012], Martin and Kok [2018], and Zhu et al. [2019]. We thank Michael Church and three anonymous reviewers for their critical reviews and numerous insightful comments and suggestions. T.P. acknowledges support from grant National Natural Science Foundation of China (No. 11750410687 Q Sediment transport rate [kg/(ms)] Θ ≡ τ/((ρ p − ρ f )gd) Shields number or Shields parameter s ≡ ρ p /ρ f Particle-fluid-density ratio Re ≡ U b H/ν f Reynolds number Re * ≡ u * d/ν f Shear Reynolds number Ga ≡ (s − 1)gd 3 /ν f Galileo number (also called Yalin parameter) St ≡ s|v r |d/(9ν f ) Stokes number, where |v r | is the relative velocity between two particles just before they collide M * ≡ M/(ρ p d) Nondimensionalized sediment transport load Q * ≡ Q/(ρ p d (ρ p /ρ f − 1)gd) Nondimensionalized sediment transport rate I ≡ γd/ P/ρ p Inertial number J ≡ ρ f ν f γ/P Viscous number K ≡ J + c K I 2 Viscoinertial number, where c K is an order-unity fit parameter Pe ≡ γd/

√
T Péclet number C m = 1/2 Added mass coefficient κ = 0.4 von Kármán constant ψ Pocket angle ψ Y Pocket angle for particles resting within the deepest pockets of the bed surface L arm Lever arm length [m] α Bed slope angle ∆Z Critical dimensionless vertical particle displacement required for entrainment ∆X Critical dimensionless horizontal particle displacement required for entrainment µ C Effective Coulomb friction coefficient encoding the combined effects of sliding and rolling friction in entrainment µ ≡ −τ p /P Ratio between particle shear stress and particle pressure (bulk friction coefficient) µ g Surface friction coefficient of granular particle µ s Static friction coefficient of granular bulk (yield stress ratio) µ b Bulk friction coefficient at the interface between bed and transport layer. In contrast to µ and µ s , µ b includes contributions from stresses associated with the particle fluctuation motion in addition to contributions from intergranular contacts. ξ ∝ |µ − µ s | −ν Correlation length associated with the yielding transition, where ν = 0.5 is the critical exponent F, F D , F L , F t , F n , F e Instantaneous force applied by the fluid on a particle [kgm/s 2 ]. Subscript (D, L, t, n, e) refers to nature of force (drag, lift, tangential, normal, effective). T , T D , T L , T t , T n , T e Duration of turbulent fluctuation event [s]. Subscript (D, L, t, n, e) refers to nature of applied fluid force (drag, lift, tangential, normal, effective). I f Impulse of turbulent fluctuation event [kgm/s] E f Energy of turbulent fluctuation event [kgm 2 /s 2 ] F c Force resisting initial particle motion [kgm/s 2 ] u c Critical instantaneous local flow velocity associated with resisting forces [m/s] I f c Critical impulse required for fluid entrainment [kgm/s] W c Critical work done by flow event required for fluid entrainment [kgm 2 /s 2 ] C eff Energy transfer coefficient, describing the fraction of energy transferred from flow to target particle during turbulent fluctuation event C ≡ α −1 f f (G) √ sd/δ Inverse dimensionless boundary layer thickness α f ≡ u m /u Ratio between the characteristic flow velocity u m associated with the largest turbulent fluctuations and the local mean flow velocity u T max Maximal duration of turbulent fluctuation events [s] f (G) Factor that encodes information about particle shape, orientation, and the pocket geometry v i Impact velocity [m/s] v r (v 2D r ) Rebound velocity (projected into incident plane) [m/s] v e (v 2D e ) Ejection velocity (projected into incident plane) [m/s] θ i Impact angle θ r (θ 2D r ) Rebound angle (projected into incident plane) θ e (θ 2D e ) Ejection angle (projected into incident plane) E i Impact energy [kgm 2 /s 2 ] E e (E 2D e ) Ejection energy (projected into incident plane) [kgm 2 /s 2 ] N e Average number of ejected particles P r Rebound probability (e 2D ≡ |v 2D r |/|v i |) e ≡ |v r |/|v i | (Projected) rebound restitution coefficient e z ≡ −v r z /v iz Vertical rebound restitution coefficient A, B, A 2D , B 2D , χ, r, r 2D , n 0 , ζ Dimensionless parameters appearing in empirical or semiempirical relations describing the collision process between an incident bead and a granular packing α r Normal rebound restitution coefficient in the impact plane β r Tangential rebound restitution coefficient in the impact plane Restitution coefficient for binary particle collision V b Effective value of the local particle velocity averaged over elevations near the bed surface [m/s] f in Particle feeding frequency at flume entrance [1/s] f out Frequency of particles passing an illuminated window near the flume exit [1/s] f Q Fraction of active aeolian saltation transport v ↑ Initial particle velocity in thought experiment in section 4.2.1 [m/s] E ↑ Initial particle energy in thought experiment in section 4.2.1 [kgm 2 /s 2 ] E c Critical energy that E ↑ must exceed for particle to continuously rebound along the surface [kgm 2 /s 2 ] n tr /n tot Number of transported particles relative to the total number of bed surface particles n v t Number of particles that are faster than a certain velocity threshold v t z h Hop height [m] t h Hop time [s] f (Re * , z/d) Function given by equation (A1) v s Settling velocity [m/s] U x ≡ u x / √ sgd Dimensionless transport layer-averaged fluid velocity V x ≡ v x / √ sgd Dimensionless transport layer-averaged horizontal particle velocity sgd Dimensionless transport layer-averaged vertical particle velocity Z ≡ z/d Dimensionless transport layer thickness Z ∆ = 0.7 Dimensionless average elevation of the particles' center during particle-bed rebounds z o Surface roughness [m] c 1 , c 2 , c 3 Model constants in equations (33c)-(33e) Bed sediment entrainment Mobilization of bed sediment Fluid entrainment Entrainment caused by the action of flow forces Incipient motion Initiation of sediment transport by fluid entrainment Impact entrainment Entrainment caused by the impacts of transported particles onto the bed Sediment transport Sediment motion caused by the shearing of an erodible sediment bed by flow of a Newtonian fluid Aeolian sediment transport Wind-driven sediment transport Fluvial sediment transport Liquid-driven sediment transport (despite its name, not limited to fluvial environments) Nonsuspended sediment transport Sediment transport in which the fluid turbulence is unable to support the submerged particle weight Saltation transport Nonsuspended sediment transport with comparably large transport layers (h d) Bedload transport Nonsuspended sediment transport with comparably small transport layers (h ∼ d) Transport capacity (or saturation) Loosely, the maximum amount of sediment a given flow can carry without causing net sediment deposition at the bed. More precisely, in the context of nonsuspended transport of nearly monodisperse sediment, it is defined as a steady transport state at which any further net entrainment of bed sediment into the transport layer would weaken the mean turbulent flow to a degree at which it is no longer able to compensate the average energy loss of particles rebounding with the bed by their energy gain during their trajectories via fluid drag acceleration. The so defined transport capacity obeys equation (29). Creeping A superslow granular motion, usually in the form of intermittent local particle rearrangements within the sediment bed (not limited to the bed surface), that occurs below a macroscopic yield criterion Θ t (τ t ) Shields number (fluid shear stress) at a nonspecified transport threshold. For specifications, see below. t conv ∝ |Θ − Θ t | −β Time scale for transport property to converge in the steady state near Θ t , where β is a positive exponent Shields diagram (Shields curve) Diagram compiling measurements of Θ t as a function of Re * (the Shields curve Θ t (Re * )) for fluvial bedload transport conditions Θ max t (τ max t ) Viscous yield stress. The upper limit of the threshold Shields number (fluid shear stress) in the Shields diagram, which is associated with viscous bedload transport. For Θ Θ max t , a sediment bed subjected to a laminar flow at low Shear Reynolds number Re * may temporarily fail but will eventually rearrange itself into a more stable packing that resists the applied fluid shear stress. For Θ Θ max t , a sediment bed subjected to a laminar flow can no longer find packing geometries that are able to resist the applied fluid shear stress. Θ In t (τ In t ) Initiation threshold. Shields number (fluid shear stress) at which the probability of fluid entrainment of bed particles exceeds zero (which, for turbulent fluvial bedload transport, occurs much below the Shields curve). For sediment beds subjected to turbulent flows, a critical fluid shear stress does no longer describe the fluid entrainment of individual particles. However, one can still define a Shields number (Θ In t ) below which fluid entrainment does never occur. Like for Θ max t , transient behavior associated with the flow temporarily pushing particles from less stable to more stable pockets is excluded in the definition of Θ In t . Θ In t (τ In t ) Rocking initiation threshold. Shields number (fluid shear stress) above (below) which there is a nonzero (zero) probability that peaks of flow forces associated with turbulent fluctuation events acting on bed particles exceed resisting forces. That is, there is a nonzero probability that particles rock (or wobble or oscillate) within their bed pockets. Rocking may (Θ In t = Θ In t ) or may not (Θ In t < Θ In t ) lead to complete entrainment depending on the maximal duration of the strongest possible turbulent fluctuation events. Θ Rb t (τ Rb t ) Rebound threshold. Shields number (fluid shear stress) above which the mean turbulent flow is able to compensate the average energy loss of transported particles rebounding with the bed by their energy gain during their trajectories via fluid drag acceleration, giving rise to a long-lasting rebound motion. In general, this threshold is unrelated to the entrainment of bed sediment. It is also the threshold that appears in most threshold shear stress-based sediment transport expressions. Θ Rb * t (Θ Rb * * t ) Modeled rebound threshold. Values of Θ Rb t from models that consider (neglect) that the near-surface flow can assist rebounding particles in escaping the bed surface. Θ ImE t (τ ImE t ) Impact entrainment threshold. Shields number (fluid shear stress) above which entrainment of bed sediment by impacts of transported particles onto the bed is able to compensate captures of long-lasting rebounders (see Θ Rb t above) by the bed. This threshold is arguably also the threshold of continuous nonsuspended sediment transport. Θ Rb |ImE t (τ Rb |ImE t ) Modeled hybrid between rebound and impact entrainment threshold