The apparent surface roughness of moving sand transported by wind

We present a comprehensive analytical model of aeolian sand transport in saltation. It quantifies the momentum transfer from the wind to the transported sand by providing expressions for the thickness of the saltation layer and the apparent surface roughness. These expressions are for the first time entirely derived from basic physical principles. The model further predicts the sand transport rate (mass flux) and the impact threshold shear velocity. We show that the model predictions are in very good agreement with experiments and numerical state of the art simulations of aeolian saltation. PACS numbers: 45.70.Mg, 47.27.-i, 47.55.Kf, 92.10.Wa, 92.40.Gc, 92.60.Gn Submitted to: New J. Phys.


Introduction
Saltation is the dominant mechanism of aeolian sand transport on Earth's deserts under turbulent wind flow. Unidirectional wind accelerates sand grains, which perform hops of typical shapes. During their hops the wind continuously transfers momentum to the grains. Therefore the wind momentum decreases, resulting in reduced wind velocities not only within, but also above the saltation layer. Above the saltation layer the average horizontal wind velocity profile u(z) follows the well known law [1], where u * is the wind shear velocity, κ = 0.4 is the von Kármán constant, and z * o is the apparent roughness of the moving saltation layer. In the absence of sand transport z * o becomes z o , which is the surface roughness of a quiescent sand bed. In the presence of sand transport the magnitude of z * o depends on how much momentum is absorbed by the saltation layer. It is crucial for the development of sand transport models, but also for landscape modelers and coastal managers to know z * o as function of u * and z o . For instance in aeolian dune models z * o is a key quantity in the computation of the wind field over a non-flat topography, in which the shear velocity u * varies with the spatial position [2,3,4,5]. The main purpose of this paper is to derive a novel prediction of z * o , which is entirely based on physical principles. Deriving a scaling law for z * o was also approached by previous studies. First Owen [6] suggested where g is the gravity constant. (2) is also known as the Charnock relation, since Charnock [7] derived (2) for the roughness of a wind-blown water surface. Owen [6] based his formula on the assumptions that the average lift-off velocity v l , with which a grain leaves the bed, is proportional to u * , that opposing drag forces can be neglected, and that z * o is proportional to the average hop height of grains h, also called saltation height. However the author's assumptions fail to agree with measurements. Experimental studies [8,9,10,11] found v l and consequently h to be almost independent of u * within the measured range. Furthermore z * 0 cannot be proportional to h, since z * o was found to strongly vary with u * in experiments [12,13,14] in contrast to h. In addition Sherman [15] found that (2) leads to strong discrepancies with experiments close to the impact threshold u t , which is the threshold shear velocity at which saltation can be sustained through the splash process. On Earth u t is below the fluid entrainment shear velocity, needed to entrain sand grains from the soil by fluid lift [1,16]. Sherman [15] therefore extended (2) to the so-called modified Charnock relation, which ensures that at the threshold u * = u t , where no particles are moving, the roughness is unchanged, z * o = z o . Although (3) was successfully validated with the data set of Sherman and Farrell [14] for z * o , it shares the same lack of physical justification as (2).
A much more physical approach was presented by Raupach [17]. From the mixing length approximation [18], the author derived where z s is the decay height of the grain shear stress profile τ g (z), which the author assumed to be exponentially decreasing, τ g (z) = τ go e −z/zs , where τ go = τ g (0). τ g (z) describes how much momentum is transferred, at each height z, from the fluid to the grains per unit soil area and time. The difficulty in the usage of (4) is the undetermined quantity z s . Raupach [17] therefore assumed z s to be proportional to the saltation height h, which he in turn assumed to be proportional to u 2 * /g like before Owen [6]. However, as we discussed before, Owen's assumption is in disagreement with experiments. Relations very similar to (4) were also obtained in two further studies [19,20]. These studies achieved agreement with the experimental data of Rasmussen et al. [12], by introducing an ad-hoc fit relation for z s . Andreotti [19] found that the data set can be well fitted if z s scales with √ d, where d is the mean particle diameter, and he therefore suggested where s = ρ s /ρ w is the ratio of sand density ρ s and fluid density ρ w , and µ is the kinematic viscosity. A scaling law very similar to (6) was also used by Duran and Herrmann [20]. However (6) is very weakly founded on physics. Its only justification is the resulting agreement of (4) with the data set of Rasmussen et al. [12]. Therefore there is a great necessity to either validate (6) or to derive a new expression for z s from physical principles. Within this study we do the latter. We present a comprehensive analytical model of aeolian saltation, which aims to significantly improve previous analytical models [1,6,21,22,23]. It for the first time provides expressions for z s and z * o entirely derived from physical principles. Our analysis will reveal that z s is a measure for the thickness of the saltation layer and not proportional to the average hop height h. The model furthermore incorporates expressions for other important sand transport quantities, such as the mass flux Q and the impact threshold u t . The model is based on the concept of mean motion, meaning that average quantities are used for its description. In our model we separately consider the horizontal and vertical transport of grains. For each we analytically balance the average force and work rate per unit soil area applied by the wind on a grain during a trajectory versus the respective amounts applied by the soil on the grains during an impact. This results in a model parameter α, describing the ratio between the average vertical and horizontal force per unit soil area, and another model parameter β, describing the ratio between the average work rate per unit soil area in the vertical motion and the horizontal motion. With theoretical arguments it is shown that α and β are nearly independent of u * and the atmospheric conditions and only slightly varying with the buoyancy-reduced gravityg and the particle diameter d. The model further contains a third parameter γ, defined as the ratio between z s and the effective height of the mean motion z m . The final relations for z s , z * o , and Q are functions of α, β, γ, and u t . We afterwards extend our model in such a way that also u t can be computed. Thereby a fourth parameter η comes into play, describing the ratio between the average particle velocity, reduced by the particle slip velocity, and the average wind velocity. In our model we use the assumption that the grain shear stress profile is exponentially decaying, equivalent to what was assumed in previous studies [17,20] (see (5)). We validate our model with the apparent roughness data of Rasmussen et al. [12] for five different grain sizes, with combined impact threshold data of several studies [24,25,26], and with mass flux data of Creyssels et al. [10]. Furthermore, through simulations with the numerical state of the art model of Kok and Renno [27] we support our derived expressions and our statement that the model parameters are nearly independent of u * , the atmospheric conditions, as well asg and d.
The manuscript is structured in the following way. It starts with a comprehensive model description in Section 2, which is followed by the model validation in Section 3 and a discussion of the results in Section 4. The appendices incorporate long calculations and side information. There is also a glossary at the end of the manuscript, which helps to keep track of the mathematical symbols.

Model description
It is the main purpose of our paper to derive a novel expression for the apparent roughness z * o during aeolian sand transport in steady state. The main focus of our model lies therefore in the analytical description of the momentum and energy transfer from the wind to the grains. Momentum and energy transfer are the main causes for the increase of the surface roughness z o of a quiescent sand bed to the apparent roughness z * o of a moving saltation layer. In detail we use Newton's law to obtain equations, which balance the average force and work rate per unit soil area applied during a grains trajectory with the respective amounts applied during an impact. The ratio of the average force (work rate) per unit soil area for the vertical motion and force (work rate) per unit soil area for the horizontal motion is the definition of our model parameter α (β). After applying the balance laws we show that the decay height z s of the grain shear stress profile τ g (z) and subsequently z * o as well as the mass flux Q can be calculated from the impact threshold u t and our model parameters. Afterwards the model is extended in such a way that also u t can be calculated as function of the model parameters. As previous studies [17,20] we assume an exponentially decreasing τ g (z) (see (5)) and extensively use it in our calculations. This assumption is therefore discussed in a separate paragraph.
For the balance laws, we only consider wind drag and gravity as driving forces, but neglect turbulent lift forces, the Magnus force, electrostatic forces, and momentum as well as energy changes by mid-air collisions between grains for simplicity reasons and because gravity and drag dominate the sand transport [27]. Furthermore we simplify the description, by only considering average quantities, which implies that we neglect turbulent fluctuations of the wind velocities. Further simplifications are the use of monodisperse, spherical sand grains, being transported above a horizontal sand bed. The probably most crucial of all these simplifications is the negligence of mid-air collisions between saltating grains. The effect of such collisions has only rarely been subject of scientific studies [28,29,30,31], because for a comparison between saltation with and without mid-air collisions, one has to turn off mid-air collisions, what is possible (and common) in numerical simulations, but impossible in experiments. According to the most recent numerical study ( Figure 5 in Ren and Huang [31]), the change of the mass flux due to mid-air collisions is less than 10% for u * ≈ 3.5u t . This is below the typical measurement error of mass flux measurements (> 10%). According to this study the neglegance of mid-air collisions and therefore our model simplifications are acceptable up to at least u * ≈ 3.5u t .
This section is separated in several subsections. It starts with the presentation of notations and definitions, which are used for the description of our model, in Section 2.1. In Section 2.2 follows a short discussion of our main model assumption, the exponentially decreasing grain shear stress profile. After that the balance laws are applied, first for the force per unit soil area in Section 2.3 and then for the work rate per unit soil area in Section 2.4. Subsequently we discuss the invariance of the model parameters α and β in Section 2.5. Afterwards we obtain a novel relation for z s in Section 2.6, which is further discussed in Section 2.7. Then in Section 2.8 relations for z * o and Q as a function of u t and the model parameters are obtained. In Section 2.9 the model is extended, in order to allow for the computation of u t as well.

Notations and definitions
For the coming analytical calculations, we henceforth use the following notations: An index x refers to the horizontal component of a given quantity, which coincides with the direction of the wind, an index z to the vertical direction, whereby z is also the height above the sand bed. Furthermore we differentiate between the upward and downward part of a grain's trajectory by indices ↑ and ↓, respectively. Quantities evaluated at the sand bed z = 0 incorporate an additional index o. In particular quantities, which refer to a grain's impact, consist of the indices o and ↓, if the quantity is evaluated before the impact, and the indices o and ↑, if the quantity is evaluated after the impact.
In order to keep the manuscript simple, it is advantageous to predefine quantities, which are used in the following calculations. One quantity is the average particle mass per unit volume ρ(z), transported at height z. ρ(z) integrated over the whole saltation layer describes the mass M of transported sand per unit soil area.
Since we differentiate between upward and downward movement ρ(z) can be divided in the mass of upward and downward moving particles per unit volume Other important quantities are the average vectorial wind velocity profile, u(z), whose z-component is zero u(z) := u x (z) = |u(z)|, and the average vectorial particle velocity profile for the upward (downward) part of the trajectory v ↑(↓) (z). The difference between both velocities is denoted as Based on these definitions, we further define the following velocity differences by and ∆v 2 as well as the local vertical mass flux φ(z) by where we used that the vertical upward flux must exactly compensate the downward flux in steady state. Note that v z↓ (z) and thus ∆v z (z) are negative. Using (8), (14) can be rewritten as With these definitions the average gain of horizontal and vertical momentum of a transported grain per unit soil area and time between the two times it crosses height z can be written as for the horizontal and for the vertical momentum gain per unit soil area and time. τ g (z) is also known as the grain shear stress profile [17,20,23] and p g (z) can be seen as a grain normal stress (grain pressure) profile. Note that, by inserting (15) in (17), one obtains which is identical to the definition of the granular pressure in previous studies [10,32], if vertical drag is neglected.

Grain shear stress profile
In this section we discuss and motivate our main model assumption of an exponentially decreasing grain shear stress profile (see also (5)) (19) is justified in the following manner. First, an approximately exponentially decreasing mass density profile ρ(z) was measured in wind tunnels [10,11]. Although not necessarily identical, the profiles τ g (z) and ρ(z) should at least behave in a similar manner. Therefore it is very reasonable that also τ g (z) decreases approximately exponentially. Second, τ g (z) has been obtained from numerical simulations [27,33], which indeed showed an approximately exponential decrease. This is shown in Figure  1 for simulation results with the numerical model of Kok and Renno [27]. It should be noted that the mass density profile ρ(z) strongly deviates from the exponential shape at very small heights in the simulations. Such a deviation is also present for the grain shear stress profile τ g (z) (see Figure 1 at heights very close to zero), however to a much lesser extent.

Force balance
As already pointed out, the description of the momentum transfer from the wind to the grains is a key ingredient towards a description of the feedback effect of sand transport on the wind profile and thus a first step towards a prediction of the apparent roughness z * o . Newton's second law for grains moving in a particular trajectory, indicated by a lower index '1', can be written as (dz = v 1z↑(↓) dt) for the upward (downward) part of this trajectory, where f 1x↑(↓) and f 1z↑(↓) are the horizontal and vertical components of the total average force f 1↑(↓) per unit volume acting on the grain in the upward (downward) part of this trajectory. For the single- [34]. Summing the upward and downward part of (20) and (21), respectively, followed by averaging over all trajectories and integration over height therefore yields where we approximated the trajectory average of products as the product of trajectory averages in φ∆v x and φ∆v z . Here ∆v xo = ∆v x (0) and ∆v zo = ∆v z (0), and f x (z) and f z (z) are the averages of f 1x = f 1x↑ + f 1x↓ and f 1z = f 1z↑ + f 1z↓ over all trajectories, respectively. The terms on the left hand side of (22) and (23) describe the average horizontal and vertical force per unit soil area applied by the soil on the grains during an impact and the right hand side the average horizontal and vertical force per unit soil area applied by the wind on the grains during a trajectory, respectively. In the next steps we evaluate the integrals in (22) and (23). Therefore we first need an expression for the total trajectory-averaged force f. Since we neglect the Magnus force, turbulent lift forces, and momentum transfer through collisions, f is only composed of the drag and the gravity force. Further approximating the trajectory average of products as the product of trajectory averages, f can be written as where e z is the unit vector in z-direction,g = s−1 s g with s = ρ s /ρ w is the buoyancyreduced gravity (for most atmospheresg ≅ g), and v r↑(↓) = |v r↑(↓) |. C d is the drag coefficient, which is a function of the particle Reynolds number and therefore of v r↑(↓) . For many drag laws in the literature e.g. [35,36] the dependency of C d on a velocity difference V can be described by a law of the type where C o and C ∞ = C d (∞) are dimensionless parameters. Such a drag law strongly simplifies f z , which becomes where (14), (16), and v rz↑(↓) = −v z↑(↓) were used. If a drag law of another type than (25) was used e.g. [37], (26) would still be valid in very good approximation. On the other hand we rewrite f x as where denotes a weighted average of a quantity f between the upward and downward movement, f = (ρ ↑ f ↑ + ρ ↓ f ↓ )/ρ. Now we can evaluate the integral in (22) using the expression we derived for f x . We obtain where the overbar and the capital letters denote the average of a quantity f over height, We further used (7), the approximation V r ≅ V rx , which is reasonable since in aeolian saltation the horizontal motion dominates the vertical one, and we approximated the height average of the products by the products of the height On the other hand p go can now be calculated from our expression for f z . It becomes where we used our assumption (19). With the evaluation of the integrals, we can now define α ′ and our first model parameter α as where we used (16) and (17) as well as the notations ∆v xo = ∆v x (0) and ∆v zo = ∆v z (0). The advantage of this definition of α lies in the fact that α ′ is almost independent of u * , atmospheric conditions,g, and d as we show later in a separate chapter. For many conditions, we can approximate since α ′ is typically much larger than C ∞ z s /(sd), as will be verified later. It mainly means that the gravity force is large in comparison to the vertical drag force, f z ≈ −gρ.
We can thus formulate relevant sand transport quantities as a function of a constant α.
For instance from (31) we obtain a direct relation between the average velocity difference V r and α, writing and further, using (29), (30), and (32), a direct relation between the grain shear stress τ go at the bed and the mass of transported grains per unit soil area M, writing

Work rate balance
The second important ingredient towards a description of the feedback effect of the grain motion on the wind profile and towards a prediction of z * o is the description of the energy transfer from the fluid to the grains. Since we discuss a purely Newtonian problem, we separate the horizontal and vertical motion. The work rate balance with respect to the horizontal (vertical) motion can be obtained by multiplying (20) ( (21)) with v 1x↑(↓) (v 1z↑(↓) ), summing the upward and downward part, integrating over height, and averaging over all trajectories. It yields where ∆v 2 xo = ∆v 2 x (0), ∆v 2 zo = ∆v 2 z (0), and we approximated the trajectory average of products as the product of trajectory averages in φ(0)∆v 2 xo and φ(0)∆v 2 zo . The terms on the left hand side of (35) and (36) describe the average work rate during an impact and the right hand side the average work rate during a trajectory for the horizontal and vertical motion, respectively. Analogous to (24) and (27) we can now write where we used (14) in (38). Analogous to (28), integration now approximately yields where V = V x describes the average particle velocity. Note that (40) . Evaluating the integrals, we now define β ′ as (41) means that the average granular temperature V 2 z is proportional to V r V . This is different from Creyssels et al. [10] who found V 2 z to be approximately equal to v 2 z (0). The main reason for this difference is that the authors neglected vertical drag, whereas our description considers it (see (40)). As before for α ′ , the advantage of (41) lies in the fact that β ′ is almost independent of u * , atmospheric conditions,g, and d as we will show in the following.

Invariance of α ′ and β ′
Since our model relations, including the final relation for z * o , will be expressed as functions of α ′ and β ′ , it is important to discuss, how these parameters change with varying conditions. As can be seen from (30) and (41), both parameters are ratios of certain velocity differences evaluated at the soil z = 0 and therefore related to the splash-entrainment process. The splash-entrainment process dominates the entrainment of bed grains in aeolian steady state saltation, because the entrainment by wind is small due to a strong reduction of the wind velocity close to the sand bed, which even leads to decreasing near-surface velocities with increasing u * [20]. Since fluid-entrainment is not relevant, each impacting grain must exactly lead to one grain leaving the surface (rebound or ejection of new grains) on average for steady state sand transport. The average number grains leaving the surface per impacting grain can however only depend on the average impact velocity v i , angle θ i , and the relevant bed propertiesg and d. Furthermore, for given values ofg and d, the average velocity v l and angle θ l of a grain leaving the surface after an impact, called lift-off velocity and angle, can only depend on v i and θ i . This means, for given values ofg, d, and θ i , there are unique values v i , v l , and θ l , which fulfill that one impacting grain makes one grain leave the surface on average. Another necessary condition, which must hold, is that v zo↓ = v i sin θ i must be smaller than v zo↑ = v l sin θ l due to friction with the air, it would be equal in the absence of drag. The validity of this condition was observed in saltation experiments [38]. Further, from collision experiments Oger et al. [39] found that this condition is only fulfilled for small impact angles θ i 15 • . The authors also found that the number of ejected particles significantly decreases with decreasing θ i , meaning that the range of θ i -values, in which saltation can be sustained, should be rather narrow. In fact, it has been measured in experiments [8] that the average impact angle is approximately constant, θ i ≈ 11 • , with increasing u * , and thus the other splash quantities as well stay approximately constant with varying u * . Here and henceforth we refer to situations above the impact threshold u * ≥ u t and within our model limits (not too large u * ), when mentioning dependencies on u * . This means that in particular the parameters α ′ and β ′ are approximately independent of u * and atmospheric properties, and thus only functions ofg and d.
In order to get an idea of how α ′ and β ′ behave as functions ofg and d, we use the model of Kok [40], with which v zo↓ and v zo↑ and thus α ′ and β ′ can be computed. The model is briefly described in Appendix A. The model results in functions α ′ (gd) and β ′ (gd) plotted in Figure 2. For the plots the gravities of Earth,g = 9.81m/s 2 , and Mars,g = 3.71m/s 2 , were fixed and d varied. As can be seen, even as functions ofg and d the variance of α ′ and β ′ is small according to this model. Furthermore α ′ is of the order of unity, what justifies the approximation (32) for many conditions. For instance for Earth conditions the neglected term can be estimated as being of the order of 1/25, since z s is the same order of magnitude as the saltation height h, which has been measured as being equal to about 40d [11], C ∞ ≈ 1, and s ≈ 1000 on Earth.

A novel relation for z s
The definitions of the parameters α and β ′ obtained from the momentum and energy balances can now be used to express the decay height z s of the grain shear stress profile τ g (z) as a function of α and β ′ . As already pointed out z s is the key quantity towards a prediction of the apparent roughness z * o . For the calculation of z s we use where we inserted (14). We further approximate the arithmetic average of |v z↑ | and |v z↓ | by their geometric average what is reasonable since |v z↑ |/|v z↓ | < |v zo↑ |/|v zo↓ | and |v zo↑ |/|v zo↓ | is about 1.5 as measurements indicate [39], which means that the error of this approximation is less than 3%. Using (14), (16), (42) and (43), we express τ g as This allows us to rewrite (22) as where we used (19). Since (45) is valid for all heights z, it must be particularly valid for the average over height. We therefore approximately calculate z s as where β = β ′ ∆V x /(2V ) and we used (31) and (41). β is our second model parameter and like β ′ approximately constant, since we expect that ∆v x is in leading order proportional to V . (46) is the main contribution of our paper, since it is, to our knowledge, the first physically based prediction of z s and therefore the most important part towards a novel physically based prediction of z * o . If the values of the model parameters α and β are known, V remains the only undetermined quantity in (46), since V r can be calculated by (33) as a function of α. There is evidence that z s does not only describe the decay of τ g (z), but also the decay of ρ(z) for large z. This can be seen from (45), which says that τ g (z) decays in the same way as ρ C d (v r )v r v rx . Since C d (v r )v r v rx is only slightly decaying with height for large z, the decaying behaviors of ρ(z) and τ g (z) are very similar to each other. This is shown in Figure 3 for simulations with the numerical model of Kok and Renno [27]. Since the profile ρ(z)/M describes the hop height distribution of saltating grains, z s is related to the saltation height h, which is the average hop height of the grains. This is discussed in detail in the following section.

Physical meaning of z s
It was assumed in previous studies that z s and the saltation height h are proportional to each other [17,19,20]. This very natural assumption is however not valid for the saltation simulated with the numerical model of Kok and Renno [27], since the normalized profiles ρ(z)/ρ o for u * = 0.3m/s and u * = 0.8m/s in Figure 3 almost coincide with each other at small heights. This means, although z s is larger for u * = 0.8m/s, the hop height of grains transported close to the surface is almost the same and hence h increases weaker with u * than z s . It is therefore reasonable to interpret z s as the height of high-energy saltons and the height z r up to which the profiles ρ(z)/ρ o coincide as the height of low-energy saltons [19], which remains unchanged with u * . In the following we explain the reason for this behavior of ρ(z).
In our model ρ(z) decays approximately exponentially, if and only if v 2 z does not vary much with height z. This can be seen from (23), (26), and (32), which allow us to write the differential equation, using p g = −ρ v 2 z (analogous to (44)), whose solution is an exponential decrease, if and only if v 2 z is constant with z. However, there can be a huge difference between the value v 2 z (0) at the soil, which is fixed by the splash-entrainment process, and the value of v 2 z at larger heights, which is proportional to z sg . Figure 3 shows indeed a strong deviation of ρ(z) from the exponential shape at very small heights within the low energy layer z ≪ z r . In the light of our analysis, this deviation corresponds to a strong increase of v 2 z (z) from v 2 z (0) towards about z sg at larger heights. Note that wind tunnel studies [10,11] did not notice such a deviation from the exponential shape in their measurements of ρ(z). A possible cause is that their lowest measurement points, z ≈ 20d and z ≈ 40d, respectively, were already too high, and they measured only in a region, where v 2 z (z) was not varying much anymore. Note further that a very similar reason, namely that − v x v z (0) is fixed by the splash-entrainment process, is the probable cause for the very slight deviation of τ g (z) from the exponential shape at very small heights (see Figure 1). Finally note that mass flux profiles ρ(z) v x (z), which are often measured in experiments [8,9,41,42,43], decay slightly weaker than ρ(z), since the particle velocity v x (z) increases weakly with height.
After explaining ρ(z), we now calculate h and z r . First, using (41) and (47) as well as partial integration, h can be calculated as Since z r does not depend on u * and since it must be of the same structure as z s and h, namely proportional to v 2 z (z)/g at a typical height z, and since z = 0 is the only height where v 2 z (z)/g does not depend on u * due to the splash-entrainment process, z r must write where we used (14). Our analysis confirms the picture of Andreotti [19], who hypothesized that one can essentially distinguish two species in aeolian saltation, low-energy saltons (reptons) slowly moving in small hops, and high-energy saltons moving fast in huge hops. The author hypothesized that z r is a measure for the height of the focal-region, a region in which steady state wind profiles for different shear velocities u * intersect.
In this section we showed that the saltation layer can be characterized by three heights. The height of low-energy saltatons z r , the saltation height h, and the height of high-energy saltons z s , which is also a measure for the thickness of the saltation layer. In contrast to the first height, the latter two change with V and thus with u * . The prediction of V is therefore subject of the following section.

Calculation of z * o
The last step towards the calculation of z * o is to derive an expression for V , the last undetermined quantity in (46), our relation for z s . z s is the main parameter in our final relation for z * o , which will be of a similar structure as (4), the relation of Raupach [17]. For deriving V we use the following strategy. We first approximately calculate U as the wind velocity at an height z m , which denotes the height of the mean motion, U = u(z m ), with the mixing length approximation [18]. Then we compute V by where we use that V r does not change with u * (see (33)).
Following the outlined strategy, we calculate U from the mixing length approximation [17,18,20,23,33] as with In the absence of sand transport, τ g (z) = 0, the mixing length approximation yields the undisturbed logarithmic velocity profile ((1) with z * o = z o ). In the presence of sand transport, τ g (z) = 0, the velocity profile deviates from the logarithmic shape. In Appendix B u(z) as well as z * o are calculated based on our model assumption of an exponentially decreasing grain shear stress profile, (19), and the calculation approximately yields and for z > 0.1z s where u b is the reduced wind shear velocity at the bed, defined by and the exponential integral E 1 (x) as well as G(x) are defined by Note that (54) is only an approximation of u(z) for z > 0.1z s . In Appendix B one can also find an approximation for z < z s , which however has a much more complicated structure. For the coming calculations, we are only interested in the value U = u(z m ), for whose calculation both approximations perform similarly well, since z m is between 0.1z s and 0.4z s as we show later. In (53) and (54) the shear velocity at the bed u b remains undetermined. It was shown by Duran and Herrmann [20] that u b must decrease with u * starting from u b (u t ) = u t , however with a small slope. The reason is that one observes a focal region in aeolian steady state saltation, called Bagnold-focus, below which the wind velocities decrease with increasing u * [1,44]. Andreotti [19] argued that this strong decrease of wind velocities, the presence of the Bagnold-focus, and consequently the value of u b are mainly caused by grains transported in the low-energy layer. As we explained in Section 2.7, at very small heights within the low-energy layer z ≪ z r the mass density profile ρ(z) and also the grain shear stress profile τ g (z) deviate from the exponential decrease. We however used this exponential decrease of τ g (z) (see (19)), our main model assumption, for the computation of the apparent roughness and the wind profile in (53) and (54). This does not mean that (53) and (54) are wrong, because the deviation of τ g (z) from the exponential shape is very small (almost invisible in Figure 1). But it means that we cannot use the real value of u b , which is influenced by the low-energy layer. Instead we must use a value of u b , which corresponds to the extrapolation of the exponential shape of τ g (z) above z = z r to the height z = 0. Such a value of u b would be larger than the real value, because ρ(z) and thus τ g (z) deviate from the exponential shape towards higher values within the low-energy layer. We therefore propose that this value is close to u t , This is supported by simulations using the numerical model of Kok and Renno [27]. Figure 4 shows that wind profiles calculated by (54) with the hypothesis (59) are much closer to the simulated profiles than those in which the simulated values of u b were used, and this although this hypothesis eliminates the Bagnold-focus. Especially in the region which we are interested in, z > 0.1z s , (53) and (54) provide an excellent approximation of the simulated wind profile, if using (59). Note that (59) is also known as Owen's second hypothesis, and it has been used in many previous models e.g. [6,17,21,22,23]. From (53), (54), and (59) we now obtain and where γ is the third model parameter and defined by The first terms on the right hand side of (53) and (60) are identical to the relations of previous studies [17,20] (see also (4)). The second term appears, because we did not approximate the right hand side of (51) before the integration as done in these studies. Our solution is therefore more precise. Note that the height of the mean motion z m should be in leading order proportional to the decay height z s of the grain shear stress profile, because the mean motion is dominated by the motion of high-energy saltons [19], and z s is a measure for the height of high-energy saltons (see Section 2.7). Consequently γ is in leading order constant. Having obtained a relation for U, we can calculate V with (50), insert V in (46) to obtain z s , and insert z s in (60) to obtain z * o . The calculation of z * o can therefore be summarized as where we used (60) and (62) for (63a), (46), (50), and (62) for (63b), and (63d) and (63c) are the same as (33) and (61), respectively. (63b) and (63c) can be solved iteratively for U and z m , and (63d) can be solved iteratively given a certain drag law C d (V ). This is our novel prediction of z * o in the most general version. If the impact threshold u t is known, z * o can be calculated using (63a-63d) as function of the model parameters α, β, and γ. It is furthermore possible to compute the mass flux Q as function of the same parameters. For this purpose we first compute the mass of transported sand per unit soil area M from (34) and (59) as Then the mass flux Q becomes where we inserted (64) for M and (50) for V . U and V r can be computed by (63b-63d).
At the moment all model equations are functions of the model parameters and u t . In order to close the model, u t must be calculated as function of the model parameters as well. This is done in the following with the help of a closing assumption.

2.9.
Closing the model -a relation for u t (63a-63d) are already a novel expression for z * o . Like previous expressions in the literature (see (2)-(4) and (6)) it needs the impact threshold u t as an input parameter. We therefore close the model by deriving a relation for u t in this section. The strategy is as follows. We first motivate a simple expression for the particle velocity at the threshold, V t = V (u t ). We then compute V r = U t − V t , where U t = U (u t ), and combine the result with our previous expression for V r , (63d). The resulting equation can be rearranged to compute u t .
As outlined before, we motivate the following closing expression, where z mt = z m (u t ), V o = (ρ o↑ v xo↑ + ρ o↓ v xo↓ )/ρ o is the average particle slip velocity (i.e., the particle speed at the surface), where v xo↑(↓) = v x↑(↓) (0) and ρ o↑(↓) = ρ ↑(↓) (0), and η is the fourth model parameter. (67) means that the difference between average particle and slip velocity under threshold conditions is proportional to the average wind velocity. This is justified in the following manner. The particle slip velocity V o is a quantity that like the model parameters α and β only depends on the impactentrainment process. In particular V o is independent of the average wind velocity U t . From theoretical and experimental studies it is known that the profile of the average horizontal particle velocity v x (z) starts with V o at z = 0 and increases with height z [10,11,27]. The average increase with z mainly depends on the average wind velocity U t . Consequently, the average particle velocity V t is a function of the average wind velocity U t plus an offset V o . The simplest possible relation with such a behavior is given by (67). Thereby η describes how efficiently the wind accelerates transported grains under threshold conditions. Rearranging (67), η can be written as We calculate the particle slip velocity V o with the model of Kok [40], explained in Appendix A. The result, V o as function ofgd, is plotted in Figure 5. Now we can use (67) to calculate V r . According to (33), V r does not depend on u * . We can therefore write Rearranging and using (46), (62), and (67) finally yields an expression for u t , which can be summarized as

Model Validation
The model is validated by simulation results with the numerical state of the art model of Kok

Independence of the model constants of atmosphere and grain properties
For the verification of the independence of the model parameters we use simulation results of the numerical model of Kok and Renno [27] for conditions, which are summarized in Table ( [27] uses the drag law of Cheng [37], namely which we use to compute V r in (63d) and (70c). We evaluate (46) Table 1 versus the values of M , obtained directly from the simulation for Earth conditions with d = 500µm (blue) and Mars conditions with d = 250µm (red). Each circle belongs to a different shear velocity u * . The solid line indicates perfect agreement. The parameter u t , which appears in (64), was not calculated by our model, but directly obtained from the simulations (see Table 1). d = 250µm, Figure 6 shows M calculated using (64) with α as given in Table 1 versus the values of M, obtained directly from the simulations. Furthermore, for the same conditions, Figure 7 shows z s and ln(z * o /z o ) calculated using (46), (63a), and (63d) with β as given in Table 1 versus the values of z s and ln(z * o /z o ), obtained directly from the simulations. And finally, for the same conditions, Figure 8 shows V calculated using (50), (63c), and (63d) versus the values of V , obtained directly from the simulations. For all plots, each circle belongs to a different shear velocity u * and the solid line indicates perfect agreement. Conditions different from those plotted in Figs. (6)(7)(8) show in most cases the same good agreement. The values of α, β, and γ for all conditions are given in Table ( 1). It shows that α is between 0.9 and 1 for most of the tested conditions. A variance of α and thus M of about 10% is however small compared to the degree of uncertainty one usually faces in saltation mass(flux) measurements. Furthermore, β is between 0.12 and 0.135 for all tested conditions, the variance of β is therefore even less than the variance of α. Furthermore, Table (1) shows that γ is between 0.23 and 0.34 for all of the tested conditions. Therefore we can confirm that α, β, and to a lesser extend γ can indeed be used as approximately constant parameters for saltation simulated by the model of Kok and Renno [27] at least within the range of conditions displayed in Table  (1). We want to emphasize that in particular (46), which is the main contribution of our paper, well describes the behavior of the simulated decay heights z s (see Figure 7). Note that the slight disagreement of the Mars simulations from the perfect agreement in Figure 8 is probably due to turbulent fluctuations of the wind velocities, which are  (46), (63a), and (63d) with β as given in Table 1 versus the values of z s and ln(z * o /z o ), obtained directly from the simulations for Earth conditions with d = 500µm (blue) and Mars conditions with d = 250µm (red). Each circle belongs to a different shear velocity u * . The solid line indicates perfect agreement. The parameters u t and α are taken from Table (1) and the simulated values of V are used in (46), instead of calculating them with the parameter γ.  Table 1 versus the values of V , obtained directly from the simulations for Earth conditions with d = 500µm (blue) and Mars conditions with d = 250µm (red). Each circle belongs to a different shear velocity u * . The solid line indicates perfect agreement. The parameters u t and α are taken from Table (1). The parameter β is not used, instead the simulated values of z s are taken for the computation of V .  Table (1). considered by the numerical model of Kok and Renno [27], but not by our analytical model. These fluctuations are much more important for smallgd like Mars conditions with d = 250µm than for largegd like Earth conditions with d = 500µm.
Finally we check (69) by plotting V r +V o over U t = κ −1 u t ln(z mt /z o ) for all simulated conditions. Thereby V r is calculated by (63d) with the values of α in Table (1) and the values V o are also given in Table (1). Note that in the simulations V o does not change significantly with u * as we also stated before from a theoretical point of view. Figure 9 shows that the plotted circles, each of them corresponding to one of the conditions in Table (1), approximately lie on a straight line through the origin. This indicates an approximately universal behavior of η, because according to (69) V r + V o and κ −1 u t ln(z mt /z o ) are proportional to each other with a proportionality constant 1−η. The values of η are also given in Tab. (1).
In conclusion, the simulation results of the numerical state of the art model of [27] can be very well described by our analytical model indicating only slight variances of all four model parameters with changing conditions. Note that the simulated values of α, β, and V o do not follow the predictions of the model of Kok [40], which were plotted in Figs. (2) and (5).
We also want to emphasize that we entirely failed to fit the numerical data for the apparent roughness z * o , when using (6), the scaling for z s proposed by Andreotti [19]. Earth and Mars conditions could not be fitted simultaneously by (4) and (6), or alternatively (63a) and (6) with a single proportionality constant in (6). Fitting to good agreement with the numerical Earth data, led to a disagreement with the Mars data by almost two orders of magnitude, even after trying several modifications like for instance replacing √ sgd in (6) by u t . This strongly indicates that (6) does not describe the physics sufficiently well. Andreotti [19] mentioned himself that (6) is nothing but a guess: "We have not found any simple explanation of this scaling with ν = µ/ρ w . The only indication that we have identified the good parameter is the prefactor of order unity." 3.1.1. Explicit relation for mass flux Q of monodisperse sand, d = 250µm, on Earth Many mass flux relations in the literature are given for free field Earth conditions with d = 250µm e.g. [23]. In this section, we provide a further explicit prediction of Q for these conditions, based on the parameter values in the third row in Table ( 1). In contrast to the parameter values which we obtained from wind tunnel experiments discussed in the following section, the values given in Table (1) correspond to free field conditions, because the numerical model of Kok and Renno [27] was adjusted to such conditions. From α = 0.94 we obtain V r = 1.55m/s using (63d) and (71). From u t = 0.196m/s, β = 0.125, and γ = 0.33 we further obtain z mt = 0.0153m using (63b) evaluated at u * = u t . Inserting all values in (66) then yields where we used (63a-63c). This is however not an explicit relation for Q, since z m increases with u * as well. In order to obtain an explicit relation, we approximate z m ≅ z mt . This is reasonable, since the increase of Q with z m is only logarithmic. Further inserting z o = d/30 (see Table 1) and using that the first non-vanishing term 0.0125 ln 2(1−u 2 t /u 2 * ) 2 of G(u t /u * ) (see Appendix B) is very small compared to the other terms close to u t , G ≈ 0, we finally obtain

Experimental validation
For the validation with experiments we use the drag law (25) of Cheng [37], which we also used before. Furthermore, we need to account for the fact that the surface roughness z o of a quiescent sand bed is a function of the roughness Reynolds number. This is in particular important when comparing with experimental impact thresholds, because some experiments were made with very small particle diameters d in the aerodynamically smooth regime. The whole context is explained in Appendix C including equations, which are used to compute z o . In this section we validate our apparent roughness prediction, (63a-63d), with the experiments of Rasmussen et al. [12] for five different particle diameters d. The chosen data set has the advantage that the scatter in the data is small in comparison to other data sets [13,14]. Furthermore we use a combination of several data sets [24,25,26], in order to evaluate our impact threshold prediction, (70a-70c) and the data set of Creyssels et al. [10] for our mass flux prediction (66). The latter choice is motivated by the fact that Creyssels et al. [10] used particle tracking methods, which are more accurate than measurements of the mass flux Q with sand traps, which underestimate Q by up to 50% [41,45]. Furthermore the experiments were performed in the same wind tunnel as the experiments of Rasmussen et al. [12] and mass flux measurements typically vary from wind tunnel to wind tunnel. For instance the two recent measurements of the mass flux with particle tracking methods [10,11] show the same qualitative behavior, an approximate scaling of Q with u 2 * − u 2 t , but quite different magnitudes of Q, although they used the same sand in their experiments. We want to strongly emphasize that we use only one single set of parameters, α, β, γ, and η to fit all data sets at the same time and that the values of u t , obtained from our prediction (70a-70c), are used in (63a-63d) and (66). By fitting the model parameters to α = 1.02, β = 0.095, γ = 0.17, and η = 0.1 we obtain good to excellent agreement with all data sets. This is shown in Figs. (10)(11)(12), which present the comparison of (63a-63d) with the data of Rasmussen et al. [12], the comparison of (70a-70c) with the impact threshold data sets [24,25,26], and the comparison of (66) with the mass flux data of Creyssels et al. [10], respectively. It is important to note that the fitted values of the model parameters significantly differ from those in Table (1), which were adjusted to match the numerical data obtained from simulation with the model of Kok and Renno [27]. The very likely cause is that the numerical model of Kok and Renno [27] was adjusted to field data, whereas our model is adjusted to wind tunnel data. According to Sherman and Farrell [14], field and wind tunnel data of the apparent roughness z * o can differ by up to one order of magnitude. The reason for this is not fully understood. It was speculated by Raupach [17] that differences between wind tunnel and field data are attributed to not fully equilibrated saltation in wind tunnels. It is however likely that this is not the only reason, since even in large wind tunnels like the one used by Creyssels et al. [10] and Rasmussen et al. [12], for which we adjusted our model, the differences still exist. We propose that differences in the grain size distribution could be another cause. Sand in field experiments has typically a much broader grain size distribution than sand with the same mean diameter d typically used in wind tunnels. For instance the five different sands used by Rasmussen et al. [12] are the same as the ones used by Iversen and Rasmussen [46], which the latter described as "closely sized sand samples" and referred to as "uniform samples". Bagnold [21] argued that the average particle velocity V in aeolian saltation is found to be about two times higher for broadly in comparison to narrowly distributed grain sizes with the same d, because "elastic rebound of smaller grains off more massive ones is superimposed on the splash process", leading to higher saltation heights. According to (46), two times higher V in field experiments than in wind tunnel experiments would lead to about three times higher z s and therefore to a much higher apparent roughness z * o .

Discussion
We have presented a comprehensive analytical saltation model, which provides a new level of understanding of the change of the average wind profile during aeolian saltation and of aeolian saltation in general. The model incorporates (63a-63d), a novel expression for the apparent roughness z * o , and (46), a new expression for the thickness of the saltation layer z s . (46) can be seen as the main contribution of our study and it is to our knowledge the first relation for z s entirely derived from physical principles. It is an improvement of the relation (6), proposed by Andreotti [19], which had no physical foundations. Besides, our saltation model also provides (66), a novel expression for the mass flux Q, and (70a-70c), an expression for the impact threshold u t . All relations are extensively evaluated. They are in very good agreement with the recent numerical model of Kok and Renno [27] and with many experiments [10,12,24,25,26].
In particular our expression for the mass flux Q, (66), is in excellent agreement with the data set of Creyssels et al. [10] (see Figure 12). The authors found a scaling of Q with u 2 * − u 2 t , and it can be seen that this is also the dominant term in (66). The same scaling was also found by theoretical studies [19,34,47,48], and recently measured by Ho et al. [11], who showed that this scaling relation is a consequence of the bed being erodible instead of fixed. The splash-entrainment mechanism on erodible beds causes the particle slip velocity V o to be constant with increasing u * , what in turn hinders the average particle velocity V to increase fast with u * . The main increase of the mass flux Q = MV with u * comes instead from the average transported mass per unit soil area M, which scales with u 2 * − u 2 t , as we and other studies e.g. [23] showed. In conclusion, the facts that our model stays on sound physical foundations and that it is in very good agreement with experiments and state of the art simulations make us confident that it for the first time reliably quantifies the feedback of the sand transport on the wind momentum, even on planets different from Earth.
where the parameters are given by r = 0.02, F = 0.96, ǫ ≈ 1s/m, α R = 0.55, α ej = 0.15 g/g f and g f is an effective gravity, which incorporates the effect of cohesion at small particle diameters d, Here ζ is the dimensional cohesion parameter. (A.6) expresses that cohesive forces scale with d 1 in contrast to the gravity force, which scales with d 3 , which means that at small d cohesive forces become dominant. Shao and Lu [49] estimated the magnitude of ζ to be between about 1x10 −4 N/m and 5x10 −4 N/m, we use ζ = 5x10 −4 N/m.

Appendix B. Calculation of the average wind velocity profile and the apparent roughness
In this appendix we show how one can compute the average wind velocity profile and the apparent roughness from an exponentially decaying grain shear stress profile, The average wind velocity profile u(z) can be obtained from Prandtl's mixing length approximation [17,18,20,23,33] as where τ = ρ w u 2 * . The apparent roughness z * o , is the roughness of the asymptotic fluid velocity profileṽ(z) = v(z)| z→∞ , giving In order to integrate (B.2) analytically, we Taylor-expand the square root in (B.2) in the argument a exp(−z/z s ), where a = τ go /τ . It becomes where f j is given by Then (B.2) becomes Within the integral we transform the z-coordinate using such that the integral transforms to where E 1 (x) is called the exponential integral. It can be expressed as and Ein(x) is an analytic function with the series-expansion Using (B.9) the velocity profile finally writes To our knowledge, the sum on the right hand side is not analytically treatable with common methods. However, we found that it can be very well approximated by where both constants, 2.56 and 1.154, are fit constants, ensuring very good agreement between the infinite sum and the approximation. This approximation performs very well over the whole range of a, the relative errors being typically below 1%, as is shown in Figure B1. Now we can finally write     Figure B4 versus the exact solution of the initial boundary value problem (squares), calculated by (B.13) with an accuracy of at least 10 −5 . It can be seen that (B.29) is an excellent approximation of the exact solution of the boundary value problem for small z and (B.21) an excellent approximation for large z. Since both approximations underestimate the analytic solution, the maximum value of u(z) calculated using (B.21) and (B.29), represents an excellent approximation for the whole range of z.

Appendix C. Surface roughness of a quiescent sand bed
The surface roughness z o in the absence of saltation depends on the roughness Reynolds number where k s is the equivalent Nikuradse roughness [52,53]. k s equals the grain diameter d, if the grains of the sand bed are monodisperse, spherical, and very well arranged, meaning that the center point of each particle of the topmost layer is at the same height. However, under more natural conditions k s can be larger, depending on the grain size distribution and the arrangement of the sand bed, and the shape of the grains. A typical value, which is used by engineers, is k s = 3d 84 for water flows in pipes and flumes [53], where d 84 denotes the diameter value which is larger than 84% of the grains of the grain size distribution, or k s = 2d for wind flows [15]. However, since we validate our model with experiments, which used narrowly distributed sand [10,12], we use k s = d. Note that the value of k s does not much influence the final results in most cases. The well known and widely used roughness law is obtained for large roughness Reynolds numbers R ew > 100, which is called the aerodynamic rough regime. On the other hand, in the limit of low roughness Reynolds numbers R ew < 3, the roughness is proportional to the size of the viscous sublayer z o = µ/(9ρ w u * ), (C. 3) which is called the aerodynamic smooth regime. Most natural conditions for aeolian saltation on Earth fall between those regimes, what is called the aerodynamic transitional regime. For instance we obtain R ew ≈ 6 for wind with a shear velocity of u * = 0.4m/s over a typical sand surface with a mean diameter of d = 250µm and k s = d. The behavior of z o as function of R ew was measured by Nikuradse [52] for pipe flows and described by Cheng and Chiew [51] by the following empirical relation z o /k s is plotted in Figure C1 as function of R ew . It describes a function, which has a minimum in the transitional regime at R ew = 9.6 and converges against (C.2) for large R ew and against (C.3) for low R ew . The same behavior was measured by Dong et al. [54] in wind tunnels. Since the variance of z o /k s between the transitional and rough regime is not very large and the measurement errors of z o are large, it is appropriate for saltation models to use a constant value for z o /k s like z o = k s /30 in these regimes. However one cannot neglect the very strong increase of z o for Reynolds numbers below R ew = 3 in the smooth regime. Since we consider very small particle diameters in our model, for which R ew < 3, we use (C.4) and (C.5) to compute z o .