A two-species continuum model for aeolian sand transport

Wind-driven sand transport is the most noticeable process shaping the morphology of arid regions on the Earth, Mars and elsewhere. It is responsible for the spontaneous creation of a whole hierarchy of self-organized dynamic structures from ripples over isolated dunes to devastating fields of shifting sands. It also contributes considerably to dust proliferation, which is a major determinant of our global climate. There is thus an urgent need for mathematical models that can efficiently and accurately predict aeolian sand fluxes. However, the task is very much complicated by the complex turbulent flow of the driving medium (e.g. streaming air) and the erratic nature of the grain hopping it excites [1]. Yet, it is by now well understood that this hopping of grains accelerated by the wind has some characteristic structure [2]. Highly energetic 'saltating' grains make a dominant contribution to the overall mass transport. When impacting on the sand bed, they dissipate some of their energy in a complex process called splash, ejecting a cloud of 'reptating' grains [3, 4]. A snapshot of this splash would show a whole ensemble of trajectories corresponding to a distribution of jump lengths from short, over intermediate, to large hops. However, grain-scale studies and theoretical considerations have indicated that a reduced description in terms of only two idealized populations (sometimes called saltons and reptons) should indeed be able to provide a faithful parametrization of the complex aeolian sand transport process and the ensuing structure formation [2, 5]. A rule of thumb to say which grains in a real splash pertain to idealized population of saltating or reptating grains is whether or not their jump height exceeds a threshold of about a few grain diameters.

A drawback of the two-species description has been that it is still conceptually and computationally quite demanding. For reasons of simplicity and computational efficiency, many theoretical studies have therefore chosen to reduce the mathematical description even further, to mean-field-type 'single-trajectory' models [6]. This may be justified for certain purposes, e.g. for the mathematical modelling of aeolian sand dunes, which are orders of magnitude larger than the characteristic length scales involved in the saltation process and thus not expected to be very sensitive to the details on the grain scale. Their formation and migration is thought to predominantly depend on large-scale features of the wind field and of the saltation flux, chiefly the symmetry breaking of the turbulent flow over the dune and the delayed reaction of the sand transport to the wind [7]. Moreover, the reptating grains, although they are many, are generally thought to contribute less to the overall sand transport, because they have short trajectories and quickly get trapped in the bed again. Therefore, it seems admissible to concentrate on the saltating particles. On this basis, numerically efficient models for one effective grain species that can largely be identified with the saltating grain fraction have been constructed. The Sauermann model [8] is a popular and widely used example of such mean-field continuum models. One-species models have become a very successful means of gaining analytical insight into [9–14] and to conduct efficient large-scale numerical simulations of [15–19] the complex structure formation processes caused by aeolian transport. The reduction to a single representative trajectory makes the one-species models analytically tractable and computationally efficient. However, it is also responsible for some weaknesses concerning both the way in which the two species are subsumed into one [5] and how they feed back onto the wind [12]. These entail imperfections in the model predictions, most noticeably a systematic overestimation of the stationary flux at high wind speed (see figure 5). Therefore, the one-species models have been criticized for their lack of numerical accuracy and internal consistency [5, 20]. There is also a number of interesting phenomena that cannot be quantitatively modelled within a one-species model, because they specifically depend on one of the two species, or on their interaction, or because the two species are not merely conceptually but also physically different.

As three important representatives of such phenomena, we name dust emission from the desert, ripples and megaripples. Dust is created, exposed to the wind and emitted by aeolian sand transport, and saltating and reptating particles play quite different roles in these processes [1, 21]. For example, fragmentation processes might be driven by the bombardment of high-energy saltating grains [22], but not by the slower reptating grains. In contrast, the emission of dust hidden from the wind underneath the top layer of grains in a sand bed could arguably be linked to the absolute number of reptating grains dislodged from the bed. The reptating grain fraction and its sensitivity to the local slope of the sand bed are, moreover, held responsible for the spontaneous evolution of aeolian ripples [1, 2, 4, 23–29]. But the number of reptating particles depends on the strength of the splash caused by the saltating grains. This interdependence of the reptating and saltating species becomes even more transparent when the two species correspond to visibly different types of grains. This is the case for megaripples, which have wavelengths in the metre range and form on strongly polydisperse sand beds in a process accompanied by a pronounced grain sorting [30–32]. The mechanism underlying their formation and evolution is still not entirely understood, but it seems that the highly energetic bombardment of small saltating grains drives the creep of the larger grains armouring the ripple crests. The mentioned grain-sorting emerges due to the grain-size-dependent hop lengths.

The demand for more faithful descriptions of aeolian sand transport has recently spurred notable efforts to eliminate the deficiencies of the one-species models [11, 12, 33] or to develop a numerically efficient two-species model [5]. It also motivated the development of the model described in detail below, as we felt that the problem has not yet been cured at its root. In our opinion, many of the amendments proposed so far have targeted the symptoms rather than the disease by invoking some ad hoc assumptions. It was our aim to modify the sand transport equations from bottom-up, based on a careful analysis of the physics on the grain scale and including the feedback of the two dissimilar grain populations on the turbulent flow. In this way, we derived a conceptually simple, analytically tractable and numerically efficient two-species model with only two phenomenological fit parameters. They serve to parametrize the complicated splash process and take very reasonable values if the predicted flux law is fitted to measured data. Compared with the time when the first continuum sand transport models were formulated, we could build on more detailed experimental and theoretical knowledge of the grain-scale physics [5, 34–39] and rely on more comprehensive empirical information about the wind shear stress and sand transport to test our predictions [40–43].

The plan of the paper is as follows. We first summarize some of the pertinent basic notions of turbulent flows and aeolian sand transport and introduce the two-species formalism. Then we implement the two-species physics on the level of the sand flux in section 3.1 and on the level of the turbulent closure for the wind in section 3.2. In sections 3.3 and 3.4, we address the most interesting observables, namely the speed and the average number of hopping grains, and the saturated sand flux. We finally compare our results with other models and experiments in section 4, before concluding in section 5.

Our analysis of the two-phase flow of air and sand is similar, in spirit, to the Sauermann model [8]. The Sauermann model is a continuum or hydrodynamic model. Rather than with the positions and velocities of individual grains it deals with spatio-temporal fields of these quantities, namely the local mass density $\varrho (\vec {x},t)$ and sand transport velocity $\vec {v}(\vec {x},t)$. It is of mean-field type since it reduces both fields to the mean quantities ρ(x,t), v(x,t) for a single representative trajectory, characterizing the conditions along the wind direction (the x-direction), with the distribution of airborne grains in the vertical direction (the z-direction) integrated out. Note that ρ therefore is an area density, obtained by integrating the volume density ϱ over z. (The vertical grain distribution is a crucial element of the modified turbulent closure, however, which accounts for the feedback of the sand onto the wind.) In this paper, we deal exclusively with the saturated sand flux, i.e. the sand flux along the wind direction over a flat sand bed and under stationary wind conditions. What we call the flux q = ρv is actually a vertically integrated flux, or a sand transport rate [8]. The central aim is to derive a constitutive relation q(τ), giving the stationary flux at a given shear stress τ.

In the following, the strong mean-field approximation of the one-species models is somewhat relaxed, and reptating, and saltating particles with densities ρ^rep, ρ^sal and transport velocities v^rep, v^sal are treated separately. The flux is split up accordingly

$$\begin{equation} q = q^{\mathrm{rep}} + q^{\mathrm{sal}} = \rho^{\mathrm{sal}} v^{\mathrm{sal}} + \rho^{\mathrm{rep}} v^{\mathrm{rep}}. \end{equation} \tag{ 1 }$$

Introducing the dimensionless mass fraction φ∈[0,1] of saltating grains relative to the total number of mobile grains, the (integrated) densities of the mobile grains can be written in the form

$$\begin{equation} \rho^{\mathrm{sal}} = \varphi \rho\quad {\rm and} \quad \rho^{\mathrm{rep}} = (1-\varphi) \rho. \end{equation} \tag{ 2 }$$

Accordingly, the flux balance and the relative contribution of reptating and saltating grains to the flux read as

$$\begin{equation} q = \rho v= \rho \left[ \varphi v^{\mathrm{sal}} + (1-\varphi) v^{\mathrm{rep}} \right] \quad \mathrm{and} \quad q^{\mathrm{rep}}/q^{\mathrm{sal}} = (\varphi^{-1}-1) v^{\mathrm{rep}}/v^{\mathrm{sal}} . \end{equation} \tag{ 3 }$$

The fluxes depend on the wind velocity field, which, in turn, depends on the grain–wind feedback. To make progress in this respect, the (steady) shear stress τ in the saltation layer is split into two components carried by the airborne grains ('g') and the air ('a') itself, respectively,

$$\begin{equation} \tau = \tau_{\mathrm{g}}(z)+\tau_{\mathrm{a}}(z) = {\rm const.} \end{equation} \tag{ 4 }$$

This idea, introduced by Owen [44], which amounts to an effective two-fluid representation of the wind and the airborne sand, was used in many analytical models [5, 8, 11, 45] and confirmed by numerical simulations [6, 37–39, 46]. In the two-species model, we moreover keep track of the two species of mobile grains in (4). Similar to the flux division introduced in (1), we further split the grain-borne shear stress

$$\begin{equation} \tau_{\mathrm{g}}(z) = \tau_{\mathrm{g}}^{\mathrm{rep}}(z) + \tau_{\mathrm{g}}^{\mathrm{sal}}(z). \end{equation} \tag{ 5 }$$

Concerning the feedback of the grains on the wind, it turns out that the wind velocity profile u(z) is predominantly determined by the reptating grain fraction and thus by the functional form of τ^rep_g(z), while the number and trajectories of the saltating grains hardly affect the wind profile at all. Physically, this makes sense, since the saltating particles form a very disperse gas compared with the dense reptating layer. Also note that the reptating grains, which lose all their momentum upon impact, are responsible for the momentum loss of the flow field. The wind profile may thus be computed in a reduced one-species scheme, as explained in section 3.2 and appendix B.

Under steady conditions, a common simplifying assumption is that the airborne stress near the ground (z = 0) can be identified with the threshold shear stress τ_t required to mobilize grains from the ground [44],

$$\begin{equation} \tau_{\mathrm{a}}(0) \approx \tau_{\mathrm{t}}. \end{equation} \tag{ 6 }$$

If it were larger or smaller, an increasing or decreasing number of grains would be mobilized, respectively, resulting in an unsteady flow. This idea is indeed supported by empirical observations, which find that the number of ejected grains increases (almost linearly) with the impact speed, not their ejection velocity or jump height [34, 35, 47, 48]. The feedback of the grains on the wind thus essentially fixes the air shear stress at the ground (and, if τ ≈ τ_t, also above) to the threshold value τ_t. To be precise, τ_t is called the impact threshold, to emphasize that it is easier to lift grains from the splash cloud than to lift completely immobile grains from the ground [1]. But, for the sake of simplicity, we do not bother to make this distinction here, nor do we speculate about possible small deviations of τ_a(0) from τ_t under steady conditions [12, 33]. Combining (4) and (6), we may express the shear stress contributed by the grains near the ground as

$$\begin{equation} \tau_{\mathrm{g}}(0) = \tau_{\mathrm{g}}^{\mathrm{rep}}(0) + \tau_{\mathrm{g}}^{\mathrm{sal}}(0) =\tau-\tau_{\mathrm{t}}. \end{equation} \tag{ 7 }$$

Both stress contributions are defined as the product of the vertical sand flux ϕⁱ and the horizontal velocity difference vⁱ − vⁱ₀ upon impact, τⁱ_g(0) = (vⁱ − vⁱ₀)ϕⁱ. The vertical sand flux ϕⁱ is given by the horizontal (height integrated) flux qⁱ divided by the grains' mean hop length, ϕⁱ∝qⁱ = ρⁱvⁱ. Here, i represents the species indices 'rep' or 'sal'. Following Sauermann et al [8], we assume parabolic grain trajectories with initial horizontal and vertical velocity components vⁱ₀ and vⁱ_z0, respectively, and obtain

$$\begin{equation} \phi^i = g \rho^i /(2 v_{z0}^i) , \quad \tau_{\mathrm{g}}^i(0) = g \rho^i (v^i-v_0^i)/(2 v_{z0}^i) , \quad i=\mathrm{rep},\,\mathrm{sal}. \end{equation} \tag{ 8 }$$

Since the average ejection angle of the reptating particles is known to be about 80° [34], v^rep₀ ≈ 0 and we identify the total ejection speed of the reptating grains with its vertical component v^rep_z0. The impact angle of saltating grains is known to be almost independent of the wind strength and to take typical mean values of about 10° [49], so that we can identify the total impact speed with the horizontal speed v^sal of the impacting particle, see figure 1.

Zoom In Zoom Out Reset image size

In order to derive the sought-after constitutive equation, we have to understand the balance of the two species of mobile grains. If one considers the vertical fluxes ϕ^sal and ϕ^rep rather than the horizontal fluxes q^sal and q^rep, one has

$$\begin{equation} \phi^{\mathrm{rep}} = N \phi^{\mathrm{sal}}, \end{equation} \tag{ 9 }$$

with the average number N of reptating grains ejected by an impacting saltating grain. Together with the first relation in (8), we then have

$$\begin{equation} \varphi^{-1} = 1 + N v^{\mathrm{rep}}_{z0}/v^{\mathrm{sal}}_{z0}. \end{equation} \tag{ 10 }$$

Using also the second relation in (8) and remembering that we dismissed the horizontal component of the ejection velocity, we can thus write the overall density of mobile grains in the form

$$\begin{equation} \rho = \frac{2 \varrho_{\mathrm{air}}}{g} \frac{u_\ast^2-u_{\ast \mathrm{t}}^2}{\varphi (v^{\mathrm{sal}}-v^{\mathrm{sal}}_0)/v^{\mathrm{sal}}_{z0} + (1-\varphi) v^{\mathrm{rep}}/v^{\mathrm{rep}}_{z0}}. \end{equation} \noindent \tag{ 11 }$$

As often found in the literature, we have employed the notion of the shear velocity or friction velocity u_*, defined by τ ≡ ϱ_airu²_*, as a more intuitive measure of the shear stress τ here. The velocity ratios in the denominator of (11) are supposedly only very weakly dependent on the wind strength. The first one is recognized as an effective restitution coefficient

$$\begin{equation} \alpha \equiv v^{\mathrm{sal}}_{z0}/(v^{\mathrm{sal}}-v^{\mathrm{sal}}_0) \end{equation} \noindent \tag{ 12 }$$

of the saltating grains [8]. Because of the complexity of the splash process, which makes first-principle quantitative estimates of α forbiddingly complex, we suggest treating it as a phenomenological fit parameter. We do, however, provide a reasonable theoretical estimate of the ratio v^rep/v^rep_z0, based on the trajectory of a vertically ejected grain driven by the wind, in section 3.3 and appendix C. Anticipating the result v^rep/v^rep_z0 = 0.7(σν_airg)^1/3/u_*t, we can rewrite the density of the transported sand in the more explicit form

$$\begin{equation} \rho = \frac{2 \varrho_{\mathrm{air}}}{g} \frac{u_\ast^2-u_{\ast \mathrm{t}}^2}{\varphi/\alpha + 0.7 (1-\varphi) (\sigma \nu_{\mathrm{air}} g)^{1/3}/u_{\ast \mathrm{t}}}. \end{equation} \noindent \tag{ 13 }$$

Here ν_air denotes the kinematic viscosity of air, and since the grain–air density ratio σ ≫ 1 under terrestrial conditions, we do not bother to distinguish between σ and σ − 1 here and in the following. Note that on top of the explicit wind strength dependence of ρ via the numerator, there is an implicit, yet undetermined one via φ.

In the following, we develop the sought-after 'second generation' transport law q(u_*), based on the grain-scale physics. Some empirical input is employed to fix certain phenomenological coefficients summarized in table D.1. Our result (see table 1) accounts for essential elements of the aeolian sand transport process beyond the single-trajectory approximation and turns out to be in remarkable agreement with wind tunnel measurements.

Table 1. Scaling functions f(U) for stationary aeolian sand transport relations of the form Q(U) = (1 − U⁻²)f(U). Note that the coefficients of the various models cannot be identified even if represented by the same symbol. Explicit analytical expressions for the coefficients α₀, β₀, γ₀, a and b of the two-species model can be found in appendix D.

This work	$\displaystyle\frac {\alpha _0U + \beta _0 - \gamma _0 U^{-1}}{a U - b}$
This work, one-species limit	a−bU−1  (b<0)
Sauermann et al [8]	$a \sqrt {1+ \alpha _0 U^{-2}} + b U^{-1}$
Sørensen [11]	a+bU−1+cU−2
Durán and Herrmann [12]	a−bU−1  (a,b<0)
Pähtz et al [33]	−a+bU−1+cU−2

The two-species framework is implemented in three steps. First, we deal with the mass balance between the species, then with their feedback onto the wind and their resulting transport velocities, before we finally arrive at an improved stationary transport law q(u_*).

3.1. Sand density: two-species mass balance

We first estimate the partition of the grain population into a saltating and a reptating fraction from several evident assumptions for the splash process, based on fundamental physical principles. In our formalism, a saltating grain always rebounds upon impact and ejects N reptating grains, while a certain fraction of its kinetic energy and momentum is dissipated in the sand bed. We further assume that every reptating grain performs only a single hop after which it remains trapped in the sand bed. These mild approximations for the stationary sand transport reduce the number of free parameters of our model compared with more elaborate descriptions [5].

As argued, e.g., by Kok and Renno [37], the relevant constraint limiting the number N of ejected grains per impact is given by the conservation of the grain-borne momentum rather than the energy. In writing the momentum balance, we make a physically plausible linear-response approximation, namely that all terms scale linearly in the velocity v^sal of the impacting grain. This amounts to assuming a wind-strength-independent scattering geometry, so that we can add z- and x-components up to constant geometric factors that are later absorbed in phenomenological coefficients. Including also the bed losses proportional to v_bed in this vein, we have

$$\begin{equation} v^{\mathrm{sal}} = v^{\mathrm{sal}}_0 + v^{\mathrm{sal}}_{z0} + N v^{\mathrm{rep}}_{z0} +v_{\mathrm{bed}}. \end{equation} \tag{ 14 }$$

For the first two terms on the right-hand side of (14), the proportionality to v^sal is already implicit in the definition (12) of the restitution coefficient α. There is an important subtlety concerning the last two terms, however. Namely, the ejected grains have to gain enough momentum to jump over neighbouring grains in the bed in order to be counted as mobile grains. But they also should not themselves eject other grains upon impact; otherwise we would have to count them among the saltating fraction. Hence, the impact velocity v^rep and the ejection velocity v^rep_z0 of the reptating grains are tightly constrained and cannot, by definition, be strongly dependent on the impact velocity of the saltating grain or the wind strength. Within experimental errors, observations of the collision process of saltating grains by Rice et al [49] are indeed in reasonable agreement with a constant v^rep_z0 of the order of u_*t. In fact, they found the horizontal component of the ejection velocity to be independent of v^sal, and only a slight increase of the vertical component with v^sal. A compelling confirmation of this observation was recently obtained by better controlled laboratory experiments, in which PVC beads were injected at an impact angle of 10° onto a quiescent bed of such beads [34, 35]. There are two important consequences of the constraints on v^rep_z0. Firstly, N should grow linearly with v^sal. Secondly, since the momentum of the impacting grain is only partially transferred to the ejected grains, a critical impact velocity v^sal = v^sal_c ≫ v^rep_z0 has to be overcome to mobilize any grains at all. It can be interpreted as a constant offset in v_bed, which does not scale with v^sal. According to the collision experiments [34, 35], one can take $v^{\mathrm{sal}}_{\mathrm{c}}\approx 40\sqrt {g d}\approx 9 u_{\ast \mathrm{t}}$, where we used Bagnold's estimate u²_*t ≈ 0.01σgd [1]. A concise summary of the empirical input entering our model is provided in table D.1.

Summarizing the above discussion, we can rewrite the momentum balance (14) as

$$\begin{equation} N v^{\mathrm{rep}}_{z0} \propto v^{\mathrm{sal}} - v^{\mathrm{sal}}_{\mathrm{c}}, \end{equation} \tag{ 15 }$$

where v^sal > v^sal_c is tacitly assumed to hold throughout our discussion, and the omitted factor of proportionality should be insensitive to the wind strength. Observations of saltating grains [49], model collision experiments [34, 35, 47, 48], particle dynamics simulations [6, 50] and reduced numerical models, such as the binary collision scheme proposed by Valance and Crassous [36], all support this affine relation of the number of ejected grains, which is the cornerstone of our two-species mass balance. To fix the omitted factor in (15) and make contact with (10), we write

$$\begin{equation} N v^{\mathrm{rep}}_{z0} = (\eta/ \alpha) ( 1- v^{\mathrm{sal}}_{\mathrm{c}}/v^{\mathrm{sal}})v^{\mathrm{sal}}_{z0}. \end{equation} \noindent \tag{ 16 }$$

The constant η, which serves as the second free fit parameter of our model, determines together with v^sal_c how the momentum lost by a rebounding grain is distributed between the bed and the ejected particles. (Details of the phenomenological values of α and η obtained by fitting the complete model to the experimental flux data are given below.) Comparison with (10) now yields an explicit expression of the mass fraction

$$\begin{equation} \varphi^{-1} = 1 + (1 - v^{\mathrm{sal}}_{\mathrm{c}}/v^{\mathrm{sal}})\eta/\alpha \end{equation} \tag{ 17 }$$

in terms of v^sal alone. Thereby, the problem of specifying the stationary mass balance between saltating and reptating grains has been completely reduced to the task of finding the stationary transport velocity v^sal of the saltating particles as a function of the wind strength.

3.2. Wind velocity field: two-species turbulent closure

Before we can work out the transport velocities of the two species of grains, we need to know the height-dependent wind speed u(z) and how it is affected by the presence of the airborne grains. Unfortunately, this leads us back to the question of the grain densities that we just delegated to the calculation of the grain velocities. In the past, an elaborate self-consistent calculation has been avoided by anticipating the resulting form of the wind profile on the basis of empirical observations. In the Sauermann model, for instance, an exponential vertical decay of the grain-borne shear stress is imposed,

$$\begin{equation} \tau_{\mathrm{g}}(z) = \tau_{\mathrm{g}}(0) \exp(-z/z_{\mathrm{m}}), \end{equation} \tag{ 18 }$$

in good agreement with grain-scale simulations [6, 37, 38, 46]. Since the airborne stress τ_a(z) follows from this via (4), it reduces the task to the problem of finding the mean saltation height z_m. In appendix B, we present a refined version of this approach, adapted to the two-species approach. It accounts for the contributions of the different species to the grain-borne shear stress, separately, and also for their considerably different trajectories. However, as we detail in the appendix, the whole exercise yields essentially identical results as (18), which can in fact be interpreted as the pre-averaged two-species expression. This lends further support to the general form of (18), recommending it as a suitable basis for our further discussion of the two-species model.

From (18), we obtain, via (4), the wind speed profile u(z) by integrating Prandtl's turbulence closure

$$\begin{equation} \tau_{\mathrm{a}}(z) = \varrho_{\mathrm{air}} \kappa^2 z^2 [\partial_z u(z)]^2 . \end{equation} \tag{ 19 }$$

To proceed analytically, we use a modified secant approximation similar to what was proposed by Sørensen [11]. As shown in appendix A, the closure equation then becomes

$$\begin{equation} \partial_z u \approx \frac{u_{\ast}}{\kappa z} \left[ 1 - ( 1 - u_{\ast \mathrm{t}}/u_{\ast}) {\rm e}^{-z/z_{\mathrm{m}}} \right]. \end{equation} \tag{ 20 }$$

Upon integration from the roughness length z₀ ≪ z_m, defined by u(z₀) = 0, this yields the explicit result [11, 12]

$$\begin{equation} u(z) = \frac{u_{\ast}}{\kappa}\ln(z/z_0) - \frac{u_{\ast}-u_{\ast \mathrm{t}}}{\kappa}[\mathrm{E}_1 (z_0/z_{\mathrm{m}}) - \mathrm{E}_1(z/z_{\mathrm{m}})], \end{equation} \tag{ 21 }$$

with the integral exponential $\mathrm{E}_1(z)\equiv \int _z^\infty \!\,{\mathrm{ d}} x\,{\rm e}^{-x}/x$. For small wind speeds u_* → u_*t (u_* ⩾ u_*t), the usual logarithmic velocity profile is recovered. Inside the saltation layer, a universal asymptotic wind velocity field

$$\begin{equation} u(z) \sim \frac{u_{\ast \mathrm{t}}}{\kappa} \ln(z/z_0) \quad (z_0 \lesssim z \ll z_{\mathrm{m}}) \end{equation} \tag{ 22 }$$

emerges, which is independent of the wind strength outside the saltation layer. This rationalizes the convergence of the wind profiles for different wind strengths, found in wind tunnel measurements, without the ad hoc assumption [12] of a debatable focal point [51], which resulted in the prediction of an unphysical negative flux q(u_*) at high wind speeds [12]. A direct comparison of our above prediction for u(z) with wind tunnel data from [41, 42] can be found in figure 2. From these data (available for two grain sizes d) we determine the mean saltation height z_m. Past attempts to relate z_m to d invoked a new undetermined length scale depending solely on atmospheric properties [52]. However, the results of our fit rather support the simpler (and more natural) linear relation

$$\begin{equation} z_{\mathrm{m}} \approx 25 d \end{equation} \tag{ 23 }$$

(right panel of figure 2). This is in line with recent analytical and numerical work [33], albeit possibly with a somewhat different numerical factor dependent on specific model assumptions.

Zoom In Zoom Out Reset image size

3.3. Transport kinetics: two-species transport velocities

With the wind speed profile u(z) at hand, we can now determine the characteristic stationary transport velocity v^sal of the saltating sand fraction from a standard force-balance argument. We evaluate the wind speed at a characteristic height z^sal, at which the fluid drag on the grains is then balanced with the effective bed friction, in the usual way (for a detailed derivation and discussion, see [8, 12]). Namely, the drag force on a volume element of the saltation cloud (the force per mass acting on a single saltating grain times the density ρ^sal of the saltation cloud) is balanced with τ^sal_g(0). This yields the transport velocity of the saltating sand fraction

$$\begin{equation} v^{\mathrm{sal}} = u (z^{\mathrm{sal}}) - v^{\mathrm{sal}}_\infty , \end{equation} \tag{ 24 }$$

with the estimate [12, 53]

$$\begin{equation} v^{\mathrm{sal}}_\infty \approx \frac{ \sqrt{\sigma\,{d} g/\alpha} }{ 1.3 + 41 \, \nu_{\mathrm{air}} / \sqrt{\sigma d^3\,g / \alpha}} \end{equation} \tag{ 25 }$$

for the terminal steady-state relative velocity u − v^sal. Formally, v^sal_∞ is equal to the settling velocity asymptotically reached by a grain freely falling in air (of kinematic viscosity ν_air and with grain–air density ratio σ), with a rescaled gravitational acceleration g↦g/(2α).

We fix the u_*-independent parameter z^sal by the plausible assumption that the splash at the (impact) threshold wind speed u_*t dies out, corresponding to v^sal(u_* = u_*t) = v^sal_c. The underlying physical picture is that it is the splash that keeps the saltation process going, even below the aerodynamic entrainment threshold, because the reptating particles compensate for rebound failures [5]. Accordingly, using (24) with the logarithmic wind field, equation (22), we find that

$$\begin{equation} z^{\mathrm{sal}}=z_0 \, \exp \! \left(\kappa \frac{v^{\mathrm{sal}}_{\mathrm{c}} + v^{\mathrm{sal}}_\infty}{u_{\ast \mathrm{t}}} \right) . \end{equation} \tag{ 26 }$$

In appendix C, we use a similar argument to derive the mean reptation velocity, which turns out to be almost independent of the grain size d. It can thus be approximated by the constant

$$\begin{equation} v^{\mathrm{rep}} \approx 0.7 (\sigma \nu_{\mathrm{air}} g)^{1/3}, \end{equation} \tag{ 27 }$$

for most practical purposes. This relation corresponds to a u_*-independent reptation length 2v^repv^rep_z0/g = 0.14[σ⁵ν²_aird³/g]^1/6 in the centimetre range, while the saltation length 2v^salv^sal_z0/g∝u²_*/g is quadratic in the shear velocity u_* (see figure 3) and of the order of decimetres [1, 54, 55].

Zoom In Zoom Out Reset image size

In figure 3, the transport velocities v^sal and v^rep for the two species and the mean velocity v = φv^sal + (1 − φ)v^rep are plotted over the shear velocity u_*. An important result is that the species-averaged transport velocity v is nearly constant over a broad range of wind strengths. This is consistent with the fundamental assumption on which the model is based, namely that the number of mobilized grains is sensitive to the wind strength, but their overall transport kinetics is not.

3.4. Stationary sand flux: the two-species transport law

The stationary flux law is now obtained by collecting the above results for the saturated density ρ from (11), the saltation fraction φ from (10), the transport velocities (24) with the wind velocity (21), and inserting them into (3). We introduce reduced variables, measuring the shear velocity u_* in terms of the impact threshold u_*t and the saturated flux q in units of ϱ_airu³_*/g,

$$\begin{equation} U\equiv u_{\ast}/u_{\ast \mathrm{t}}, \quad Q\equiv q g/(\varrho_{\mathrm{air}} u_{\ast}^3) . \end{equation} \tag{ 28 }$$

The result for our two-species stationary flux relation then takes the form

$$\begin{equation} Q = (1-U^{-2}) \frac{\alpha_0U + \beta_0 - \gamma_0 U^{-1}}{a U - b}, \end{equation} \tag{ 29 }$$

with the coefficients α₀, β₀, γ₀, a and b depending on the two free parameters α and η, as summarized in appendix D. It goes without saying that this result pertains to U > 1 and that Q(U < 1) ≡ 0. Except for b, which becomes negative for small η, and in particular for η = 0, all coefficients are positive. But since γ₀ < α₀ + β₀ and b < a, the flux always remains positive.

The positivity of the flux in (29) guarantees physically reasonable results, even under transient wind conditions, as occurring, for example, in sand dune simulations. The extension of the present discussion to non-stationary conditions is a major task for future investigations, in particular with regard to the saturation length, i.e. the characteristic length scale over which the system relaxes towards the steady state [8]. Its derivation within our two-species approach is of conceptual interest, since the wind strength dependence of the saturation length has recently been the subject of debates (see, e.g., [20]).

Before we determine the free parameters α and η by comparing (29) with experimental data, we first discuss the effect of our two-species parametrization on the two related transport modes—saltation and reptation—and present two different single trajectory reductions of the two-species model. Subsequently, we compare these reduced schemes to single-species transport laws that have previously been proposed in the literature.

4.1. The two-species flux balance

Our derivation of the overall sand flux Q, given in (29), also yields the partial fluxes Q^sal ≡ q^salg/(ϱ_airu³_*) and Q^rep ≡ q^repg/(ϱ_airu³_*) of the saltating and reptating species, respectively. Their dependence on the wind velocity is illustrated in figure 4. It is frequently argued in the literature that the saltating grains dominate the sand transport, for they perform large jumps, while the reptating sand fraction might be negligible in terms of mass transport. Our model basically confirms this view for moderate winds, not too far above the impact threshold. However, it predicts that both species contribute almost equally to the total flux, under strong winds. The shift in the mass balance φ towards reptation (figure 4, right plot) reflects the dependence of the number of ejected grains on the wind strength, which follows from the splash-impact relation (15) if one inserts the transport velocities of the individual species, as given in figure 3. The variable contribution of the two species to the overall sand flux is a result that cannot be obtained from the conventional single-species models, but should be of interest for some of the more advanced applications mentioned in the introduction.

Zoom In Zoom Out Reset image size

4.2. One-species limits

It is interesting to see how the conventional one-species descriptions emerge from the two-species model upon contraction to a single effective grain species. This contraction is clearly not unique but may be performed in different ways. On the basis of the two-species flux balance, discussed in the preceding section, two natural contractions suggest themselves: one for moderate wind speeds where the flux is dominated by saltation, and one for strong winds, where the ratio of the contributions from reptating and saltating grains was found to be roughly constant.

The first single trajectory reduction, where the reptating species is dismissed, is obtained by setting φ = 1 or η = 0, which corresponds to the flux relation

$$\begin{equation} Q = (1-U^{-2})( a - b \, U^{-1} ), \end{equation} \tag{ 30 }$$

with a negative coefficient b < 0 (appendix D). Since this equation accounts only for high-energy saltating grains, the resulting absolute flux is too large compared with the original two-species description. This can be corrected by reducing the effective trajectory height or, in our formalism, by rescaling the effective height z^sal at which the fluid drag is balanced with the effective bed friction. The lower dotted line in figure 5 was obtained in this way.

Zoom In Zoom Out Reset image size

The second single trajectory reduction is motivated by the weak dependence of the mass fraction φ on u_* for strong winds u_* ≫ u_*t, which suggests replacing φ by its u_*-independent limit 1/(1 + η/α). This yields

$$\begin{equation} Q = (1-U^{-2}) ( \tilde \alpha_0 + \skew4\tilde \beta_0 \, U^{-1}), \end{equation} \noindent \tag{ 31 }$$

with

$$\begin{equation} \tilde \alpha_0 = \frac{\alpha_0}{a} \quad \mathrm{and} \quad \skew4\tilde \beta_0 = \frac{\beta_0}{2a} + \frac{\tilde \eta}{1+\tilde \eta} \frac{ v^{\mathrm{rep}}_{z0}}{u_{\ast \mathrm{t}}}, \end{equation} \noindent \tag{ 32 }$$

and the abbreviation $\tilde \eta \equiv \eta v^{\mathrm{rep}}/v^{\mathrm{rep}}_{z0}$. If one interprets the (now constant) species mass balance as a free phenomenological parameter, it may be adjusted by fine-tuning η, so that the absolute value of the sand flux predicted by this formula agrees better with experimental observations. This is how we obtained the upper dotted line in figure 5.

Both (30) and (31) have the functional form of the transport law proposed by Durán and Herrmann [12] based on the focal point assumption. Note that the effective height z^sal is small compared with the saltation height, where the wind speed obtained from the focal point assumption is very weak. As a consequence, the mean drag force is reduced, which results in a negative parameter a < 0 in (30) and a negative sand transport rate Q < 0 for large shear velocities U > |a/b|. Table 1 and figure 5 provide a summary of various stationary aeolian sand transport laws that have been proposed in the literature, in comparison with the prediction of our two-species model and its above single-trajectory contractions.

4.3. Comparison with experiments

We now fit our two-species flux relation Q(U) from (29) to the empirical data from various wind tunnel measurements [40, 43], using η and α as free fit parameters. As demonstrated in figure 6, we obtain excellent agreement with the experimental data for a grain-size-independent splash efficiency η and an effective coefficient of restitution α that increases with increasing grain size d. The transport rates for four different grain sizes collected in [40] were obtained using sand traps; the data in [43] were obtained by particle tracking techniques for a single grain size only. The different methods yield different absolute values for the sand flux, which might be due to systematic differences in the detection efficiency. In the model, these differences are reflected in the (apparently) different efficiencies with which saltating grains generate reptating grains, i.e. different values of the parameter η. We find that η = 3.8 for the data obtained by particle tracking and η = 9 for the sand trap data, independently of the grain size. But both data sets were consistent with the same formula for the bed restitution coefficient

$$\begin{equation} \alpha = 1 - d_0/d, \end{equation} \noindent \tag{ 33 }$$

with d > d₀ = 88 μm (left panel of figure 7). The rebound becomes increasingly elastic for larger grains, while α vanishes for smaller grains at a critical grain diameter d = d₀, which is in accord with the observation that saltation only occurs for sand grains larger than about 70 μm [56]. The dependence of the restitution coefficient α on the grain size is phyically plausible, for the collision of smaller grains is increasingly influenced by hydrodynamic interactions (and also by cohesive and electrostatic forces). In other words, d₀ is a characteristic grain size marking the transition from a phenomenology typical of sand to one typical of dust.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Moreover, while we also fine-tune the threshold u_*t for each grain size by hand when fitting the experimental data, the values used turn out to be in very good agreement with the expectation $u_{\ast \mathrm{t}} \approx 0.1\sqrt {\sigma g d}$, as demonstrated by the inset of figure 6. Altogether, the stationary sand transport rate observed in experiments is thus convincingly reproduced by the two-species result, equation (29), over a wide parameter range, with very plausible and physically meaningful values for the model parameters.

The observed value of d₀ can be rationalized by a hydrodynamic order of magnitude estimate. Commonly, one relates the crossover from sand-like to dust-like behaviour to the difference in transport modes of these two classes of grains. While sand is transported by grain hopping dust remains suspended in the air for a while. Balancing the settling velocity of the grains by the typical (upward) eddy currents, which are of the order of u_*t, one obtains a minimal sand grain size of about 1.5(ν²/σg)^1/3 ≈ 30 μm [1]. To match d₀ = 88 μm, the eddy velocity would have to be 4.5u_*t. Two further estimates of d₀ may be obtained as follows. The first assumes that the crossover marked by d₀ is concomitant with the crossover from a more elastic to a more viscous collision between individual grains. In this case, the relevant quantity is the Stokes number, which characterizes the particle inertia relative to the viscous forces. Using the observations of the collision experiments by Gondret et al [57] for the critical Stokes number St, below which the impacting particles do not rebound from a rigid wall, one gets a critical grain size of about 20 St^2/3ν^2/3/σ/g^1/3 ≈ 20 μm. To match d₀ = 88 μm, the critical Stokes number would have to be of the order of 190, i.e. ten times larger than the value obtained in [57], for which experimental differences might possibly be blamed. Note, however, that the most pertinent argument should be the one that gives the largest value of d₀, as it provides a lower bound for the grain size contributing to saltation. A larger estimate of d₀ is indeed obtained by observing that sand grains collide with the bed, and they are lifted from the bed by turbulent fluctuations in the wind velocity. In contrast, suspended and resting dust particles are protected from bed collisions and turbulent lift, respectively, by the laminar boundary layer that coats any solid (no slip) boundary. The characteristic length scale of this laminar coating is related to the so-called Kolmogorov dissipation scale (ν³_air/ε)^1/4 [58], where ε is the dissipation per unit mass of the driving medium (in our case air under normal conditions). For the logarithmic wind profile near the ground, equation (22), the scale-dependent dissipation takes the form ε(z) = u³_*t/κz [59]. A self-consistency argument requiring that the Kolmogorov dissipation scale at height z above the ground equals z (if no externally imposed roughness scale is available) yields z ≃ [ν³_air/ε(d₀)]^1/4. Extrapolating the logarithmic profile with $u_{\ast \mathrm{t}} = 0.1 \sqrt {\sigma g d}$ down to the scale z, we obtain z = 4.6(ν²_air/σg)^1/3 ≈ 10² μm. We find this estimate to be in excellent agreement with the data in the left panel of figure 7, which suggests

$$\begin{equation} d_0 \approx 4 (\nu_{\mathrm{air}}^2/\sigma g)^{1/3}. \end{equation} \tag{ 34 }$$

We have developed a 'second generation' continuum description of aeolian sand transport, relaxing the strong mean-field approximation inherent in the classical single-trajectory models such as [8] and introducing a two-species framework, as advocated in [5]. We started from the physics on the grain scale and corroborated our explicit analytical expressions by a comprehensive comparison with empirical data for the splash, the wind velocity field and the sand flux, covering a broad range of ambient conditions and grain sizes. While the usefulness of our general approach of contracting the complicated splash and saltation process into a drastically reduced two-species picture was strongly suggested by empirical and theoretical work, it necessarily involved some free phenomenological coefficients. For our two-species model, these are, apart from a few minor numerical coefficients listed in table D.1, the effective bed restitution coefficient α(d) for the saltating grains of diameter d and the parameter η that fixes the number of reptating grains ejected by the average impacting saltating grain. Both parameters were found to take consistent and plausible values if used as free parameters when fitting the empirical flux data obtained in wind tunnel experiments. In particular, the size dependence of the restitution coefficient fits perfectly to the phenomenological criteria commonly used to distinguish dust from sand. The predicted two-species stationary flux law (29) was found to be in excellent agreement with comprehensive data from different sources. It should provide an excellent analytical starting point for a variety of advanced applications calling for a more faithful description of the saltation process so far available—from wind-driven structure formation in the desert to saltation-driven dust production and emission.

We thank Thomas Pähtz and Jasper Kok for inspiring discussions and for making the manuscript of [33] available to us prior to its publication.

With (4) and the assumption of an exponential decay of the grain-borne shear stress with height, equation (18), the modified Prandtl turbulence closure (19) reads as

$$\begin{equation} \partial_z u = \frac{u_{\ast}}{\kappa z} \sqrt{1-\frac{\tau_{\mathrm{g}}(0)}{\tau} {\rm e}^{-z/z_{\mathrm{m}}}}. \end{equation} \noindent \tag{ A.1 }$$

To make analytical progress, we follow Sørensen [11] in approximating the square root by a secant. The right-hand side of (A.1) has the functional form $\sqrt {1 - \epsilon x}$, with $x \equiv \exp (-z/z_{\mathrm{m}})$ and  ≡ τ_g(0)/τ, which we approximate by a secant of the form ax + b(1 − x). Matching the points {0,1} and $\{1,\sqrt {1-\epsilon }\}$ yields b = 1 and $a=\sqrt {1-\epsilon }$, hence

$$\begin{equation} \sqrt{1 - \epsilon x} \approx 1- (1- \sqrt{1-\epsilon})x. \end{equation} \noindent \tag{ A.2 }$$

Under stationary conditions, using (7), we have $\sqrt {1-\epsilon }=\sqrt {(\tau -\tau _{\mathrm{g}}(0))/\tau }=u_{\ast \mathrm{t}}/u_{\ast }$, which leads to (20). This avoids an artefact of Sørensen's original approximation [11], $\sqrt {1 - \epsilon x} \approx 1 - \epsilon x$, which implies $\sqrt {\tau _{\mathrm{a}}(0)/\tau } = 1 - \tau _{\mathrm{g}}(0)/\tau$ for the shear stress at the ground, as already criticized by Durán and Herrmann [12].

In (5) of section 2, the grain-borne shear stress was split into the contributions from saltating and reptating particles, respectively. The ratio of the grain-borne shear stresses at the ground,

$$\begin{equation} \frac{\tau_{\mathrm{g}}^{\mathrm{rep}}(0)}{\tau_{\mathrm{g}}^{\mathrm{sal}}(0)} = \alpha \frac{ v^{\mathrm{rep}}}{v^{\mathrm{rep}}_{z0}} \frac{1-\varphi}{\varphi} \equiv \tilde \alpha \, \frac{1-\varphi}{\varphi}, \end{equation} \noindent \tag{ B.1 }$$

immediately follows from (8) and (12), and from our simplifying assumption that the reptating grains are ejected vertically. Assuming the exponential decay of the grain-borne shear stress τ_g with height to hold for both components, the Prandtl turbulent closure reads as

$$\begin{equation} \rho_{\mathrm{air}} \kappa^2 z^2 (\partial_z u)^2 = \tau - \tau_{\mathrm{g}}^{\mathrm{rep}}(0) {\rm e}^{-z/ z^{\mathrm{rep}}_{\mathrm{m}}} - \tau_{\mathrm{g}}^{\mathrm{sal}}(0) {\rm e}^{-z/ z^{\mathrm{sal}}_{\mathrm{m}}}. \end{equation} \noindent \tag{ B.2 }$$

We exploit the strong scale separation z^sal_m/z^rep_m ≃ 10² between the characteristic jump heights of saltating and reptating grains [1, 60], on which the two-species model is based. (The precise value turns out to be irrelevant to our discussion.) It allows the closure to be solved for two separate height ranges: (i) z < z^rep_m ≪ z^sal_m associated with reptation, where we may set $\exp (-z/ z^{\mathrm{sal}}_{\mathrm{m}}) \to 1$, and (ii) z ≫ z^rep_m, associated with saltation, where we may set $\exp (-z/ z^{\mathrm{rep}}_{\mathrm{m}}) \to 0$. Applying the secant approximation for the square root as described in appendix A, we can perform the integrations within both ranges and match the asymptotic solutions at z = z^rep_m. Using (7) and (B.1) to eliminate τ^sal_g(0) and τ^rep_g(0), this yields

$$\begin{equation} \fl u(z) = \left\{\begin{array}{@{}ll@{}} \displaystyle \frac{S u_{\ast}}{\kappa} \ln\frac{z}{z_0} - \frac{S u_{\ast} - u_{\ast \mathrm{t}}}{\kappa}\left[ \mathrm{E}_1\!\left(z_0/z^{\mathrm{rep}}_{\mathrm{m}}\right) - \mathrm{E}_1\!\left(z/z^{\mathrm{rep}}_{\mathrm{m}}\right) \right] , & z \leqslant z^{\mathrm{rep}}_{\mathrm{m}}, \\ \displaystyle u(z^{\mathrm{rep}}_{\mathrm{m}}) + \frac{u_{\ast}}{\kappa} \ln\frac{z}{z^{\mathrm{rep}}_{\mathrm{m}}}- \frac{u_{\ast} (1 - S )}{\kappa} \left[ \mathrm{E}_1\!\left(z^{\mathrm{rep}}_{\mathrm{m}}/z^{\mathrm{sal}}_{\mathrm{m}}\right) - \mathrm{E}_1\!\left(z/z^{\mathrm{sal}}_{\mathrm{m}}\right) \right] , & z > z^{\mathrm{rep}}_{\mathrm{m}}, \end{array}\right. \end{equation} \tag{ B.3 }$$

where we abbreviated

$$\begin{equation} S = \sqrt{1 - \frac{(1-\tau_{\mathrm{t}}/\tau)\varphi}{\varphi + (1-\varphi)\tilde \alpha}} . \end{equation} \tag{ B.4 }$$

Using this in the force balance, (24), we gain the implicit equation

$$\begin{equation} v^{\mathrm{sal}} = u(z^{\mathrm{sal}},v^{\mathrm{sal}}) - v^{\mathrm{sal}}_\infty , \end{equation} \tag{ B.5 }$$

where the fraction φ of saltating grains itself depends on the velocity v^sal of the saltating grains via (10). To solve it, we note that the fraction φ of saltating particles is a decreasing function of the shear stress τ, bracketed by 1 and 0 (figure 4, right panel). Hence, the product (1 − τ_t/τ)φ is a small number for all τ ⩾ τ_t. Setting it to zero in (B.3), we obtain

$$\begin{equation} \fl u(z) \sim \frac{u_{\ast}}{\kappa} {\rm ln} \left(\frac{z}{z_0}\right) - \frac{u_{\ast}-u_{\ast \mathrm{t}}}{\kappa} \times \left\{ \begin{array}{@{}ll@{}}\lbrack {\rm E}_{1}(z_{0}/z^{\rm rep}_{\rm m})-{\rm E}_{1}(z/z^{\rm rep}_{\rm m}) \rbrack , & z \leqslant z^{\rm rep}_{\rm m},\\ \lbrack {\rm E}_{1}(z_{0}/z^{\rm rep}_{\rm m}) -{\rm E}_{1}(1) \rbrack , & z > z^{\rm rep}_{\rm m}. \end{array}\right. \end{equation} \tag{ B.6 }$$

Since the exponential integral E₁(z/z_m) vanishes rapidly for increasing z > z_m, we recover the wind speed profile (21) successfully employed in the one-species models, with the characteristic decay height given by z_m ≡ z^rep_m. Inserting (B.6) into (B.5) yields an affine increase of v^sal with u_*, in good accord with the exact numerical solution of (B.5) (figure B.1, right panel).

Zoom In Zoom Out Reset image size

For the numerical solution, one has to deal with the two free parameters z^rep_m and α. While the former is directly related to the wind profile, the latter is a fit parameter of the model, determined from a comparison of the predicted sand transport rate with empirical data (figure 6). To make progress, we vary z^rep_m and take the value of α from section 3.4, where the sand flux is estimated by means of the pre-averaged approach for the wind profile. For a grain diameter d = 242 μm, we obtained α ≈ 0.63 (i.e. $\tilde \alpha=1.4$), which is consistent with collision experiments, as argued in [12]. From the numerical solution of (B.5), we find good agreement between the self-consistently gained u(z) and wind tunnel measurements [42] for z^rep_m = 25d, which is exactly the same value as obtained in section 3.2 within the pre-averaged approach. This supports our observation that (B.3) can be approximated by taking the limit φ → 0 in (B.6), which yields the pre-averaged wind profile for z_m ≡ z^rep_m. This result is almost independent of the ratio z^sal_m/z^rep_m, for (B.6) is independent of z^sal_m. Thereby, we formally confirm the intuitive expectation that the feedback of the grains on the wind profile is predominantly due to the many reptating particles and hardly affected by the saltating particles.

As to the saltating grain fraction in section 3.3, we estimate the transport velocity of the reptating grains from the grain scale physics, i.e. from the hop-averaged horizontal velocity of an individual grain. (Note that we do not introduce a new variable to distinguish the grain-scale velocity from the mean transport velocity, because the context prevents confusion.) The time-dependent velocity of a reptating grain obeys the drag relation

$$\begin{equation} \partial_t v^{\mathrm{rep}} = \frac{g}{(v^{\mathrm{rep}}_\infty)^2}\, \left| u-v^{\mathrm{rep}} \right| \, \left(u-v^{\mathrm{rep}}\right), \end{equation} \tag{ C.1 }$$

similar to that for saltating grains [8, 12]. For saltating grains, an additional friction force (besides the drag force) would appear on the right-hand side of the equation of motion, representing the mean loss of momentum upon rebound. But, since the reptating grains perform only a single hop, such a friction term does not enter (C.1). We assume that the ejection is essentially vertical with the initial velocity v^rep_z0 of the order of u_*t. A more accurate discussion would not substantially change our findings, as confirmed by the numerical solution of (C.1). Note that the nonlinearity of the drag law entails a time-dependent 'terminal settling velocity' v^rep_∞, dependent on the actual relative grain velocity u − v^rep (e.g. [53]). However, for our purpose, and in view of the low reptation trajectories, we can safely approximate the reptation velocity from (C.1) by inserting the wind speed u(z^rep) at a given reptation height and the steady-state terminal velocity

$$\begin{equation} v^{\mathrm{rep}}_\infty = \sqrt{\sigma g d} [0.95 + 20 \, \nu_{\mathrm{air}} /\sqrt{\sigma g d^3}]^{-1} \end{equation} \tag{ C.2 }$$

derived from the effective drag law proposed in [53], similar to (25) (see also [8, 12]). Neglecting moreover vertical drag forces, the maximum height of the reptation trajectory is (v^rep_z0)²/(2g) ≈ 10d < z_m. Consequently, we may insert the ground-level wind field (22) and obtain the mean reptation velocity

$$\begin{equation} v^{\mathrm{rep}} \approx \frac{u_{\ast \mathrm{t}} \ln(z^{\mathrm{rep}}/z_0)}{\kappa \left[ 1 + \kappa \; {(v^{\mathrm{rep}}_\infty)}^2/\left(2v^{\mathrm{rep}}_{z0} u_{\ast \mathrm{t}} \ln(z^{\mathrm{rep}}/z_0)\right) \right]}, \end{equation} \tag{ C.3 }$$

where we approximated the time of flight by that for a parabolic trajectory, 2v^rep_z0/g, as usual. Inserting the empirical observation v^rep_z0 ≈ u_*t as well as u²_*t = 0.01σgd and z₀ = d/20, see table D.1, results in an estimate for the reptation velocity v^rep as a function of the grain diameter d, illustrated by a solid line in figure C.1, which we may identify with the (mean field) transport velocity of the reptating grain fraction.

Zoom In Zoom Out Reset image size

The yet undetermined value of the effective reptation height z^rep is expected to be of the order of the maximum height (v^rep_z0)²/(2g) of the trajectory. Indeed, if we compare the approximate result given by (C.3) with the numerical solution v^rep(t) of (C.1) averaged over the whole trajectory,

$$\begin{equation} \overline{v^{\mathrm{rep}}} = \frac{g}{2 v^{\mathrm{rep}}_{z0}} \int_0^{2v^{\mathrm{rep}}_{z0}/g} {\rm d} t\; v^{\mathrm{rep}}(t) , \end{equation} \tag{ C.4 }$$

we obtain a good match for large grain diameters d > 800 μm for z^rep = 0.14u²_*t/(2g). Note, however, that the numerical solution yields an almost d-independent reptation velocity $\overline {v^{\mathrm{rep}}} \approx 0.5\,\mathrm{m}\,{\rm s}^{-1}$ for the relevant grain sizes $d \approx 100\,\mu \mathrm{m}\dots 1\,\mathrm{mm}$, as illustrated in figure C.1. To find an estimate for $\overline {v^{\mathrm{rep}}}$, which still captures the dependence on other parameters appearing in (C.3), we evaluate this equation for a representative grain diameter within the relevant range. A closer look at (C.3) reveals that the reptation velocity at the inflection point d ≈ 6.7ν^2/3_air(σg)^−1/3 ≈ 150 μm can be approximated by v^rep ≈ (ν_airσg)^1/3 = 0.68 m s⁻¹, which provides the wanted parameter dependence. To better match the absolute values found from the numerical solution, we insert a factor of 0.7 by hand, thus arriving at the analytical estimate given in (27).

Here we give explicit expressions for the coefficients occurring in the saturated sand flux (29):

$$\begin{equation} a = 2\alpha (1 + \tilde \eta) \frac{v^{\mathrm{sal}}_{\mathrm{c}} + v^{\mathrm{sal}}_\infty}{u_{\ast \mathrm{t}}} - 2\alpha \frac{1+\tilde \eta}{\kappa}\mathcal{E} \! \left[ z_0\exp \! \left( \kappa \frac{v^{\mathrm{sal}}_{\mathrm{c}} + v^{\mathrm{sal}}_\infty}{u_{\ast \mathrm{t}}}\right) \right], \end{equation} \tag{ D.1 }$$

$$\begin{equation} b = a - 2 \alpha \frac{v^{\mathrm{sal}}_{\mathrm{c}}}{u_{\ast \mathrm{t}}} \end{equation} \tag{ D.2 }$$

$$\begin{equation} \alpha_0 = \frac{a^2}{(1 + \tilde \eta)^2}, \end{equation} \tag{ D.3 }$$

$$\begin{equation} \beta_0 = \frac{2 \tilde \eta a}{1 + \tilde \eta} \left( \frac{v^{\mathrm{rep}}_{z0}}{u_{\ast \mathrm{t}}} + \frac{a}{1 + \tilde \eta} - \frac{b}{\tilde \eta} \right), \end{equation} \tag{ D.4 }$$

$$\begin{equation} \gamma_0 = \frac{\tilde \eta v^{\mathrm{rep}}_{z0}}{u_{\ast \mathrm{t}}} \left[\frac{\tilde \eta v^{\mathrm{rep}}_{z0}}{u_{\ast \mathrm{t}}} + 2(a-b) \right] - \left( \frac{1 + \tilde \eta}{2a} \beta_0 \right)^2 \end{equation} \tag{ D.5 }$$

with the abbreviations

$$\begin{equation} \mathcal{E}(z) = {\rm E}_1(z_0/z_{\mathrm{m}}) - {\rm E}_1(z/z_{\mathrm{m}}) \end{equation} \tag{ D.6 }$$

and

$$\begin{equation} \tilde \eta \equiv \eta \, v^{\mathrm{rep}}/v^{\mathrm{rep}}_{z0}. \end{equation} \tag{ D.7 }$$

All parameters occurring in these equations are listed in table D.1. As explained in the main text, these quantities are determined as follows. From collision experiments, we obtain the numerical values of the critical impact velocity v^sal_c and the vertical ejection speed v^rep_z0. The roughness length z₀ and the height z_m (the characteristic decay height of the air-borne sand density) are estimated by fitting experimentally observed wind velocity profiles above the saltation layer. Finally, the parameters η, α and the threshold shear velocity u_*t are determined by fitting the sand transport rate to wind tunnel measurements.

Table D.1. The parameters occurring in equations (D.1)–(D.5), which are either numerical constants or dependent on the grain diameter d, the gravitational acceleration g = 9.81 m s⁻¹, the sand–air density ratio σ = 2163 or the kinematic viscosity of air, ν_air = 1.5 × 10⁻⁵ m² s⁻¹.

Parameter	Value or formula
u*t	$0.1 \sqrt {\sigma g d}$
κ	0.4
η	9 [40], 3.8 [43]
vsalc	10u*t
vrepz0	u*t
z0	d/20
zm	25 d
α	1−d0/d,  d0=4(ν2air/σg)1/3
vsal∞	$\sqrt {\sigma g d/\alpha } [1.3 + 41 \sqrt {\alpha } \, \nu _{\mathrm{air}} /\sqrt {\sigma g d^3}]^{-1}$
vrep	0.7(σνairg)1/3