Scaling laws for planetary sediment transport from DEM-RANS numerical simulations

Abstract

We use an established discrete element method (DEM) Reynolds-averaged Navier–Stokes (RANS)-based numerical model to simulate non-suspended sediment transport across conditions encompassing almost seven orders of magnitude in the particle–fluid density ratio $s$, ranging from subaqueous transport ($s=2.65$) to aeolian transport in the highly rarefied atmosphere of Pluto ($s=10^7$), whereas previous DEM-based sediment transport studies did not exceed terrestrial aeolian conditions ($s\approx 2000$). Guided by these simulations and by experiments, we semi-empirically derive simple scaling laws for the cessation threshold and rate of equilibrium aeolian transport, both exhibiting a rather unusual $s^{1/3}$-dependence. They constitute a simple means to make predictions of aeolian processes across a large range of planetary conditions. The derivation consists of a first-principle-based proof of the statement that, under relatively mild assumptions, the cessation threshold physics is controlled by only one dimensionless control parameter, rather than two expected from dimensional analysis. Crucially, unlike existing models, this proof does not resort to coarse-graining the particle phase of the aeolian transport layer above the bed surface. From the pool of existing models, only that by Pähtz et al. (J. Geophys. Res.: Earth, vol. 126, 2021, e2020JF005859) is somewhat consistent with the combined numerical and experimental data. It captures the scaling of the cessation threshold and the $s^{1/3}$-dependence of the transport rate, but fails to capture the latter's superimposed grain size dependence. This hints at a lack of understanding of the transport rate physics and calls for future studies on this issue.

1. Introduction

Aeolian (wind-driven) transport of non-suspended grains, including sand, ice and snow, is a ubiquitous phenomenon that leads to a rich variety of multiscale bedforms on Earth and other planetary bodies (Bourke et al. 2010; Kok et al. 2012; Diniega et al. 2017). As suggested by the presence of wind streaks and dunes, it may even occur in the very rarefied atmospheres of Neptune's moon Triton (Sagan & Chyba 1990), Pluto (Telfer et al. 2018) and the comet 67P/Churyumov-Gerasimenko (Thomas et al. 2015; Jia, Andreotti & Claudin 2017).

Driven by fluid drag and gravity, most transported sand-sized and larger grains regularly interact with the bed surface as flow turbulence is too weak to suspend them. For denser fluids, such as water and most other liquids, this near-surface grain motion occurs in the form of rolling, sliding and small hops (bedload), whereas for lighter fluids, like most gases, grains move in more energetic hops (saltation). At equilibrium, the deposition of transported grains on the bed is exactly balanced by the entrainment of bed grains into the transport layer. The rate at which equilibrium aeolian transport takes place and the threshold wind speed below which it ceases constitute the two arguably most important statistical transport properties in the context of bedform formation and evolution in natural environments (Kok 2010a; Durán Vinent et al. 2019). In particular, in natural environments, topography inhomogeneities, strong turbulent fluctuations and a variety of wind-unrelated mechanisms to generate airborne grains, along with very long natural sediment fetches, can plausibly initiate transport and lead to equilibrium transport above the cessation threshold (Pähtz et al. 2020, § 3.3.3.4). This may even be true in environments where the aeolian transport initiation threshold for an idealised flat sediment bed is much larger than the cessation threshold, like potentially on Mars (Kok 2010a), Pluto (Telfer et al. 2018) and Saturn's moon Titan (Comola et al. 2022), as well as in Antarctica. In fact, although Antarctica's surface is covered by very cohesive (Pomeroy & Gray 1990) old snow and ice (cohesion increases the initiation threshold probably much more than the cessation threshold (Comola et al. 2019b, 2022; Pähtz et al. 2021; Besnard et al. 2022)), aeolian snow and ice transport occurs there even at relatively low wind speeds that are likely much below the initiation threshold (Leonard et al. 2011).

Since the highly random, collective motion of bed and transported grains eludes a rigorous analytical description, existing physical models of equilibrium aeolian transport have relied on drastically coarse-graining the particle phase of the aeolian transport layer above the bed surface (Ungar & Haff 1987; Andreotti 2004; Claudin & Andreotti 2006; Kok & Renno 2009; Kok 2010b; Durán, Claudin & Andreotti 2011; Berzi, Jenkins & Valance 2016; Berzi, Valance & Jenkins 2017; Lämmel & Kroy 2017; Pähtz & Durán 2018a, 2020; Andreotti et al. 2021; Pähtz et al. 2021; Comola et al. 2022; Gunn & Jerolmack 2022). The most common modelling approach is to represent the grain motion by a single or multiple saltation trajectories. Depending on the number and kind of considered trajectories and the assumed outcome of grain–bed collisions, such models can yield fundamentally different scaling laws for the cessation threshold and/or equilibrium transport rate, with predictions varying by about an order of magnitude when applied to Martian-pressure atmospheric conditions (Pähtz et al. 2020; Gunn & Jerolmack 2022).

One reason for the strong variability of both existing cessation threshold and equilibrium transport rate predictions is a lack of consensus on the physical picture behind the cessation threshold. On the one hand, it has been modelled as an ‘impact entrainment threshold’ (Pähtz et al. 2020), the smallest wind velocity at which random captures of saltating grains by the bed can be compensated by the splash of bed grains due to grain–bed impacts (Andreotti 2004; Claudin & Andreotti 2006; Kok & Renno 2009; Kok 2010b; Andreotti et al. 2021; Comola et al. 2022). On the other hand, it has been modelled as a ‘rebound threshold’ (Pähtz et al. 2020), the smallest wind velocity required to replenish the energy saltating grains lose when rebounding with the bed, independent of grain capture and splash (Berzi et al. 2017; Pähtz et al. 2021; Gunn & Jerolmack 2022). We previously proposed and supported the hypothesis that both these dynamic thresholds play a role in saltation dynamics: the former as the dynamic threshold of continuous and the latter as the dynamic threshold of intermittent saltation and therefore as the actual cessation threshold (Pähtz & Durán 2018a; Pähtz et al. 2020; Pähtz et al. 2021). If true, this could have the unintended consequence that measurements of one are mistaken for the other dynamic threshold. For example, Pähtz et al. (2021) proposed that the recent dynamic-threshold measurements in a low-pressure wind tunnel by Andreotti et al. (2021) may constitute data of the continuous-transport threshold, and not of the cessation threshold as the experimenters claimed. This would be problematic as these data have been used to develop new cessation threshold models and compare their predictive capabilities with those of older ones (Andreotti et al. 2021; Gunn & Jerolmack 2022).

Here, we show that, under relatively mild assumptions, one can obtain insights into the physics of the cessation threshold without resorting to coarse-graining the particle phase of the aeolian transport layer above the bed surface. In detail, if the bed surface can be considered as a flat boundary, with scale-free boundary conditions describing the outcome of grain–bed collisions, and the driving wind as a smooth inner turbulent boundary layer flow that interacts with grains via Stokes drag, then the threshold shear velocity, appropriately non-dimensionalised, is a function of only one dimensionless control parameter, rather than two expected from dimensional analysis (§ 3). We confirm this prediction, and therefore its underlying assumptions, with numerical simulations using an existing discrete element method (DEM)-based numerical model (Durán, Andreotti & Claudin 2012, introduced in § 2) of equilibrium transport of cohesionless non-suspended sediments. The simulated transport conditions encompass almost seven orders of magnitude in the particle–fluid density ratio $s$, ranging from subaqueous transport ($s=2.65$) to aeolian transport in the highly rarefied atmosphere of Pluto ($s=10^7$), whereas previous DEM-based sediment transport studies did not exceed terrestrial aeolian conditions ($s\approx 2000$). We also use the simulation data to semi-empirically derive simple scaling laws for the cessation threshold and equilibrium transport rate, and to test existing models (§ 3). The derived scaling laws are consistent with experimental data, except the dynamic-threshold measurements by Andreotti et al. (2021), in line with the aforementioned hypothesis that the latter constitute data of the continuous-transport threshold rather than the cessation threshold (discussed in more detail in § 4).

2. Numerical model

We use the numerical model of Durán et al. (2012), which couples a continuum Reynolds-averaged description of hydrodynamics with a DEM for the grain motion under gravity, buoyancy and fluid drag. The drag force is given by $\boldsymbol {F_d}=\frac {1}{8}\rho _f{\rm \pi} d^2C_d|\boldsymbol {u_r}|\boldsymbol {u_r}$, where $\rho _f$ is the fluid density, $d$ the median grain diameter, $\boldsymbol {u_r}$ the fluid–grain velocity difference and

(2.1)\begin{equation} C_d=\left(\sqrt{\frac{Re_c}{|\boldsymbol{u_r}|d/\nu}}+\sqrt{C_d^\infty}\right)^2 \end{equation}

the drag coefficient, with $\nu$ the kinematic viscosity. Most simulations are carried out using the parameter values $Re_c=24$ and $C_d^\infty =0.5$, close to those for spherical grains (Camenen 2007), while a few simulations are carried out using different values (specified when done so) to test the effect of drag modifications, which may for example occur in very-low-pressure atmospheres due to drag rarefaction (Crowe et al. 2012). Spherical grains ($10^4\unicode{x2013}10^5$) with mild polydispersity are confined in a quasi-two-dimensional domain of length $\approx 10^3d$, with periodic boundary conditions in the flow direction, and interact via normal repulsion (restitution coefficient $e=0.9$) and tangential friction (contact friction coefficient $\mu _c=0.5$). The bottom-most grain layer is glued on a bottom wall, while the top of the simulation domain is reflective but so high that it is never reached by transported grains. The Reynolds-averaged Navier–Stokes (RANS) equations are combined with a semi-empirical mixing length closure that accounts for the viscous sublayer of the turbulent boundary layer and ensures a smooth hydrodynamic transition from high to low particle concentration at the bed surface:

(2.2)\begin{equation} \frac{\mathrm{d} l_m}{\mathrm{d} z}=\kappa\left[1-\exp\left(-\sqrt{\frac{u_xl_m}{7\nu}}\right)\right], \end{equation}

where $l_m(z)$ is the height-dependent mixing length, $\kappa =0.4$ the von Kármán constant and $u_x(z)$ the mean flow velocity field. This parametrisation quantitatively reproduces measurements of $u_x(z)$ in the absence of transport. Simulations with this numerical model are insensitive to $e$ and, therefore, insensitive to viscous damping (Pähtz & Durán 2018a,b). The simulations reproduce measurements of the rate and cessation threshold of terrestrial aeolian transport, and viscous and turbulent subaqueous transport (figures 1 and 3 of Pähtz & Durán (2018a) and figure 4 of Pähtz & Durán (2020)), height profiles of relevant equilibrium transport properties (figure 2 of Pähtz & Durán (2018a) and figure 6 of Durán, Andreotti & Claudin (2014a)) and aeolian ripple formation (Durán, Claudin & Andreotti 2014b).

2.1. Average of simulated quantities

We define two types of averages of a particle property $A_p$. Based on the spatial homogeneity of the simulations, the mass-weighted average of $A_p$ over all particles within an infinitesimal vertical layer $(z,z+\mathrm {d} z)$ and all time steps (after reaching the steady state) is (Pähtz & Durán 2018b)

(2.3)\begin{equation} \langle A\rangle(z)=\left.\sum_{z_p\in(z,z+\mathrm{d} z)}m_pA_p\right/\sum_{z_p\in(z,z+\mathrm{d} z)}m_p, \end{equation}

where $m_p$ and $z_p$ are the particle mass and elevation, respectively. We also define the average of a vertical profile $\langle A\rangle (z)$ over the transport layer as (Pähtz & Durán 2018a)

(2.4)\begin{equation} \bar{A}=\left.\int_0^\infty\rho\langle A\rangle\,\mathrm{d} z\right/\int_0^\infty\rho\,\mathrm{d} z, \end{equation}

where $\rho$ is the local particle concentration. The bed surface elevation $z=0$ is defined as the elevation at which $p_g\,\mathrm {d}\langle v_x\rangle /\mathrm {d} z$ is maximal (Pähtz & Durán 2018b), where $\langle v_x\rangle$ is the average grain velocity in the streamwise direction and $p_g(z)=-\int _z^\infty \rho \langle a_z\rangle \,\mathrm {d} z^\prime$ the normal-bed granular pressure, with $\boldsymbol {a}$ the acceleration of grains by non-contact forces.

2.2. Calculation of transport rate and cessation threshold

We calculate the sediment transport rate $Q$ as (Pähtz & Durán 2018b)

(2.5)\begin{equation} Q=\int_{-\infty}^\infty\rho\langle v_x\rangle\,\mathrm{d} z. \end{equation}

When $Q$ vanishes, the grain-borne shear stress at the bed surface $\tau _g(0)$ also vanishes, with $\tau _g(z)=\int _z^\infty \rho \langle a_x\rangle \,\mathrm {d} z^\prime$ the grain-borne shear stress profile. We therefore extrapolate the cessation threshold value $\tau _t$ of the fluid shear stress $\tau$ at which $Q$ vanishes using the approximate relation (Pähtz & Durán 2018b)

(2.6)\begin{equation} \tau_g(0)=\tau-\tau_t, \end{equation}

where we treat $\tau _t$ as a fit parameter.

2.3. Dimensionless control parameters and rescaling of physical quantities

The average properties of equilibrium sediment transport are mainly determined by a few grain and environmental parameters: the grain and fluid density ($\rho _p$ and $\rho _f$, respectively), median grain diameter ($d$), kinematic fluid viscosity ($\nu$), fluid shear velocity ($u_\ast \equiv \sqrt {\tau /\rho _f}$) and gravitational constant ($g$) or its buoyancy-reduced value $\tilde {g}\equiv (1-\rho _f/\rho _p)g$ (for air, $\tilde {g}\simeq g$). Physical quantities with a superscript ‘$+$’ are rescaled using units of $\rho _p$, $\tilde {g}$ and $\nu$. For example,

(2.7)\begin{gather} d^+=\tilde{g}d/(\tilde{g}\nu)^{2/3}, \end{gather}

(2.8)\begin{gather}u^+_\ast=u_\ast{/}(\tilde{g}\nu)^{1/3}, \end{gather}

(2.9)\begin{gather}Q^+=Q/(\rho_p\nu). \end{gather}

As we show, this rescaling is well suited to describe the relevant physical processes underlying the cessation threshold scaling. A given environmental condition is fully determined by the values of three dimensionless numbers (Pähtz & Durán 2020):

(2.10)\begin{gather} s\equiv\rho_p/\rho_f=1/\rho_f^{+}, \end{gather}

(2.11)\begin{gather}Ga\equiv\sqrt{s\tilde{g}d^3}/\nu=\sqrt{s}d^{{+}3/2}, \end{gather}

(2.12)\begin{gather}\varTheta\equiv u_\ast^2/(s\tilde{g}d)=u_\ast^{{+}2}/(sd^+). \end{gather}

Numerical simulations are carried out for various combinations of the particle–fluid density ratio $s$ and Galileo number $Ga$, exceeding previously simulated conditions by almost four orders of magnitude in $s$ (table 1), and for Shields numbers $\varTheta$ ranging from weak conditions near its cessation threshold value $\varTheta _t$ to intense conditions far above $\varTheta _t$.

Table 1. Simulated particle–fluid density ratios $s$ and Galileo numbers $Ga$. $^\ast$The condition $s=2.5\times 10^5$, $Ga=1$ corresponds to a typical transport environment on Mars ($d\approx 100\ \mathrm {\mu } \mathrm {m}$) and $s=10^7$, $Ga=0.2$ to a hypothetical transport environment on Pluto ($d\approx 200\ \mathrm {\mu } \mathrm {m}$). Simulations with significantly larger respective values of $Ga$ are unstable for these large-$s$ conditions. We have been unable to fix this issue and do not know whether it has numerical or physical causes. The asterisk symbol, ${\dagger}$, indicates conditions simulated in our previous studies (Pähtz & Durán 2018a, 2020).

2.4. Sediment transport regimes for near-threshold conditions

Since the mixing length-based Reynolds-averaged description of hydrodynamics used in the numerical model neglects turbulent fluctuations around the mean turbulent flow, simulated sediment transport is always non-suspended. Near the cessation threshold (subscript $t$), non-suspended transport occurs as either bedload or saltation (see the introduction), which we distinguish through the criterion (Pähtz & Durán 2018a)

(2.13)\begin{equation} \text{Transport regime}= \begin{cases} \text{bedload} & \text{if}\ \overline{v_z^2}_t/\tilde{g}< d ,\\ \text{saltation} & \text{if}\ \overline{v_z^2}_t/\tilde{g}\geqslant d. \end{cases} \end{equation}

The quantity $\overline {v_z^2}/\tilde {g}$ describes the contribution of hopping grains to the characteristic transport height of all transported grains $\bar {z}$, where the latter also include those that role and slide. In particular, for saltation near the cessation threshold, $\overline {v_z^2}_t/\tilde {g}\simeq \bar {z}_t$, whereas $\overline {v_z^2}_t/\tilde {g}$ is significantly smaller than $\bar {z}_t$ for bedload transport (figure 1). Henceforth, $\overline {v_z^2}/\tilde {g}$ and $\bar {z}$ are termed hop height and transport layer thickness, respectively, for simplicity.

Figure 1. Transport layer thickness $\bar {z}_t$ versus hop height $\overline {v_z^2}_t/\tilde {g}$, both relative to the grain size $d$. Symbols correspond to numerical simulations near the cessation threshold for various combinations of the density ratio $s$ and Galileo number $Ga$ (see table 1), with open and filled symbols indicating bedload and saltation conditions, respectively.

3. Results

This section is organised as follows. First, it shows the data and scaling laws of the cessation threshold and equilibrium transport rate obtained from the simulations for the saltation regime (§ 3.1). Second, it presents semi-empirical physical justifications of these laws, including a first-principle-based proof of the statement that, under relatively mild assumptions, the rescaled cessation threshold $u_{\ast t}^+$ is a function of only one dimensionless control parameter (§ 3.2). Third, it tests existing models from the literature against the numerical data (§ 3.3). Fourth, it provides semi-empirical generalisations of the scaling laws that bridge between the saltation and bedload regimes (§ 3.4) and shows how they are affected by modifications of the drag law (§ 3.5), which may occur, for example, in highly rarefied atmospheres due to drag rarefaction.

3.1. Simulation data and scaling laws for saltation

3.1.1. Cessation threshold

Of the physical parameters affecting the shear velocity at the cessation threshold $u_{\ast t}$, the surface air pressure $P$ varies most strongly with the planetary environment. Furthermore, for a given planetary environment, the grain size $d$ is the most strongly varying relevant physical parameter. To isolate the effect of $P$ on $u_{\ast t}$, we normalise $u_{\ast t}$ in terms of relevant parameters that do neither depend on $P$ nor on $d$, $U_{\ast t}\equiv u_{\ast t}/(\mu g/\rho _p)^{1/3}$ (using that the dynamic viscosity $\mu =\rho _f\nu$ does not depend on $P$), and compare it with the density ratio $s$, which incorporates the effect of $P$ isolated from that of $d$.

For saltation, the simulations reveal a lower bound for $U_{\ast t}$ scaling as $s^{1/3}$ (figure 2, filled circles). This is distinct from the classical scaling of the saltation initiation threshold with $s^{1/2}$ (Greeley et al. 1976, 1980; Iversen & White 1982; Greeley et al. 1984; Burr et al. 2015, 2020; Swann et al. 2020) (figure 2, gray crosses), which follows from a balance between flow-induced and resisting forces or torques acting in bed surface grains (Pähtz et al. 2020). Roughly the same $s^{1/2}$-scaling was also found for the dynamic-threshold measurements by Andreotti et al. (2021) carried out in a low-pressure wind tunnel (figure 2, black crosses). As mentioned in the introduction and discussed in more detail in § 4, these measurements may constitute data of the continuous-transport threshold rather than the cessation threshold.

Figure 2. Cessation threshold shear velocity normalised using air-pressure- and grain-size-independent natural units $U_{\ast t}\equiv u_{\ast t}/(\mu g/\rho _p)^{1/3}$ versus density ratio $s$. Symbols that appear in the legend correspond to initiation (Greeley et al. 1976, 1980; Iversen & White 1982; Greeley et al. 1984; Burr et al. 2015, 2020; Swann, Sherman & Ewing 2020) and cessation threshold measurements for aeolian transport of quartz (Bagnold 1937; Martin & Kok 2018; Zhu et al. 2019), clay loam (Chepil 1945) and snow at sea level (Sugiura et al. 1998) and high altitude (Clifton, Rüedi & Lehning 2006, HA). The dynamic-threshold measurements by Andreotti et al. (2021) may constitute data of the continuous-transport threshold rather than the cessation threshold (discussed in § 4). Symbols that do not appear in the legend correspond to numerical simulations for various combinations of $s$ and the Galileo number $Ga$ (see table 1 and figure 1), with open and filled symbols indicating bedload and saltation conditions, respectively (see figure 1 for the definition). The solid line corresponds to $U_{\ast t}\propto s^{1/3}$ and represents the lower bound for cessation and initiation thresholds of saltation.

In addition to its $s^{1/3}$-scaling, $U_{\ast t}$ varies with the normalised median grain diameter $D_\ast \equiv \sqrt {s}d^+=\sqrt {s}d\tilde {g}/(\tilde {g}\nu )^{2/3}$, described by the following relationship between the rescaled cessation threshold $u_{\ast t}^+$ (note that $u_{\ast t}^+=U_{\ast t}/(s-1)^{1/3}$) and $D_\ast$:

(3.1)\begin{equation} u^+_{{\ast} t}=u^{{+}min}_{{\ast} t}\max\left[\left(\frac{D_\ast}{D^{min}_\ast}\right)^{{-}1/2},\left(\frac{D_\ast}{D^{min}_\ast}\right)^{1/2}\right]. \end{equation}

It contains the parameters $D^{min}_\ast$ and $u^{+min}_{\ast t}$, which denote the location and magnitude, respectively, of the minimum of the function $u_{\ast t}^+(D_\ast )$, corresponding to the lower bound of $U_{\ast t}$ for saltation in figure 2. Equation (3.1) is consistent with the simulations (figure 3a) and experiments (figure 3b) for the saltation regime, though with slightly different parameter values: $(D^{min}_\ast,u^{+min}_{\ast t})=(16,1.6)$ versus $(D^{min}_\ast,u^{+min}_{\ast t})=(18,2.3)$, respectively. The associated relative change of $u_{\ast t}^+$ by $2.3/1.6\simeq 1.4$ is well within the typical systematic uncertainty of cessation threshold measurements. For example, Creyssels et al. (2009) reported $\varTheta _t=0.009$ for their terrestrial wind tunnel experiments ($d=242\ \mathrm {\mu } \mathrm {m}$), obtained from extrapolating transport rate measurements to vanishing transport using the transport rate model of Ungar & Haff (1987), whereas Pähtz & Durán (2020) reported $\varTheta _t=0.0035$ for the very same data using a different transport rate model for the extrapolation, resulting in a relative change of $\sqrt {0.009/0.0035}\simeq 1.6$.

Figure 3. Rescaled cessation threshold shear velocity $u^+_{\ast t}$ versus normalised median grain diameter $D_\ast \equiv \sqrt {s}d^+$. Symbols in (a) correspond to numerical simulations of saltation (see figure 1 for the definition) for various combinations of the density ratio $s$ and Galileo number $Ga$ (see table 1 and figure 1). Symbols in (b) correspond to experimental cessation threshold data (see legend) for terrestrial aeolian saltation of quartz (Bagnold 1937; Martin & Kok 2018; Zhu et al. 2019), clay loam (Chepil 1945) and snow at sea level (Sugiura et al. 1998) and high altitude (Clifton et al. 2006, HA). The solid lines correspond to (3.1), with $(D_\ast ^{min},u^{+min}_{\ast t})=(16,1.6)$ in (a) and $(D_\ast ^{min},u^{+min}_{\ast t})=(18,2.3)$ in (b).

3.1.2. Equilibrium transport rate

The simulations of saltation and experiments reasonably collapse on the master curve (figure 4)

(3.2)\begin{equation} Q^+{/}d^{{+}3/2}=1.7s^{1/3}(\varTheta-\varTheta_t)+12s^{1/3}(\varTheta-\varTheta_t)^2 \end{equation}

if $Ga\sqrt {s}>81$. The vast majority of planetary transport occurring in nature and most of the simulated saltation conditions satisfy this criterion. Note that $Ga\sqrt {s}$ can be interpreted as a Stokes-like number (Berzi et al. 2016), encoding the importance of grain inertia relative to viscous drag forcing, and controls the transition to viscous bedload (Pähtz et al. 2021).

Figure 4. Normalised sediment transport rate $s^{-1/3}Q^+/d^{+3/2}$ versus Shields number in excess of the cessation threshold $\varTheta -\varTheta _t$. Symbols in (a) correspond to numerical simulations of saltation (see figure 1 for the definition) for various combinations of the density ratio $s$ and Galileo number $Ga$ (see table 1 and figure 1) with $Ga\sqrt {s}>81$, and Shields number $\varTheta$. Symbols in (b) correspond to measurements for different grain sizes (indicated in the legend) for terrestrial aeolian saltation of minerals (Creyssels et al. 2009; Ho et al. 2011; Ho 2012; Martin & Kok 2017; Ralaiarisoa et al. 2020) and snow (Sugiura et al. 1998). The values of $\varTheta _t$ in (b) for a given experimental data set are obtained from extrapolating (3.2) to vanishing transport. Note that Ralaiarisoa et al. (2020) reported that transport may not have been completely in equilibrium in their experiments. The solid lines correspond to (3.2).

3.2. Physical justifications of saltation scaling laws

3.2.1. First-principle-based proof that $u_{\ast t}^+$ is a function of only $D_\ast$

In general, the shear velocity at the cessation threshold $u_{\ast t}$ is a function of the five control parameters $\rho _p$, $\rho _f$, $\nu$, $\tilde {g}$ and $d$ (Claudin & Andreotti 2006). These parameters involve three units (mass, length and time). According to the $\varPi$ theorem (Barenblatt 1996), the physical system, and therefore any dimensionless system property such as $u_{\ast t}^+$, is then controlled by two dimensionless numbers, for example the density ratio $s$ and the normalised median grain diameter $D_\ast$:

(3.3)\begin{equation} u^+_{{\ast} t}=f(s,D_\ast). \end{equation}

To determine the function $f$ in (3.3), existing cessation threshold models have made various idealisations of the fluid–particle system (Claudin & Andreotti 2006; Kok 2010b; Berzi et al. 2016, 2017; Pähtz & Durán 2018a; Andreotti et al. 2021; Pähtz et al. 2021; Gunn & Jerolmack 2022). In particular, they all drastically coarse-grain the particle phase of the aeolian transport layer above the bed surface, either by representing the entire grain motion by identical periodic saltation trajectories (Claudin & Andreotti 2006; Kok 2010b; Berzi et al. 2016, 2017; Pähtz & Durán 2018a; Andreotti et al. 2021; Pähtz et al. 2021; Gunn & Jerolmack 2022) or by an average motion behaviour (Kok 2010b; Pähtz & Durán 2018a) that is mathematically equivalent to an identical periodic trajectory representation (Pähtz et al. 2020).

Here, in contrast to previous models, we do not resort to any such coarse-graining. Instead, we idealise the system in the following comparably mild manner.

1. (i) We consider only buoyancy and Stokes drag as fluid–grain interactions, neglecting form drag contributions. This would be justified if relatively fast saltating grains dominated the near-threshold grain dynamics, since comparably faster grains exhibit comparably lower fluid–particle velocity differences and, thus, comparably less form drag relative to Stokes drag.

2. (ii) Due to the typically relatively small shear Reynolds numbers associated with planetary transport near the cessation threshold, $Ga\sqrt {\varTheta _t}\lesssim 10$, we consider a smooth inner turbulent boundary layer mean flow velocity profile $u_x(z)$, neglecting hydrodynamically rough contributions (and turbulent fluctuations, which are also neglected in the numerical simulations).

3. (iii) Since vanishingly few grains are in motion sufficiently close to the cessation threshold, we neglect the feedback of the grain motion on the flow.

4. (iv) Since saltation trajectories are typically much larger than the grain size, we consider an idealised flat bed and assume that the zero level of the flow velocity coincides with the grain elevation at grain–bed impact ($z=0$), neglecting the effect of the flow very near the bed surface to the overall grain motion.

5. (v) While we do not specify the distribution of grain lift-off velocities $f_\uparrow$ and grain impact velocities $f_\downarrow$, we assume that the boundary conditions mapping $f_\downarrow$ to $f_\uparrow$ in the steady state are scale-free, as for grain–bed rebounds (Beladjine et al. 2007), neglecting the potential effect of $\sqrt {\tilde {g}d}$ on grain–bed collisions. Most grains ejected by the splash of a grain impacting the bed with velocity $\boldsymbol {v}_\downarrow$ exhibit a velocity on the order of $\sqrt {\tilde {g}d}$ and only the few grains corresponding to the upper-tail end of the distribution exhibit an ejection velocity proportional to $|\boldsymbol {v}_\downarrow |$ (Lämmel et al. 2017). Hence, this assumption effectively means that grain–bed rebounds and/or rare extreme ejection events dominate the saltation dynamics relevant for the cessation threshold scaling.

Under the above assumptions, the equations of motion for a given grain are (Pähtz et al. 2021)

(3.4)\begin{gather} \dot{v}^+_z={-}1-v^+_z/v^+_s, \end{gather}

(3.5)\begin{gather}\dot{v}^+_x=(u^+_x-v^+_x)/v^+_s, \end{gather}

(3.6)\begin{gather}u^+_x=u^+_\ast f_u(u^+_\ast z^+), \end{gather}

where $\boldsymbol {v}^+$ is the rescaled grain velocity, $v_s^+=4sd^{+2}/(3Re_c)$ the rescaled Stokes settling velocity (obtained from the high-viscosity limit of (2.1)), and $f_u(X)$ denotes a function describing $u_x/u_\ast$ for an undisturbed smooth inner turbulent boundary layer. It obeys $f_u(X)=X$ within the viscous sublayer of the turbulent boundary layer ($X\lesssim 5$) and $f_u(X)\simeq \kappa ^{-1}\ln (9X)$ within its log-layer ($X\gtrsim 30$). Extrapolated into the transitional buffer layer in between, both profiles would intersect at about $X=11$, which is why $\delta _\nu =11\nu /u_\ast$ is termed viscous-sublayer thickness.

Parametrised by $v^+_s$ and $u^+_\ast$, (3.4)–(3.6) map $f_\uparrow$ to $f_\downarrow$. Combined with the scale-free boundary conditions, mapping $f_\downarrow$ back to $f_\uparrow$, they imply that the grain motion is fully determined by $v^+_s$ and $u^+_\ast$. For a given $v^+_s$, the cessation threshold $u^+_{\ast t}$ then corresponds to the smallest value of $u^+_\ast$ for which a solution of the combined system exists (Pähtz et al. 2021). This implies that there is a function $f$ mapping $D_\ast =\sqrt {18v^+_s}$ (valid for $Re_c=24$, the standard case of non-rarefied drag) to $u^+_{\ast t}$:

(3.7)\begin{equation} u^+_{{\ast} t}=f(D_\ast). \end{equation}

In summary, the above assumptions simplify the general two-parametric dependence of $u_{\ast t}^+$ in (3.3) to the one-parametric dependence in (3.7), in agreement with (3.1).

3.2.2. Semi-empirical model of cessation threshold scaling

While the above analysis explains why $u_{\ast t}^+=f(D_\ast )$ in (3.1), it does not yield the function $f$ itself. Here, we derive the expression for $f$ in (3.1) guided by the simulations. The latter show that the minimum $u^{+min}_{\ast t}$ for saltation occurs when the hop height $\overline {v_z^2}_t/\tilde {g}$ is about equal to the viscous-sublayer thickness $\delta _{\nu t}=11\nu /u_{\ast t}$ near the cessation threshold (figure 5a). This can be explained using the empirical, yet physically reasonable, simulation-supported proportionality between the average fluid velocity $\overline {u_x}_t$ and $(\overline {v_z^2})^{1/2}_t$ near the cessation threshold (figure 5b). In fact, averaging (3.6) over all grain trajectories and the transport layer, using the approximation $\overline {f_u(u^+_{\ast t}z^+)}_t\simeq f_u(u^+_{\ast t}\overline {z^+}_t)$, and using this proportionality approximately yields for saltation ($\overline {z^+}_t\simeq \overline {v_z^{+2}}_t$, see figure 1):

(3.8)\begin{equation} u^+_{{\ast} t}\propto\left[\frac{u^+_{{\ast} t}\overline{v_z^{{+}2}}_t}{f_u^2\left(u^+_{{\ast} t}\overline{v_z^{{+}2}}_t\right)}\right]^{1/3}. \end{equation}

Within the viscous sublayer ($u^+_{\ast t}\overline {v_z^{+2}}_t\lesssim 5$), this relation simplifies to $u^+_{\ast t}\propto (u^+_{\ast t}\overline {v_z^{+2}}_t)^{-1/3}$ and within the log-layer approximately to $u^+_{\ast t}\propto (u^+_{\ast t}\overline {v_z^{+2}}_t)^{1/3}$, neglecting the logarithmic term. The crossover between the two power laws occurs about at $u^+_{\ast t}\overline {v_z^{+2}}_t=11$, that is, when the hop height exceeds the viscous-sublayer thickness ($\overline {v_z^2}_t/\tilde {g}=\delta _{\nu t}$). Hence, the parabolic law

(3.9)\begin{equation} u^+_{{\ast} t}=u^{{+}min}_{{\ast} t}\max\left[\left(\frac{\overline{v_z^2}_t}{\tilde{g}\delta_{\nu t}}\right)^{{-}1/3},\left(\frac{\overline{v_z^2}_t}{\tilde{g}\delta_{\nu t}}\right)^{1/3}\right] \end{equation}

fits the saltation data reasonably well (solid line in figure 5a).

Figure 5. (a) Rescaled threshold shear velocity $u^+_{\ast t}$ versus ratio between hop height $\overline {v_z^2}_t/\tilde {g}$ and viscous-sublayer thickness $\delta _{\nu t}=11\nu /u_{\ast t}$ near the cessation threshold. (b) Rescaled transport layer-averaged fluid velocity $\overline {u^+_x}_t$ versus $(\overline {v_z^{+2}})^{1/2}_t$ near the cessation threshold. (c) Plot of $u^+_{\ast t}\overline {v_z^{+2}}_t$ versus $v^+_s$. Symbols correspond to numerical simulations of saltation (see figure 1 for the definition) for various combinations of the density ratio $s$ and Galileo number $Ga$ (table 1 and figure 1). The solid lines in (a), (c) and (b) correspond to (3.9), (3.10) and $\overline {u^+_x}_t=6(\overline {v_z^{+2}})^{1/2}_t$, respectively.

Following from the analysis we have used to deduce (3.7), the grain kinematics near the cessation threshold, and thus $\overline {v_z^{+2}}_t$, should be controlled by $u^+_{\ast t}$ or $v^+_s$. Indeed, the simulations of saltation suggest the empirical relation (figure 5c)

(3.10)\begin{equation} u^+_{{\ast} t}\overline{v_z^{{+}2}}_t=1.5v_s^{{+}3/4}, \end{equation}

which leads to (3.1) with $D_\ast ^{min}=\sqrt {18}(11/1.5)^{2/3}\simeq 16$.

According to the above model, the grain size scaling of $u^+_{\ast t}$ in (3.1), despite being mathematically equivalent to the well-known cohesive ($u_{\ast t}\sim d^{-1/2}$, left branch) and cohesionless ($u_{\ast t}\sim d^{1/2}$, right branch) limits of the saltation initiation threshold (Shao & Lu 2000), follows purely from hydrodynamics rather than the onset of cohesion at small grain sizes.

3.2.3. Physics behind equilibrium transport rate scaling

Analytical, physical models of the equilibrium transport rate $Q$ for aeolian saltation typically separate it into the mass of transported sediment per unit area of the bed $M$ and its average streamwise velocity $V$ through $Q=MV$. In most models, it is reasoned that the scaling of $V$ is in one way or another linked to grain–bed collisions, and since the average outcome of grain–bed collisions should be roughly independent of the wind speed at equilibrium, $V$ is taken as equal to its near-threshold value $V_t$ (Ungar & Haff 1987; Durán et al. 2011; Kok et al. 2012; Berzi et al. 2016). However, it has been shown that, for sufficiently intense saltation, midair collisions significantly disturb grain trajectories (Carneiro et al. 2013; Pähtz & Durán 2020; Ralaiarisoa et al. 2020), leading to an additional additive term increasing as $M^+/d^+$ (Pähtz & Durán 2020):

(3.11)\begin{equation} Q^+=M^+V^+_t(1+c_MM^+{/}d^+), \end{equation}

where $c_M$ is a constant parameter. It is not trivial to evaluate the scalings of $M^+$ and $V^+_t$ with the simulation data, since extracting $M$ and $V$ from DEM-based numerical transport simulations is ambiguous (Durán et al. 2012; Pähtz & Durán 2018b). One possible way is to define $M$ as the mass $M_0$ of grains moving above the bed surface ($z=0$) per unit bed area and $V$ as their average streamwise velocity (Pähtz & Durán 2018b):

(3.12)\begin{gather} M\equiv\int_0^\infty\rho\,\mathrm{d} z=M_0, \end{gather}

(3.13)\begin{gather}V\equiv\frac{\displaystyle\int\nolimits_0^\infty\rho\langle v_x\rangle\,\mathrm{d} z}{ \displaystyle\int\nolimits_0^\infty\rho\,\mathrm{d} z}=\overline{v_x}. \end{gather}

This definition uses that most (but not all) sediment transport occurs at elevations $z>0$, especially for saltation and, therefore, $M_0\overline {v_x}=\int _0^\infty \rho \langle v_x\rangle \,\mathrm {d} z\simeq \int _{-\infty }^\infty \rho \langle v_x\rangle \,\mathrm {d} z=Q$ (Pähtz & Durán 2018b). Alternatively, one can define $V$ as the mass flux-weighted average $\overline {v_x}^q$ of the streamwise velocity of all grains and $M_q$, the associated value of $M$, as $M_q\equiv Q/\overline {v_x}^q$ (Durán et al. 2012):

(3.14)\begin{gather} M\equiv\frac{\left(\displaystyle\int\nolimits_{-\infty}^\infty\rho\langle v_x\rangle\,\mathrm{d} z\right)^2}{ \displaystyle\int\nolimits_{-\infty}^\infty\rho\langle v_x^2\rangle\,\mathrm{d} z}=M_q, \end{gather}

(3.15)\begin{gather}V\equiv\frac{\displaystyle\int\nolimits_{-\infty}^\infty\rho\langle v_x^2\rangle\,\mathrm{d} z}{ \displaystyle\int\nolimits_{-\infty}^\infty\rho\langle v_x\rangle\,\mathrm{d} z}=\overline{v_x}^q, \end{gather}

where $\bar {{\cdot }}^q\equiv ({1}/{Q})\int _{-\infty }^\infty \rho \langle v_x{\cdot }\rangle \,\mathrm {d} z$.

For the above two definitions of $M$ and $V$, the simulations are roughly described by scaling laws in which a comparably small part of the $s^{1/3}$-scaling factor in (3.2) goes into $M^+/d^+$ and a comparably large part into $V^+_t/\sqrt {d^+}$ (figure 6). However, the exact partitioning of $s^{1/3}$ depends on the chosen definition (figures 6(a) and 6(b) versus figures 6(c) and 6(d)):

(3.16a,b)\begin{gather} M^+_0\propto s^{1/12}d^+(\varTheta-\varTheta_t),\quad\overline{v^+_x}_t\propto s^{1/4}\sqrt{d^+}, \end{gather}

(3.17a,b)\begin{gather}M^+_q\propto d^+(\varTheta-\varTheta_t),\quad\overline{v^+_x}^q_t\propto s^{1/3}\sqrt{d^+}. \end{gather}

The latter scaling is consistent with the prediction $M^+\propto d^+(\varTheta -\varTheta _t)$ from physical models (Ungar & Haff 1987; Durán et al. 2011; Berzi et al. 2016; Pähtz & Durán 2020) and with (3.2) when combined with (3.11). However, it means that $V^+_t\propto s^{1/3}\sqrt {d^+}$, which is a highly unusual scaling, different from the existing models $V^+_t\propto \sqrt {d^+}$ (Ungar & Haff 1987; Berzi et al. 2016) and $V^+_t\propto u^+_{\ast t}$ (Durán et al. 2011; Kok et al. 2012; Pähtz & Durán 2020).

Figure 6. (a,c) Normalised transport loads $s^{-1/12}M^+_0/d^+$ and $M^+_q/d^+$, using the definitions (3.12) and (3.14), respectively, of $M$; and (b,d) normalised average streamwise grain velocities $s^{-1/4}\overline {v^+_x}/\sqrt {d^+}$ and $s^{-1/3}\overline {v^+_x}^q/\sqrt {d^+}$, using the definition (3.13) and (3.15), respectively, of $V$ versus Shields number in excess of the cessation threshold $\varTheta -\varTheta _t$. Symbols correspond to numerical simulations of saltation (see figure 1 for the definition) for various combinations of the density ratio $s$ and Galileo number $Ga$ (see table 1 and figure 1) with $Ga\sqrt {s}>81$, and Shields number $\varTheta$. The solid lines in (a,b) correspond to the left equations in (3.16a,b) and (3.17a,b), respectively.

3.3. Test of existing models against simulations of saltation

3.3.1. Test of cessation threshold models

The most important assumption that led to the simulation-supported statement that the rescaled cessation threshold $u_{\ast t}^+$ is solely controlled by the normalised median grain diameter $D_\ast$ in § 3.2.1 is that of scale-free boundary conditions. The only existing cessation threshold model with scale-free boundary condition is that of Pähtz et al. (2021), which we here compare with the most recent alternative, that of Gunn & Jerolmack (2022). The latter's most important feature is that it superimposes a $Ga$-dependent damping on the scale-free laws describing grain–bed rebounds, where the damping function is essentially fitted to agreement with experimental cessation threshold data. We find that, while the model of Pähtz et al. (2021) captures the simulation data very well, the model of Gunn & Jerolmack (2022), with its drag and lift laws being modified to those employed in the simulations (i.e. (2.1) and no lift) for a fair comparison, is in very strong disagreement (figure 7). This is discussed in § 4.

Figure 7. Evaluation of the cessation threshold models of (a) Pähtz et al. (2021) and (b) Gunn & Jerolmack (2022), where the latter's drag and lift laws have been modified to those employed in the simulations for a fair comparison. Rescaled cessation threshold $u_{\ast t}^+$ versus normalised median grain diameter $D_\ast \equiv \sqrt {s}d^+$. Symbols correspond to numerical simulations of saltation (see figure 1 for the definition) for various combinations of the density ratio $s$ and Galileo number $Ga$ (see table 1 and figure 1). The solid lines indicate the respective model predictions. Their color characterises $s$ in accordance with the symbol color.

3.3.2. Test of equilibrium transport rate models

The simulations of saltation are not or not well captured by the two most widely used physical models of the equilibrium aeolian transport rate: the model of Ungar & Haff (1987) and others (e.g. Jenkins & Valance 2014; Berzi et al. 2016), $Q^+/d^{+3/2}=f_1(\varTheta -\varTheta _t)$ (figure 8a) and the model of Durán et al. (2011) and others (Kok et al. 2012; Pähtz & Durán 2020), $Q^+/(d^+u^+_{\ast t})=f_2(\varTheta -\varTheta _t)$ (figure 8b).

Figure 8. Evaluation of the physically based functional relationships for the sediment transport rate by Ungar & Haff (1987) and Durán et al. (2011). Normalised sediment transport rate (a) $Q^+/d^{+3/2}$ and (b) $Q^+/(d^+u^+_{\ast t})$ versus Shields number in excess of the cessation threshold $\varTheta -\varTheta _t$. Symbols correspond to numerical simulations of saltation (see figure 1 for the definition) for various combinations of the density ratio $s$ and Galileo number $Ga$ (see table 1 and figure 1) with $Ga\sqrt {s}>81$, and Shields number $\varTheta$.

3.4. Generalised scaling laws across saltation and turbulent bedload

It is possible to semi-empirically generalise (3.1) to also include turbulent bedload conditions, defined by $s\lesssim 10$ and $D_\ast \gtrsim D_\ast ^{min}$ (equivalent to $Gas^{1/4}\gtrsim 64$, which ensures that transported grains significantly penetrate the log-layer; Pähtz & Durán 2020). Turbulent bedload includes not only hopping grains but also rolling grains. The threshold shear velocity required to sustain a pure, very slow rolling motion along the bed surface scales as $u_{\ast t}\propto \sqrt {s\tilde {g}d}$ (Pähtz et al. 2021), which corresponds to $u^+_{\ast t}\propto s^{1/4}$ at the cessation threshold minimum $D_\ast =D_\ast ^{min}$. We find that the empirical relation $u^+_{\ast t}=\sqrt {f_s}u^{+min}_{\ast t}$, with $f_s\equiv (1+\sqrt {10/s})^{-1}$, captures the transition from $u^+_{\ast t}\propto s^{1/4}$ for $s\ll 10$ to $u^+_{\ast t}=u^{+min}_{\ast t}$ for $s\gg 10$ at $D_\ast =D_\ast ^{min}$. The resulting generalised cessation threshold equation is

(3.18)\begin{equation} u^+_{{\ast} t}=\sqrt{f_s}u^{{+}min}_{{\ast} t}\max\left[\left(\frac{D_\ast}{D^{min}_\ast}\right)^{{-}1/2},\left(\frac{D_\ast}{D^{min}_\ast}\right)^{1/2}\right]. \end{equation}

It is consistent with the simulations and experiments across aeolian and fluvial transport conditions with $Ga\sqrt {s}\gtrsim 81$ (figure 9).

Figure 9. Bedload-corrected rescaled cessation threshold shear velocity $u^+_{\ast t}/\sqrt {f_s}$ versus normalised median grain diameter $D_\ast \equiv \sqrt {s}d^+$. Symbols in (a) correspond to numerical simulations for various combinations of the density ratio $s$ and Galileo number $Ga$ (see table 1 and figure 1) with $Ga\sqrt {s}>81$, where open and filled symbols indicate bedload and saltation conditions, respectively (see figure 1 for the definition). Symbols in (b) correspond to experimental cessation threshold data (see the legend) for terrestrial aeolian saltation of quartz (Bagnold 1937; Martin & Kok 2018; Zhu et al. 2019), clay loam (Chepil 1945) and snow at sea level (Sugiura et al. 1998) and high altitude (Clifton et al. 2006, HA), and a compilation of experimental threshold data for subaqueous bedload (Buffington & Montgomery 1997). Only data with $Ga\sqrt {s}>81$ are shown. The solid lines correspond to (3.1), with $(D_\ast ^{ min},u^{+min}_{\ast t})=(16,1.6)$ in (a) and $(D_\ast ^{min},u^{+min}_{\ast t})=(18,2.3)$ in (b).

Furthermore, an empirical generalisation of (3.2) to turbulent bedload conditions is given by

(3.19)\begin{equation} Q^+{/}d^{{+}3/2}=1.7s^{1/3}(\varTheta-\varTheta_t)+13f_ss^{1/3}(\varTheta-\varTheta_t)^2, \end{equation}

consistent with the simulations and experiments across aeolian and fluvial transport conditions with $Ga\sqrt {s}\gtrsim 81$ (figure 10).

Figure 10. Bedload-corrected normalised sediment transport rate $f_ss^{-1/3}Q^+/d^{+3/2}$ versus bedload-corrected Shields number in excess of the cessation threshold $f_s(\varTheta -\varTheta _t)$. Symbols in (a) correspond to numerical simulations for various combinations of the density ratio $s$ and Galileo number $Ga$ (see table 1 and figure 1) with $Ga\sqrt {s}>81$, and Shields number $\varTheta$, where open and filled symbols indicate bedload and saltation conditions, respectively (see figure 1 for the definition). Symbols in (b) correspond to measurements for different grain sizes (indicated in the legend) for terrestrial aeolian saltation of minerals (Creyssels et al. 2009; Ho et al. 2011; Ho 2012; Martin & Kok 2017; Ralaiarisoa et al. 2020) and snow (Sugiura et al. 1998), and subaqueous bedload (Meyer-Peter & Müller 1948; Smart & Jaeggi 1983; Capart & Fraccarollo 2011). We corrected the raw laboratory data by Smart & Jaeggi (1983) and Capart & Fraccarollo (2011) for sidewall drag using the method described in § 2.3 of Guo (2015) and for steep bed slopes using $u_\ast ^2|_{corrected}=u_\ast ^2/f_\alpha$, with $f_\alpha \equiv 1-\tan \alpha /0.63$ (Pähtz et al. 2021). The values of $\varTheta _t$ in (b) for a given experimental data set are obtained from extrapolating (3.19) to vanishing transport. Note that Ralaiarisoa et al. (2020) reported that transport may not have been completely in equilibrium in their experiments. The solid lines correspond to (3.19).

Put together, (3.18) and (3.19) can be used to predict the equilibrium transport rate for arbitrary combinations of the density ratio $s$, Galileo number $Ga$ and Shields number $\varTheta$ with $Ga\sqrt {s}\gtrsim 81$ for non-rarefied drag. When compared with the simulations, these equations perform significantly better than the unified model of the cessation threshold and equilibrium transport rate of Pähtz et al. (2021) (figure 11). While the latter captures the $s^{1/3}$-dependence of $Q^+$, it fails to capture the $d^+$-dependence of $Q^+/s^{1/3}$ observed in the simulations.

Figure 11. Predicted versus simulated sediment transport rate $Q^+$. (a) Predictions by (3.1) with $(D^{min}_\ast,u^{+min}_{\ast t})=(16,1.6)$ and (3.19). (b) Predictions by the model of Pähtz et al. (2021). Symbols correspond to numerical simulations for various combinations of the density ratio $s$ and Galileo number $Ga$ (see table 1 and figure 1) with $Ga\sqrt {s}>81$, where open and filled symbols indicate bedload and saltation conditions, respectively (see figure 1 for the definition). The solid lines indicate perfect agreement. The dashed lines indicate a deviation by a factor of two.

3.5. Effect of drag law and generalisation to drag in rarefied atmospheres

The analysis in § 3.2.1 suggests that the normalised median grain diameter $D_\ast \equiv \sqrt {s}d^+$ in (3.1) and (3.18) should be redefined as $D_\ast \equiv \sqrt {18v_s^+}=\sqrt {24s/Re_c}d^+$ (from (2.1)), which is equal to $\sqrt {s}d^+$ only in the case of non-rarefied drag ($Re_c=24$). To test this prediction as well as the effect of the form drag coefficient $C_d^\infty$, we carried out additional simulations using $C_d^\infty =0$ and $Re_c=[6,24,96]$ for a few saltation conditions. We find that these simulations, indeed, still satisfy (3.18) and therefore (3.1) when the redefined $D_\ast$ is used (figure 12a). They also still satisfy (3.19) and therefore (3.2), which are not affected by the redefinition of $D_\ast$ (figure 12b).

Figure 12. (a) Bedload-corrected rescaled cessation threshold shear velocity $u^+_{\ast t}/\sqrt {f_s}$ versus normalised median grain diameter, here redefined as $D_\ast \equiv \sqrt {18v_s^+}=\sqrt {24s/Re_c}d^+$. (b) Bedload-corrected normalised sediment transport rate $f_s s^{-1/3} Q^+/d^{+3/2}$ versus bedload-corrected Shields number in excess of the cessation threshold $f_s (\varTheta - \varTheta_t)$. Non-yellow symbols correspond to numerical simulations, carried out using the standard drag law parameters $Re_c=24$ and $C_d^\infty =0.5$, for various combinations of the density ratio $s$ and Galileo number $Ga$ (see table 1 and figure 1) with $Ga\sqrt {s}>81$ and Shields number $\varTheta$, where open and filled symbols indicate bedload and saltation conditions, respectively (see figure 1 for the definition). Yellow symbols correspond to additional simulations with modified drag law parameters as indicated in the legend.

In rarefied atmospheres, the mean free path $\lambda$ of the air molecules becomes comparable to the median grain diameter $d$, or the Knudsen number $Kn\equiv \lambda /d=\sqrt {{\rm \pi} k/2}s/(c^+d^+)$ (Crowe et al. 2012), with $c^+=c/(\tilde {g}\nu )^{1/3}$ the rescaled speed of sound and $k$ the adiabatic exponent, comparable to unity. This leads to a $Kn$-dependent correction $f_{Kn}\equiv 1+Kn[2.49+0.84\exp (-1.74/Kn)]$ (Crowe et al. 2012) of Stokes drag via $Re_c=24/f_{Kn}$. Note that typically $f_{Kn}\simeq 1$ for $s\lesssim 10^6$. Hence, the results in figure 12 support that the following generalised definition of $D_\ast$ should be used for highly rarefied atmospheres ($s\gtrsim 10^6$):

(3.20)\begin{equation} D_\ast\equiv\sqrt{f_{Kn}}\sqrt{s}d^+. \end{equation}

4. Discussion

4.1. Choice of dynamic-threshold measurements for evaluation of cessation threshold models

Equilibrium saltation becomes intermittent below the continuous-transport threshold, characterised by alternating periods of equilibrium saltation and periods of rest (Martin & Kok 2018). The cessation threshold is therefore the wind strength at which equilibrium saltation would cease if extrapolated from the continuous-transport regime, that is, as the zero-point of equilibrium transport equations such as (2.6) or (3.2). It is also the threshold of intermittent saltation (Martin & Kok 2018). To evaluate the cessation threshold scaling law in (3.1), we have therefore chosen exclusively measurements that either extrapolate continuous saltation in some manner to vanishing transport (Clifton et al. 2006; Zhu et al. 2019) or directly measure the cessation of intermittent saltation (Bagnold 1937; Chepil 1945; Sugiura et al. 1998; Martin & Kok 2018; Zhu et al. 2019). Both methods require that equilibrium transport conditions can at least temporarily be established during the experiments (Pähtz et al. 2020), usually by feeding sufficient sediment when the test section is too short for transport to reach equilibrium. This requirement was probably not satisfied in all of the above-cited measurements. The snow drift wind tunnel by Clifton et al. (2006), who did not feed snow at the tunnel entrance, was probably too short to establish equilibrium conditions for their beds of old and therefore cohesive snow, since cohesion can dramatically increase the fetch required to reach equilibrium (Comola et al. 2019a). For this reason, we have only compared with their data for freshly fallen snow.

Unfortunately, many other studies have not employed the same criteria when choosing measurements to evaluate their cessation threshold models (Claudin & Andreotti 2006; Kok 2010b; Berzi et al. 2017; Andreotti et al. 2021; Gunn & Jerolmack 2022). This has largely been driven by the belief that there is only one dynamic threshold, implying that any kind of dynamic-threshold measurement is at least a proxy for the cessation threshold. However, we have presented evidence for the hypothesis that the continuous-transport threshold is a second kind of dynamic threshold with an underlying physics different from that of the cessation threshold (Pähtz & Durán 2018a; Pähtz et al. 2020; Pähtz et al. 2021). An important example for a potential misinterpretation of measured dynamic thresholds as cessation thresholds is the study by Andreotti et al. (2021) for the following reasons.

1. (i) Andreotti et al. (2021), who carried out their measurements in a pressurised-wind tunnel, explicitly mentioned that they were only able to establish equilibrium transport for air pressures relatively close to ambient pressure ($P\gtrsim 30\,000\ \mathrm {Pa}$) but not for the vast majority of studied pressure conditions (down to $P\approx 200\ \mathrm {Pa}$): ‘below $300\ \mathrm {hPa}$ [the erosional zone] encompasses the entire bed’.

2. (ii) Andreotti et al. (2021) explicitly defined threshold conditions ‘as the transition between saltation of groups of particles (bursts) to intermittent saltation of single particles (at high pressure) or no transport (at low pressure)’. For high-pressure conditions, the so measured threshold is, by definition, larger than the cessation threshold (i.e. the threshold of intermittent saltation; Martin & Kok 2018). For low-pressure conditions, the measurements are difficult to interpret due to the lack of equilibrium transport.

3. (iii) Andreotti et al. (2021) accompanied their direct threshold measurements with indirect measurements obtained from extrapolating to vanishing transport. However, since they have not established equilibrium (for most pressure conditions), this extrapolation does not necessarily yield the cessation threshold.

4. (iv) Pähtz et al. (2021) hypothesised that the continuous-transport threshold is the smallest wind shear stress at which an average grain ejected by an impacting grain can be accelerated into a steady trajectory. A modification of their trajectory-based model based on this hypothesis captured the measurements by Andreotti et al. (2021), suggesting that their employed experimental method yields a threshold akin to the continuous-transport threshold.

The potential misinterpretation of the measurements by Andreotti et al. (2021) as cessation threshold measurements is highly relevant, since it led Gunn & Jerolmack (2022) to introduce a Galileo number ($Ga$)-dependent viscous damping of grain–bed rebounds in their cessation threshold model in an attempt to capture these data. However, this rebound damping is the very reason for the very strong disagreement between their model and the here presented numerical data of the cessation threshold (§ 3.3.1). Note that, from a physical perspective, rebound damping should not depend on $Ga$ but on the Stokes number associated with the grain's impact velocity $St=\rho _p|\boldsymbol {v_\downarrow }|d/\mu$ (Berzi et al. 2016, 2017; Andreotti et al. 2021), which is experimentally known to control the viscous damping of frontal grain collisions with a flat plate (Gondret, Lance & Petit 2002). Since typical values of $|\boldsymbol {v}_{\downarrow }|/\sqrt {\tilde {g}d}$ for Martian saltation are at the very least comparable to, if not much larger than, those for terrestrial saltation (because of $V_t\propto s^{1/3}\sqrt {\tilde {g}d}$, see § 3.2.3), and since $\rho _p$, $d$ and $\mu$ are of the same order of magnitude on Earth and Mars, typical values of $St$ on Mars are many orders of magnitude too large for viscous damping to play a meaningful role. In addition, even if there was a strong damping of frontal grain–plate collisions, this would not necessarily translate into a strong damping of grain–bed collisions. In fact, we previously reported only slight differences between DEM-RANS simulations of saltation for undamped (normal restitution coefficient $e=0.9$) and nearly fully damped ($e=0.01$) frontal grain–grain collisions (Pähtz & Durán 2018a). Even for $e=0.01$, grains can saltate in large hops on the order of $100d$ high (Movie S3 of Pähtz & Durán 2018a).

4.2. Recommendations for how to reliably measure the saltation cessation threshold for low-pressure atmospheric conditions

As described in the previous section, a reliable wind tunnel measurement of the cessation threshold for a given low-pressure atmospheric condition requires that equilibrium transport conditions can be established, at least temporarily. Since we are currently unable to predict with confidence the fetch distance saltation needs to reach equilibrium as a function of the atmospheric pressure, and since the required fetch could potentially be very large, it makes sense to design an experimental set-up that allows for adjustable sand feeding. However, this may be challenging given the closed-conduit nature of pressurised-wind tunnels. Once equilibrium transport is established in one way or another, we recommend to carry out measurements in the continuous-transport regime of the equilibrium transport rate $Q$ (or a proxy thereof) as a function of the shear velocity $u_\ast$ and then extrapolate these measurements to vanishing transport using $Q=c_1(u_\ast ^2-u_{\ast t}^2)+c_2(u_\ast ^2-u_{\ast t}^2)^2$ (consistent with (3.2)), where $c_1$, $c_2$ and $u_{\ast t}$ are treated as fit parameters. The resulting cessation threshold $u_{\ast t}$ should be substantially smaller than the dynamic-transport threshold. In fact, for the terrestrial wind tunnel measurements by Creyssels et al. (2009), this extrapolation method yields the value $u_{\ast t}\simeq 0.13\ \mathrm {m}\ \mathrm {s}^{-1}$ (Pähtz & Durán 2020), which is nearly a factor of $2$ smaller than the smallest wind shear velocity ($u_\ast \simeq 0.24\ \mathrm {m}\ \mathrm {s}^{-1}$) for which Creyssels et al. (2009) reported continuous equilibrium transport.

5. Conclusions

Guided by simulations with a well-established DEM-based numerical model (Durán et al. 2012) and existing experimental data, we have semi-empirically derived the scaling behaviours of the cessation threshold shear velocity $u_{\ast t}$ and rate $Q$ of equilibrium sediment transport across almost seven orders of the particle–fluid density ratio $s$, ranging from subaqueous transport ($s\approx 2.65$) to aeolian transport in the highly rarefied atmosphere on Pluto ($s\approx 10^7$). For saltation transport, occurring in planetary aeolian environments, they are

(5.1)\begin{gather} u_{{\ast} t}=2.3(\tilde{g}\nu)^{1/3}\max\left[(D_\ast{/}18)^{{-}1/2},(D_\ast{/}18)^{1/2}\right], \end{gather}

(5.2)\begin{gather}Q=1.7s^{1/3}\rho_p(d/\tilde{g})^{1/2}(u_\ast^2-u_{{\ast} t}^2)+12s^{{-}2/3}\rho_p(\tilde{g}^3d)^{{-}1/2}(u_\ast^2-u_{{\ast} t}^2)^2, \end{gather}

where $\rho _p$ is the particle density, $\nu$ the kinematic fluid viscosity, $\tilde {g}\equiv (1-1/s)g$ the buoyancy-reduced gravity, $d$ the median grain diameter and $D_\ast \equiv \sqrt {s}\tilde {g}d/(\tilde {g}\nu )^{2/3}$ its normalised value. In highly rarefied atmospheres ($s\gtrsim 10^6$), $D_\ast$ is calculated by the more general (3.20), accounting for drag rarefaction effects. Put together, (5.1) and (5.2) constitute a simple means to make predictions of aeolian processes across a large range of planetary conditions.

The derivation of (5.1) consists of a first-principle-based proof of the statement that $u_{\ast t}/(\tilde {g}\nu )^{1/3}$ is a function of only $D_\ast$ (§ 3.2.1). In contrast to existing cessation threshold models, this proof does not resort to coarse-graining the particle phase of the aeolian transport layer above the bed surface, but requires comparably much milder assumptions. Its arguably most critical underlying assumption is that scale-free boundary conditions describe the outcome of grain–bed collisions. The validation of the above statement with our extensive simulation data set therefore indicates that the characteristic velocity scale $\sqrt {\tilde {g}d}$ of grains ejected by the splash of an impacting grain plays no important role for the physics behind the cessation threshold. Instead, grain–bed rebounds and/or splash ejection events associated with the upper-tail end of the ejection velocity distribution are seemingly the physical processes that need to be considered.

The left and right term of the right-hand side of (5.2) are consistent with the saltation limit and collisional limit, respectively, of the $Q$-scaling derived by Pähtz & Durán (2020), with a threshold mean grain velocity scaling as $V_t\propto s^{1/3}\sqrt {\tilde {g}d}$. This scaling strongly deviates from the previous physical transport laws by Ungar & Haff (1987, $V_t\propto \sqrt {\tilde {g}d}$) and Durán et al. (2011, $V_t\propto u_{\ast t}$). For example, the law by Ungar & Haff (1987), which has been adjusted to Earth conditions, underestimates the sediment transport rate for the simulated Martian conditions by a factor of about $5$. Only the recent model of Pähtz et al. (2021) comes somewhat close to reproducing this scaling. It captures the $s^{1/3}$-dependence of $V_t$, but fails to capture its proportionality to $\sqrt {\tilde {g}d}$. This hints at a quite fundamental lack of understanding of the transport rate physics and calls for future studies on this issue.

For Martian atmospheric conditions, the cessation threshold values predicted by the numerical simulations and (5.1) are much smaller than the recent dynamic-threshold measurements by Andreotti et al. (2021). This is particularly odd given that both the numerical simulations and (5.1) are in agreement with terrestrial experimental data. If the simulations were fundamentally wrong, one would expect them to fail for all conditions, not just for Martian conditions. In § 4, we have therefore presented arguments for why the experimental methods used by Andreotti et al. (2021) may have yielded a threshold different from $u_{\ast t}$. This issue needs to be resolved in future studies, since knowing the ‘true’ value of $u_{\ast t}$ is crucial for understanding the time evolution of Martian landscapes.