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. Reference Bourke, Lancaster, Fenton, Parteli, Zimbelman and Radebaugh2010; Kok et al. Reference Kok, Parteli, Michaels and Karam2012; Diniega et al. Reference Diniega, Kreslavsky, Radebaugh, Silvestro, Telfer and Tirsch2017). 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 Reference Sagan and Chyba1990), Pluto (Telfer et al. Reference Telfer2018) and the comet 67P/Churyumov-Gerasimenko (Thomas et al. Reference Thomas2015; Jia, Andreotti & Claudin Reference Jia, Andreotti and Claudin2017).
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 Reference Kok2010a; Durán Vinent et al. Reference Durán Vinent, Andreotti, Claudin and Winter2019). 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. Reference Pähtz, Clark, Valyrakis and Durán2020, § 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 Reference Kok2010a), Pluto (Telfer et al. Reference Telfer2018) and Saturn's moon Titan (Comola et al. Reference Comola, Kok, Lora, Cohanim, Yu, He, McGuiggan, Hörst and Turney2022), as well as in Antarctica. In fact, although Antarctica's surface is covered by very cohesive (Pomeroy & Gray Reference Pomeroy and Gray1990) old snow and ice (cohesion increases the initiation threshold probably much more than the cessation threshold (Comola et al. Reference Comola, Kok, Chamecki and Martin2019b, Reference Comola, Kok, Lora, Cohanim, Yu, He, McGuiggan, Hörst and Turney2022; Pähtz et al. Reference Pähtz, Liu, Xia, Hu, He and Tholen2021; Besnard et al. Reference Besnard, Dupont, Ould El Moctar and Valance2022)), aeolian snow and ice transport occurs there even at relatively low wind speeds that are likely much below the initiation threshold (Leonard et al. Reference Leonard, Tremblay, Thom and MacAyeal2011).
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 Reference Ungar and Haff1987; Andreotti Reference Andreotti2004; Claudin & Andreotti Reference Claudin and Andreotti2006; Kok & Renno Reference Kok and Renno2009; Kok Reference Kok2010b; Durán, Claudin & Andreotti Reference Durán, Claudin and Andreotti2011; Berzi, Jenkins & Valance Reference Berzi, Jenkins and Valance2016; Berzi, Valance & Jenkins Reference Berzi, Valance and Jenkins2017; Lämmel & Kroy Reference Lämmel and Kroy2017; Pähtz & Durán Reference Pähtz and Durán2018a, Reference Pähtz and Durán2020; Andreotti et al. Reference Andreotti, Claudin, Iversen, Merrison and Rasmussen2021; Pähtz et al. Reference Pähtz, Liu, Xia, Hu, He and Tholen2021; Comola et al. Reference Comola, Kok, Lora, Cohanim, Yu, He, McGuiggan, Hörst and Turney2022; Gunn & Jerolmack Reference Gunn and Jerolmack2022). 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. Reference Pähtz, Clark, Valyrakis and Durán2020; Gunn & Jerolmack Reference Gunn and Jerolmack2022).
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. Reference Pähtz, Clark, Valyrakis and Durán2020), 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 Reference Andreotti2004; Claudin & Andreotti Reference Claudin and Andreotti2006; Kok & Renno Reference Kok and Renno2009; Kok Reference Kok2010b; Andreotti et al. Reference Andreotti, Claudin, Iversen, Merrison and Rasmussen2021; Comola et al. Reference Comola, Kok, Lora, Cohanim, Yu, He, McGuiggan, Hörst and Turney2022). On the other hand, it has been modelled as a ‘rebound threshold’ (Pähtz et al. Reference Pähtz, Clark, Valyrakis and Durán2020), 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. Reference Berzi, Valance and Jenkins2017; Pähtz et al. Reference Pähtz, Liu, Xia, Hu, He and Tholen2021; Gunn & Jerolmack Reference Gunn and Jerolmack2022). 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 Reference Pähtz and Durán2018a; Pähtz et al. Reference Pähtz, Clark, Valyrakis and Durán2020; Pähtz et al. Reference Pähtz, Liu, Xia, Hu, He and Tholen2021). 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. (Reference Pähtz, Liu, Xia, Hu, He and Tholen2021) proposed that the recent dynamic-threshold measurements in a low-pressure wind tunnel by Andreotti et al. (Reference Andreotti, Claudin, Iversen, Merrison and Rasmussen2021) 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. Reference Andreotti, Claudin, Iversen, Merrison and Rasmussen2021; Gunn & Jerolmack Reference Gunn and Jerolmack2022).
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 Reference Durán, Andreotti and Claudin2012, 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. (Reference Andreotti, Claudin, Iversen, Merrison and Rasmussen2021), 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. (Reference Durán, Andreotti and Claudin2012), 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
\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 Reference Camenen2007), 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. Reference Crowe, Schwarzkopf, Sommerfeld and Tsuji2012). 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:
\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 Reference Pähtz and Durán2018a,Reference Pähtz and Duránb). 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 (Reference Pähtz and Durán2018a) and figure 4 of Pähtz & Durán (Reference Pähtz and Durán2020)), height profiles of relevant equilibrium transport properties (figure 2 of Pähtz & Durán (Reference Pähtz and Durán2018a) and figure 6 of Durán, Andreotti & Claudin (Reference Durán, Andreotti and Claudin2014a)) and aeolian ripple formation (Durán, Claudin & Andreotti Reference Durán, Claudin and Andreotti2014b).
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 Reference Pähtz and Durán2018b)
\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 Reference Pähtz and Durán2018a)
\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 Reference Pähtz and Durán2018b), 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 Reference Pähtz and Durán2018b)
\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 Reference Pähtz and Durán2018b)
\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,
\begin{gather} d^+=\tilde{g}d/(\tilde{g}\nu)^{2/3}, \end{gather}
\begin{gather}u^+_\ast=u_\ast{/}(\tilde{g}\nu)^{1/3}, \end{gather}
\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 Reference Pähtz and Durán2020):
\begin{gather} s\equiv\rho_p/\rho_f=1/\rho_f^{+}, \end{gather}
\begin{gather}Ga\equiv\sqrt{s\tilde{g}d^3}/\nu=\sqrt{s}d^{{+}3/2}, \end{gather}
\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 Reference Pähtz and Durán2018a, Reference Pähtz and Durán2020).
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 Reference Pähtz and Durán2018a)
\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. Reference Greeley, White, Leach, Iversen and Pollack1976, Reference Greeley, Leach, White, Iversen and Pollack1980; Iversen & White Reference Iversen and White1982; Greeley et al. Reference Greeley, Iversen, Leach, Marshall, White and Williams1984; Burr et al. Reference Burr, Bridges, Marshall, Smith, White and Emery2015, Reference Burr, Sutton, Emery, Nield, Kok, Smith and Bridges2020; Swann et al. Reference Swann, Sherman and Ewing2020) (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. Reference Pähtz, Clark, Valyrakis and Durán2020). Roughly the same 
$s^{1/2}$-scaling was also found for the dynamic-threshold measurements by Andreotti et al. (Reference Andreotti, Claudin, Iversen, Merrison and Rasmussen2021) 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. Reference Greeley, White, Leach, Iversen and Pollack1976, Reference Greeley, Leach, White, Iversen and Pollack1980; Iversen & White Reference Iversen and White1982; Greeley et al. Reference Greeley, Iversen, Leach, Marshall, White and Williams1984; Burr et al. Reference Burr, Bridges, Marshall, Smith, White and Emery2015, Reference Burr, Sutton, Emery, Nield, Kok, Smith and Bridges2020; Swann, Sherman & Ewing Reference Swann, Sherman and Ewing2020) and cessation threshold measurements for aeolian transport of quartz (Bagnold Reference Bagnold1937; Martin & Kok Reference Martin and Kok2018; Zhu et al. Reference Zhu, Huo, Zhang, Wang, Pähtz, Huang and He2019), clay loam (Chepil Reference Chepil1945) and snow at sea level (Sugiura et al. Reference Sugiura, Nishimura, Maeno and Kimura1998) and high altitude (Clifton, Rüedi & Lehning Reference Clifton, Rüedi and Lehning2006, HA). The dynamic-threshold measurements by Andreotti et al. (Reference Andreotti, Claudin, Iversen, Merrison and Rasmussen2021) 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$:
\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. (Reference Creyssels, Dupont, Ould El Moctar, Valance, Cantat, Jenkins, Pasini and Rasmussen2009) 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 (Reference Ungar and Haff1987), whereas Pähtz & Durán (Reference Pähtz and Durán2020) 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 Reference Bagnold1937; Martin & Kok Reference Martin and Kok2018; Zhu et al. Reference Zhu, Huo, Zhang, Wang, Pähtz, Huang and He2019), clay loam (Chepil Reference Chepil1945) and snow at sea level (Sugiura et al. Reference Sugiura, Nishimura, Maeno and Kimura1998) and high altitude (Clifton et al. Reference Clifton, Rüedi and Lehning2006, 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)
\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. Reference Berzi, Jenkins and Valance2016), encoding the importance of grain inertia relative to viscous drag forcing, and controls the transition to viscous bedload (Pähtz et al. Reference Pähtz, Liu, Xia, Hu, He and Tholen2021).
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. Reference Creyssels, Dupont, Ould El Moctar, Valance, Cantat, Jenkins, Pasini and Rasmussen2009; Ho et al. Reference Ho, Valance, Dupont and Ould El Moctar2011; Ho Reference Ho2012; Martin & Kok Reference Martin and Kok2017; Ralaiarisoa et al. Reference Ralaiarisoa, Besnard, Furieri, Dupont, Ould El Moctar, Naaim-Bouvet and Valance2020) and snow (Sugiura et al. Reference Sugiura, Nishimura, Maeno and Kimura1998). 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. (Reference Ralaiarisoa, Besnard, Furieri, Dupont, Ould El Moctar, Naaim-Bouvet and Valance2020) 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 Reference Claudin and Andreotti2006). These parameters involve three units (mass, length and time). According to the 
$\varPi$ theorem (Barenblatt Reference Barenblatt1996), 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$:
\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 Reference Claudin and Andreotti2006; Kok Reference Kok2010b; Berzi et al. Reference Berzi, Jenkins and Valance2016, Reference Berzi, Valance and Jenkins2017; Pähtz & Durán Reference Pähtz and Durán2018a; Andreotti et al. Reference Andreotti, Claudin, Iversen, Merrison and Rasmussen2021; Pähtz et al. Reference Pähtz, Liu, Xia, Hu, He and Tholen2021; Gunn & Jerolmack Reference Gunn and Jerolmack2022). 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 Reference Claudin and Andreotti2006; Kok Reference Kok2010b; Berzi et al. Reference Berzi, Jenkins and Valance2016, Reference Berzi, Valance and Jenkins2017; Pähtz & Durán Reference Pähtz and Durán2018a; Andreotti et al. Reference Andreotti, Claudin, Iversen, Merrison and Rasmussen2021; Pähtz et al. Reference Pähtz, Liu, Xia, Hu, He and Tholen2021; Gunn & Jerolmack Reference Gunn and Jerolmack2022) or by an average motion behaviour (Kok Reference Kok2010b; Pähtz & Durán Reference Pähtz and Durán2018a) that is mathematically equivalent to an identical periodic trajectory representation (Pähtz et al. Reference Pähtz, Clark, Valyrakis and Durán2020).
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.
(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.
(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).(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.
(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.(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. Reference Beladjine, Ammi, Oger and Valance2007), 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. Reference Lämmel, Dzikowski, Kroy, Oger and Valance2017). 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. Reference Pähtz, Liu, Xia, Hu, He and Tholen2021)
\begin{gather} \dot{v}^+_z={-}1-v^+_z/v^+_s, \end{gather}
\begin{gather}\dot{v}^+_x=(u^+_x-v^+_x)/v^+_s, \end{gather}
\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. Reference Pähtz, Liu, Xia, Hu, He and Tholen2021). 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}$:
\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):
\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
\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)
\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 Reference Shao and Lu2000), 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 Reference Ungar and Haff1987; Durán et al. Reference Durán, Claudin and Andreotti2011; Kok et al. Reference Kok, Parteli, Michaels and Karam2012; Berzi et al. Reference Berzi, Jenkins and Valance2016). However, it has been shown that, for sufficiently intense saltation, midair collisions significantly disturb grain trajectories (Carneiro et al. Reference Carneiro, Araújo, Pähtz and Herrmann2013; Pähtz & Durán Reference Pähtz and Durán2020; Ralaiarisoa et al. Reference Ralaiarisoa, Besnard, Furieri, Dupont, Ould El Moctar, Naaim-Bouvet and Valance2020), leading to an additional additive term increasing as 
$M^+/d^+$ (Pähtz & Durán Reference Pähtz and Durán2020):
\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. Reference Durán, Andreotti and Claudin2012; Pähtz & Durán Reference Pähtz and Durán2018b). 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 Reference Pähtz and Durán2018b):
\begin{gather} M\equiv\int_0^\infty\rho\,\mathrm{d} z=M_0, \end{gather}
\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 Reference Pähtz and Durán2018b). 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. Reference Durán, Andreotti and Claudin2012):
\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}
\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)):
\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}
\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 Reference Ungar and Haff1987; Durán et al. Reference Durán, Claudin and Andreotti2011; Berzi et al. Reference Berzi, Jenkins and Valance2016; Pähtz & Durán Reference Pähtz and Durán2020) 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 Reference Ungar and Haff1987; Berzi et al. Reference Berzi, Jenkins and Valance2016) and 
$V^+_t\propto u^+_{\ast t}$ (Durán et al. Reference Durán, Claudin and Andreotti2011; Kok et al. Reference Kok, Parteli, Michaels and Karam2012; Pähtz & Durán Reference Pähtz and Durán2020).
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. (Reference Pähtz, Liu, Xia, Hu, He and Tholen2021), which we here compare with the most recent alternative, that of Gunn & Jerolmack (Reference Gunn and Jerolmack2022). 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. (Reference Pähtz, Liu, Xia, Hu, He and Tholen2021) captures the simulation data very well, the model of Gunn & Jerolmack (Reference Gunn and Jerolmack2022), 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. (Reference Pähtz, Liu, Xia, Hu, He and Tholen2021) and (b) Gunn & Jerolmack (Reference Gunn and Jerolmack2022), 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 (Reference Ungar and Haff1987) and others (e.g. Jenkins & Valance Reference Jenkins and Valance2014; Berzi et al. Reference Berzi, Jenkins and Valance2016), 
$Q^+/d^{+3/2}=f_1(\varTheta -\varTheta _t)$ (figure 8a) and the model of Durán et al. (Reference Durán, Claudin and Andreotti2011) and others (Kok et al. Reference Kok, Parteli, Michaels and Karam2012; Pähtz & Durán Reference Pähtz and Durán2020), 
$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 (Reference Ungar and Haff1987) and Durán et al. (Reference Durán, Claudin and Andreotti2011). 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 Reference Pähtz and Durán2020). 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. Reference Pähtz, Liu, Xia, Hu, He and Tholen2021), 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
\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 Reference Bagnold1937; Martin & Kok Reference Martin and Kok2018; Zhu et al. Reference Zhu, Huo, Zhang, Wang, Pähtz, Huang and He2019), clay loam (Chepil Reference Chepil1945) and snow at sea level (Sugiura et al. Reference Sugiura, Nishimura, Maeno and Kimura1998) and high altitude (Clifton et al. Reference Clifton, Rüedi and Lehning2006, HA), and a compilation of experimental threshold data for subaqueous bedload (Buffington & Montgomery Reference Buffington and Montgomery1997). 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
\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. Reference Creyssels, Dupont, Ould El Moctar, Valance, Cantat, Jenkins, Pasini and Rasmussen2009; Ho et al. Reference Ho, Valance, Dupont and Ould El Moctar2011; Ho Reference Ho2012; Martin & Kok Reference Martin and Kok2017; Ralaiarisoa et al. Reference Ralaiarisoa, Besnard, Furieri, Dupont, Ould El Moctar, Naaim-Bouvet and Valance2020) and snow (Sugiura et al. Reference Sugiura, Nishimura, Maeno and Kimura1998), and subaqueous bedload (Meyer-Peter & Müller Reference Meyer-Peter and Müller1948; Smart & Jaeggi Reference Smart and Jaeggi1983; Capart & Fraccarollo Reference Capart and Fraccarollo2011). We corrected the raw laboratory data by Smart & Jaeggi (Reference Smart and Jaeggi1983) and Capart & Fraccarollo (Reference Capart and Fraccarollo2011) for sidewall drag using the method described in § 2.3 of Guo (Reference Guo2015) 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. Reference Pähtz, Liu, Xia, Hu, He and Tholen2021). 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. (Reference Ralaiarisoa, Besnard, Furieri, Dupont, Ould El Moctar, Naaim-Bouvet and Valance2020) 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. (Reference Pähtz, Liu, Xia, Hu, He and Tholen2021) (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. (Reference Pähtz, Liu, Xia, Hu, He and Tholen2021). 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. Reference Crowe, Schwarzkopf, Sommerfeld and Tsuji2012), 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. Reference Crowe, Schwarzkopf, Sommerfeld and Tsuji2012) 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$):
\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 Reference Martin and Kok2018). 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 Reference Martin and Kok2018). 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. Reference Clifton, Rüedi and Lehning2006; Zhu et al. Reference Zhu, Huo, Zhang, Wang, Pähtz, Huang and He2019) or directly measure the cessation of intermittent saltation (Bagnold Reference Bagnold1937; Chepil Reference Chepil1945; Sugiura et al. Reference Sugiura, Nishimura, Maeno and Kimura1998; Martin & Kok Reference Martin and Kok2018; Zhu et al. Reference Zhu, Huo, Zhang, Wang, Pähtz, Huang and He2019). Both methods require that equilibrium transport conditions can at least temporarily be established during the experiments (Pähtz et al. Reference Pähtz, Clark, Valyrakis and Durán2020), 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. (Reference Clifton, Rüedi and Lehning2006), 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. Reference Comola, Gaume, Kok and Lehning2019a). 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 Reference Claudin and Andreotti2006; Kok Reference Kok2010b; Berzi et al. Reference Berzi, Valance and Jenkins2017; Andreotti et al. Reference Andreotti, Claudin, Iversen, Merrison and Rasmussen2021; Gunn & Jerolmack Reference Gunn and Jerolmack2022). 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 Reference Pähtz and Durán2018a; Pähtz et al. Reference Pähtz, Clark, Valyrakis and Durán2020; Pähtz et al. Reference Pähtz, Liu, Xia, Hu, He and Tholen2021). An important example for a potential misinterpretation of measured dynamic thresholds as cessation thresholds is the study by Andreotti et al. (Reference Andreotti, Claudin, Iversen, Merrison and Rasmussen2021) for the following reasons.
(i) Andreotti et al. (Reference Andreotti, Claudin, Iversen, Merrison and Rasmussen2021), 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’.(ii) Andreotti et al. (Reference Andreotti, Claudin, Iversen, Merrison and Rasmussen2021) 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 Reference Martin and Kok2018). For low-pressure conditions, the measurements are difficult to interpret due to the lack of equilibrium transport.
(iii) Andreotti et al. (Reference Andreotti, Claudin, Iversen, Merrison and Rasmussen2021) 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.
(iv) Pähtz et al. (Reference Pähtz, Liu, Xia, Hu, He and Tholen2021) 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. (Reference Andreotti, Claudin, Iversen, Merrison and Rasmussen2021), 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. (Reference Andreotti, Claudin, Iversen, Merrison and Rasmussen2021) as cessation threshold measurements is highly relevant, since it led Gunn & Jerolmack (Reference Gunn and Jerolmack2022) 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. Reference Berzi, Jenkins and Valance2016, Reference Berzi, Valance and Jenkins2017; Andreotti et al. Reference Andreotti, Claudin, Iversen, Merrison and Rasmussen2021), which is experimentally known to control the viscous damping of frontal grain collisions with a flat plate (Gondret, Lance & Petit Reference Gondret, Lance and Petit2002). 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 Reference Pähtz and Durán2018a). Even for 
$e=0.01$, grains can saltate in large hops on the order of 
$100d$ high (Movie S3 of Pähtz & Durán Reference Pähtz and Durán2018a).
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. (Reference Creyssels, Dupont, Ould El Moctar, Valance, Cantat, Jenkins, Pasini and Rasmussen2009), this extrapolation method yields the value 
$u_{\ast t}\simeq 0.13\ \mathrm {m}\ \mathrm {s}^{-1}$ (Pähtz & Durán Reference Pähtz and Durán2020), 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. (Reference Creyssels, Dupont, Ould El Moctar, Valance, Cantat, Jenkins, Pasini and Rasmussen2009) reported continuous equilibrium transport.
5. Conclusions
Guided by simulations with a well-established DEM-based numerical model (Durán et al. Reference Durán, Andreotti and Claudin2012) 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
\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}
\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 (Reference Pähtz and Durán2020), 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 (Reference Ungar and Haff1987, 
$V_t\propto \sqrt {\tilde {g}d}$) and Durán et al. (Reference Durán, Claudin and Andreotti2011, 
$V_t\propto u_{\ast t}$). For example, the law by Ungar & Haff (Reference Ungar and Haff1987), 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. (Reference Pähtz, Liu, Xia, Hu, He and Tholen2021) 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. (Reference Andreotti, Claudin, Iversen, Merrison and Rasmussen2021). 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. (Reference Andreotti, Claudin, Iversen, Merrison and Rasmussen2021) 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.
Funding
T.P. acknowledges support from the National Natural Science Foundation of China (no. 12272344). O.D. acknowledges support from the Texas A&M Engineering Experiment Station.
Declaration of interests
The authors report no conflict of interest.
