A two-species continuum model for aeolian sand transport

Starting from the physics on the grain scale, we develop a simple continuum description of aeolian sand transport. Beyond popular mean-field models, but without sacrificing their computational efficiency, it accounts for both dominant grain populations, hopping (or"saltating") and creeping (or"reptating") grains. The predicted stationary sand transport rate is in excellent agreement with wind tunnel experiments simulating wind conditions ranging from the onset of saltation to storms. Our closed set of equations thus provides an analytically tractable, numerically precise, and computationally efficient starting point for applications addressing a wealth of phenomena from dune formation to dust emission.


Introduction
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. But 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. Onespecies models have become a very successful means of gaining analytical insight into [9,10,11,12,13,14] and to conduct efficient large-scale numerical simulations of [15,16,17,18,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 both concerning the way how 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 noticeable a systematic overestimation of the stationary flux at high wind speed (see figure 5, below). Therefore, the one-species models have been criticized for their lack of numerical accuracy and internal consistency [5,20]. There are 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,4,2,23,24,25,26,27,28,29]. But the amount 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,31,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 different hop lengths of the grains of different sizes.
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 a more detailed experimental and theoretical knowledge of the grain-scale physics [5,34,35,36,37,38,39] and rely on more comprehensive empirical information about the wind shear stress and sand transport to test our predictions [40,41,42,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 Figure 1. The two-species picture of aeolian sand transport maps the ensemble of mobile grains onto two representative species. High-energy ("saltating") grains accelerated by the wind (arrows in the wind-speed profile indicating the wind velocity) rebound upon impact and eject some low-energy ("reptating") grains in a splash. Sand transport is quantified in terms of the transport velocities, v sal and v rep , and mean densities, ρ sal and ρ rep , of the two species. The momentum balance of the splash process (inset, see section 3.1) is effectively encoded in the impact velocity v sal , the horizontal and vertical rebound velocities of the saltating grain, v sal 0 and v sal z0 , respectively, and the number N and velocity v rep z0 = v rep 0 of the ejected reptating grains (neglecting a small horizontal velocity component, for simplicity).
in section 5.

Basic notions of aeolian sand transport and the two-species parametrization
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, viz. the local mass density ( x, t) and sand transport velocity v( 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 from integrating the volume density over z. (The vertical grain distribution is a crucial element of the modified turbulent closure, though, which accounts for the feedback of the sand onto the wind.) In this contribution, 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 Introducing the dimensionless mass fraction ϕ ∈ [0, 1] of saltating grains relative to the total amount of mobile grains, the (integrated) densities of the mobile grains can be written in the form ρ sal = ϕρ and ρ rep = (1 − ϕ)ρ .
Accordingly, the flux balance and the relative contribution of reptating and saltating grains to the flux read as 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 up into two components carried by the airborne grains ("g") and the air ("a") itself, respectively, 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,38,39,46]. In the twospecies 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 up the grain-borne shear stress 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], 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 Both stress contributions are defined as the product of the vertical sand flux φ i and the horizontal velocity 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 i 0 and v i z0 , respectively, and obtain Since the average ejection angle of the reptating particles is known to be about 80 • [34], v rep 0 ≈ 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.
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 with the average number N of reptating grains ejected by an impacting saltating grain.
Together with the first relation in (8), we then have 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 As often found in the literature, we have employed the notion of the shear velocity or friction velocity u * , defined by τ ≡ air u 2 * , 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 of the saltating grains [8]. Because of the complexity of the splash process, which makes first-principle quantitative estimates of α forbiddingly complex, we suggest to treat 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(σν air g) 1/3 /u * t , we can rewrite the density of the transported sand in the more explicit form 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 D1. Our result (see table 1) accounts for essential elements of the aeolian sand transport process beyond the singletrajectory approximation, and turns out to be in remarkable agreement with wind tunnel measurements.

Implementing the two-species framework
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 * ).

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 grainborne 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 windstrength-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 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 of 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 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 sal where we used Bagnold's estimate u 2 * t ≈ 0.01σgd [1]. A concise summary of the empirical input entering our model is provided in table D1.
Summarizing the above discussion, we can rewrite the momentum balance (14) as 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 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 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.

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, 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 get, via (4), the wind speed profile u(z) by integrating Prandtl's turbulence closure 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 Upon integration from the roughness length z 0 z m , defined by u(z 0 ) = 0, this yields the explicit result [11,12] with the integral exponential E 1 (z) ≡ ∞ z dx 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 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 (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.

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 forcebalance 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 with the estimate [12,53] v sal ∞ ≈ σdg/α 1.3 + 41 ν air / σd 3 g/α 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 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 v rep ≈ 0.7(σν air g) 1/3 , for most practical purposes. This relation corresponds to a u * -independent reptation length 2v rep v rep z0 /g = 0.14[σ 5 ν 2 air d 3 /g] 1/6 in the centimetre range, while the saltation Figure 3. The transport velocities of the saltating (dashed) and reptating (dotted) sand fraction versus the rescaled shear velocity U ≡ u * /u * t, as predicted by (24) and (27), respectively. The solid line represents the average transport velocity in a reduced effective single-trajectory description.
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.

Stationary sand flux: 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 air u 3 * /g, U ≡ u * /u * t , Q ≡ qg/( air u 3 * ) .
The result for our two-species stationary flux relation then takes the form with the coefficients α 0 , β 0 , γ 0 , 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 γ 0 < α 0 + β 0 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]).

Discussion
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.

The two-species flux balance
Our derivation of the overall sand flux Q, given in (29), also yields the partial fluxes Q sal ≡ q sal g/( air u 3 * ) and Q rep ≡ q rep g/( air u 3 * ) 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 splashimpact 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.

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 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. 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 to replace ϕ by its u * -independent limit 1/(1 + η/α). This yields and the abbreviationη ≡ ηv rep /v 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. Table 1. Scaling functions f (U ) for stationary aeolian sand transport relations of the form Q(U ) = 1 − U −2 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 α 0 , β 0 , γ 0 , a and b of the two-species model can be found in Appendix D.

This work
This work, one-species limit  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.

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 Figure 6. The rescaled saturated flux Q = qg/( air u 3 * ) for various grain sizes (d = 125 µm, 170 µm, 242 µm, 320 µm, 242 µm from bottom to top): data from wind tunnel measurements using particle tracking methods [43] (stars) and sand traps [40] (other symbols) are compared with the prediction by the twospecies model, equation (29) (solid curves). The restitution coefficient α served as a global but grain-size-dependent fit parameter (see figure 7), while η = 9 and η = 3.8 was fixed for the sand trap and the particle tracking data, respectively. (The open symbols represent several rescaled data sets, obtained for a variety of bed slopes.) with d > d 0 = 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 0 , 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 0 is a characteristic grain size marking the transition from a phenomenology typical of sand to one typical of dust.
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 * t ≈ 0.1 √ σgd, 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 0 can be rationalized by a hydrodynamic order of magnitude estimate. Commonly, one relates the crossover from sand-like to dustlike behaviour to the difference in transport modes of these two classes of granular. 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(ν 2 /σg) 1/3 ≈ 30 µm [1]. To match d 0 = 88 µm, the eddy velocity would have to be 4.5u * t . Two further estimates of d 0 may be obtained as follows. The first assumes that the crossover marked by d 0 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 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 0 , as it provides a lower bound for the grain size contributing to saltation. A larger estimate of d 0 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 (ν 3 air /ε) 1/4 [58], where ε is the dissipation per unit mass of the driving medium (in our case air at normal conditions). For the logarithmic wind profile near the ground, equation (22), the scale-dependent dissipation takes the form ε(z) = u 3 * 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 [ν 3 air /ε(d 0 )] 1/4 . Extrapolating the logarithmic profile with u * t = 0.1 √ σgd down to the scale z, we obtain z = 4.6(ν 2 air /σg) 1/3 ≈ 10 2 µm. We find this estimate to be in excellent agreement with the data in the left panel of figure 7, which suggests

Conclusions
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 D1, 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. needed to eject any grains, over the shear velocity: self-consistent numerical solution of (B.5) (symbols) and its analytical prediction in the limit (1−τt/τ )ϕ → 0, for which u(z) is given by (B.6), (line).
borne shear stress τ g with height to hold for both components, the Prandtl turbulent closure reads We exploit the strong scale separation z sal m /z rep m 10 2 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 to solve the closure for two separate height ranges: (i) z < z rep m z sal m associated with reptation, where we may set exp(−z/z sal m ) → 1, and (ii) z z rep m , associated with saltation, where we may set exp(−z/z rep m ) → 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 where we abbreviated Using this in the force balance, (24), we gain the implicit equation 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 Since the exponential integral E 1 (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 B1, right panel).
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. α = 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.

Appendix C. Reptation velocity
As to the saltating grain fraction in 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 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 site 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  (27), respectively. The latter can be understood as an average over the relevant grain sizes (white background).
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 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 ) 2 /(2g) ≈ 10d < z m . Consequently, we may insert the ground-level wind field (22) and obtain the mean reptation velocity v rep ≈ u * t ln(z rep /z 0 ) 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 2 * t = 0.01σgd and z 0 = d/20, see table D1, 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 C1, which we may identify with the (mean field) transport velocity of the reptating grain fraction.
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 ) 2 /(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, we obtain a good match for large grain diameters d > 800 µm for z rep = 0.14u 2 * t /(2g). Note, however, that the numerical solution yields an almost d-independent reptation velocity v rep ≈ 0.5 m s −1 for the relevant grain sizes d ≈ 100 µm . . . 1 mm, as illustrated in figure C1. To find an estimate for v 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 −1 , which provides the wanted parameter dependence. To better match the absolute values found from the numerical solution, we insert a factor 0.7 by hand, thus arriving at the analytical estimate given in (27).
All parameters occurring in these equations are listed in table D1. 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 0 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.