The physics of wind-blown sand and dust

The Saffman lift force is caused by the sharp gradient in the fluid velocity above the particle bed which, due to Bernoulli's principle, creates a lower pressure above the particle than below it (Saffman 1965). Measurements indicate that, at roughness Reynolds numbers (defined in section 2.1.4.1) of ∼10²–10⁴, the lift force is of the order of magnitude of the drag force (Einstein and El-Samni 1949, Chepil 1958, Apperley 1968, Bagnold 1974, Brayshaw et al 1983, Dwivedi et al 2011). On the basis of these measurements, Greeley and Iversen (1985) and Carneiro et al (2011) argued that the lift force plays an important role in initiating saltation. However, at lower roughness Reynolds numbers of ∼1–10, characteristic for saltation in air, the theoretical treatments of Saffman (1965) and McLaughlin (1991) predict a lift force that is substantially smaller than the drag force. Measurements of the lift force at low Reynolds number are difficult to make, but available numerical simulations and measurements are nonetheless consistent with these theories (Loth 2008). The Saffman lift force might thus play a relatively minor role in the initiation of saltation, consistent with the negligible role it plays in determining the trajectories of saltating particles (see section 2.1.2). Further research is required to settle this question.

In addition to the Saffman lift force, some experiments suggest that a lift force due to the pressure deficit at the core of dust devils and other vortices might aid in lifting particles (Greeley et al 2003, Neakrase and Greeley 2010a) (see section 5.3.2).

When granular particles such as dust and sand are brought into contact, they are subjected to several kinds of cohesive interparticle forces, including van der Waals forces, water adsorption forces and electrostatic forces (e.g. Castellanos 2005). Van der Waals forces are themselves a collection of forces that arise predominantly from interatomic and intermolecular interactions, for example from dipole–dipole interactions (Krupp 1967, Muller 1994). For macroscopic objects, the van der Waals forces scale linearly with particle size (Hamaker 1937). Electrostatic forces arise from net electric charges on neighboring particles (Krupp 1967). These charges can be generated through a large variety of mechanisms (Lowell and Rose-Innes 1980, Baytekin et al 2011), including through exchange of ion or electrons between contacting particles (McCarty and Whitesides 2008, Liu and Bard 2008, Forward et al 2009, Kok and Lacks 2009). Finally, water adsorption forces are caused by the condensation of liquid water on particle surfaces, which creates attractive forces both through bonding of adjacent water films and through the formation of capillary bridges between neighboring particles, whose surface tension causes an attraction between these particles (see section 4.1.1.1) (Herminghaus 2005, Nickling and McKenna Neuman 2009). Experiments using atomic force microscopy have found a strong increase of cohesive forces with relative humidity for hydrophilic materials, but not for hydrophobic materials (Jones et al 2002). Mineralogy and surface contaminants affecting hydrophilicity can thus be expected to affect the dependence of cohesive forces on relative humidity.

Although the behavior of the different interparticle forces for aspherical and rough sand and dust is poorly understood (e.g. McKenna Neuman and Sanderson 2008), these forces can be expected to scale with the area of interaction. By combining this assumption with the Hertzian theory of elastic contact (Hertz 1882), Johnson, Kendall and Roberts (1971) developed a theory predicting that the force required to separate two spheres scales linearly with the particle size. This JKR theory has been verified by a range of experiments (e.g. Horn et al 1987, Chaudhury and Whitesides 1992), including experiments that explicitly confirmed the linear dependence of the interparticle force on particle size (Ando and Ino 1998, Heim et al 1999). Therefore, absent a strong correlation of particle morphology or other surface properties with particle size, the interparticle force is expected to scale with the soil particle size (e.g. Shao and Lu 2000).

Iversen and White (1982) used a combination of theory and curve-fitting to wind tunnel measurements (Iversen et al 1976) to derive a semi-empirical expression for the saltation fluid threshold. Their expression is more general than that of Bagnold (equation (2.5)) because it includes the effects of lift and interparticle forces. Specifically, Iversen and White (1982) expressed the dimensionless threshold parameter A_ft in equation (2.5) in a series of semi-empirical equations valid for different ranges of the friction Reynolds number $R_{\ast_{\rm t}}$ (see also Greeley and Iversen 1985):

$$\begin{eqnarray} &&\fl A_{{\rm ft}} =0.2\sqrt {\frac{( {1+0.006/\rho _{\rm p} gD_{\rm p}^{{\rm 2.5}}} )}{1+2.5R_{{\rm \ast_t}}}},\quad {\rm for}\ 0.03\le R_{{\rm \ast_t}} \le 0.3,\nonumber\\ &&\fl A_{{\rm ft}} =0.129\sqrt {\frac{( {1+0.006/\rho _{\rm p} gD_{\rm p}^{{\rm 2.5}}} )}{1.928R_{{\rm \ast_t}}^{{\rm 0.092}} -1}},\nonumber\\ &&\quad{\rm for}\ 0.3\le R_{{\rm \ast_t}} \le 10,\tqs {\rm and}\nonumber\\ &&\fl A_{{\rm ft}} =0.120\sqrt {( {1+0.006/\rho _{\rm p} gD_{\rm p}^{{\rm 2.5}}} )} \nonumber\\ &&\hspace{6pt}\times\left\{{1-0.0858\exp \left[ {-0.0671\left( {R_{{\rm \ast_t}} -10} \right)} \right]} \right\},\nonumber\\ &&\quad{\rm for}\ R_{{\rm \ast_t}} \ge 10, \end{eqnarray} \tag{ 2.6 }$$

where

$$\begin{equation} R_{{\rm \ast_t}} =\rho _{\rm a} u_{{\rm \ast_{ft}}} D_{\rm p} /\mu , \end{equation} \tag{ 2.7 }$$

and where μ is the dynamic viscosity.

Although the semi-empirical expression of Iversen and White (1982) is in good agreement with a large number of wind tunnel measurements for varying conditions (e.g. Iversen et al 1976), Shao and Lu (2000) argued that equation (2.6) can be simplified substantially. In particular, they showed that the dependence of A_ft on the Reynolds number appears limited, such that equation (2.6) can be simplified by assuming that A_ft is independent of the Reynolds number. Shao and Lu (2000) furthermore assumed that the interparticle forces scale with the particle diameter, which is supported by both theory and experiments (section 2.1.1.2). By balancing the entraining forces (aerodynamic drag and lift) against the stabilizing forces (gravity and interparticle forces), Shao and Lu then obtained a relatively straightforward and elegant expression for the saltation fluid threshold

$$\begin{equation} u_{\ast_{\rm ft}} =A_{\rm N} \sqrt {\frac{\rho _{\rm p} -\rho _{\rm a}}{\rho _{\rm a}}gD_{\rm p} +\frac{\gamma}{\rho _{\rm a} D_{\rm p}}} , \end{equation} \tag{ 2.8 }$$

where the dimensionless parameter A_N = 0.111, which is close to the value originally obtained by Bagnold (equation (2.5)). The parameter γ scales the strength of the interparticle forces; for dry, loose dust or sand on Earth, Shao and Lu (2000) recommended a value of γ = 1.65 × 10⁻⁴–5.00 × 10⁻⁴ N m⁻¹. Kok and Renno (2006) obtained γ = 2.9 × 10⁻⁴ N m⁻¹, consistent with this proposed range, by fitting equation (2.8) to the threshold electric field required to lift loose dust and sand particles.

The above expressions of Iversen and White (1982) and Shao and Lu (2000) are compared to measurements of the fluid threshold in figure 5. A critical feature of the threshold curves is the occurrence of a minimum in the fluid threshold for particle sizes of approximately ∼75–100 µm. At this particle size, the interparticle forces are approximately equal in magnitude to the gravitational force, such that the lifting of smaller and larger particles is impeded by increases in the strength of the interparticle and gravitational forces, respectively, relative to the aerodynamic forces. The occurrence of this minimum in the fluid threshold is critical in understanding the physics of dust emission (section 4). Indeed, it shows that saltation is initiated at wind speeds well below that required to aerodynamically lift dust aerosols of a typical size of ∼1–10 µm. As a consequence, dust aerosols are only rarely lifted directly by wind (Loosmore and Hunt 2000). Instead, dust aerosols are predominantly emitted by the impacts of saltating particles on the soil surface (Gillette et al 1974, Shao et al 1993a). However, several plausible dust lifting mechanisms that do not directly involve saltation have been proposed to occur on Mars (see section 4.2).

Zoom In Zoom Out Reset image size

Note that the expressions for the saltation fluid threshold reviewed in this section are for ideal soils: loose, dry sand without roughness elements such as pebbles, rocks or vegetation. As such, the applicability of these expressions is limited to dry sand dunes and beaches. In contrast, the emission of dust aerosols often occurs from soils containing roughness elements or soil moisture. Several corrections are thus required to calculate the fluid threshold for saltation (and thus dust emission) for these soils. The corrections for soil moisture content and for the presence of roughness elements are discussed in sections 4.1.1.1 and 4.1.1.2, respectively.

Most early numerical models of saltation neglected the effect of turbulence on particle trajectories because of the large computational cost (e.g. Ungar and Haff 1987, Anderson and Haff 1991). The steady increase in computational resources has allowed more recent numerical models to explicitly account for the effect of turbulence on particle trajectories (Almeida et al 2006, Kok and Renno 2009a). The results of Kok and Renno (2009a) indicate that, for saltation on Earth, turbulence substantially affects the trajectories of particles smaller than ∼200 µm (figure 7). However, the characteristics of turbulence in the saltation layer remain uncertain, particularly because the effects of saltation on turbulence intensity and the Lagrangian timescale over which the fluctuating component of the fluid velocity is correlated (van Dop et al 1985) have been poorly studied. Indeed, Kok and Renno (2009a) reported discrepancies between the simulated and measured vertical flux profiles of ∼100–200 µm particles, which could indicate inaccuracies in the treatment of turbulence in numerical saltation models. In particular, some wind tunnel measurements indicate that the wakes shed by saltating particles can increase the turbulence intensity (Taniere et al 1997, Nishimura and Hunt 2000, Zhang et al 2008, Li and McKenna Neuman 2012), which is not accounted for in most numerical saltation models.

Zoom In Zoom Out Reset image size

The balance between saltators lost through failure to rebound and saltators gained through splash largely determines the characteristics of steady-state saltation (Kok 2010a). Consequently, the probability that a particle does not rebound is a critical parameter in physically based saltation models. Unfortunately, there have been very few studies of the probability that an impacting saltator fails to rebound, which contrasts with the many published studies on particle splash (e.g. Nalpanis et al 1993, Beladjine et al 2007). The most important exception is the numerical study of Anderson and Haff (1991), which suggested that the probability that a saltating particle will rebound upon impact can be approximated by

$$\begin{equation} P_{{\rm reb}} =B[ {1-\exp ( {-\gamma v_{{\rm imp}}})} ], \end{equation} \tag{ 2.10 }$$

where v_imp is the saltator's impact speed. Mitha et al (1986) determined the parameter B to be ∼0.94 for experiments with 4 mm steel particles, whereas Anderson and Haff (1991) determined a similar value of B ≈ 0.95 for ∼250 µm sand particles. That is, even at large impact speeds not all saltating particles rebound from the surface, because some will dig into the particle bed (Rice et al 1996). The parameter γ has not been experimentally determined, although the numerical experiments of Anderson and Haff (1991) indicate that it is of order 2 s m⁻¹. This value is roughly consistent with the value of γ = 1 s m⁻¹ that Kok and Renno (2009a) deduced from numerically reproducing the impact threshold, which is heavily dependent on equation (2.10) and the splash function. Note that the dimensionality of γ suggests that it could depend on the particle diameter, the gravitational constant, the spring constant and other relevant physical parameters.

If the saltator does rebound, the average restitution coefficient of the collision (i.e. the fraction of the impacting momentum retained by the rebounding particle) has been determined to lie in the range of ∼0.5–0.6 for a loose sand bed by a wide range of numerical and laboratory experiments (Willetts and Rice 1985, 1986, Anderson and Haff 1988, 1991, McEwan and Willetts 1991, Nalpanis et al 1993, Rice et al 1995, 1996, Dong et al 2002, Wang et al 2008, Gordon and Neuman 2011). The mean restitution coefficient appears to be relatively invariant to changes in particle size (Rice et al 1995), but decreases with the impact angle (Beladjine et al 2007, Oger et al 2008). Moreover, recent measurements by Gordon and Neuman (2009) indicate that the restitution coefficient is a declining function of the saltator impact speed, which is consistent with both viscoelastic (Brilliantov et al 1996) and plastic (Lu and Shao 1999) deformation of the soil bed during impact. However, experiments with plastic beads found that the restitution coefficient remains constant or decreases only slightly with impact speed (Rioual et al 2000, Beladjine et al 2007, Oger et al 2008), such that further research is required to resolve this issue.

Measurements indicate that the mean angle from horizontal of the rebounding particle's velocity is around ∼30°–45° (McEwan and Willetts 1991, Nalpanis et al 1993, Rice et al 1995, Dong et al 2002).

In steady-state saltation, the loss of saltating particles to the soil bed through failure to rebound (equation (2.10)) must be balanced by the creation of new saltating particles through splash. Numerical simulations (Anderson and Haff 1988, 1991, Oger et al 2005, 2008), laboratory experiments (Werner 1990, Rioual et al 2000, Beladjine et al 2007) and theory (Kok and Renno 2009a) indicate that the number of splashed particles (N) scales with the impacting momentum. That is,

$$\begin{equation} N\overline {m_{{\rm spl}}} \overline {v_{{\rm spl}}} =\alpha m_{{\rm imp}} v_{{\rm imp}} , \end{equation} \tag{ 2.11 }$$

where $\overline{m_{{\rm spl}}}$ and $\overline{v_{{\rm spl}}}$ are the average mass and speed of the particles splashed by the impacting saltator with mass m_imp. The experiments of Rice et al (1995) indicate that the fraction α of the average impacting momentum spent on splashing surface particles is of the order of ∼15% for a bed of loose sand particles. (Note that the ejection of dust aerosols from the soil, which is impeded by energetic interparticle bonds, probably scales with the impacting energy instead (Shao et al 1993a, Kok and Renno 2009a).)

Since the number of splashed particles N is discrete, equation (2.11) produces two distinct regimes. For collisions with N ∼ 1, an increase in the impacting particle momentum will result in an increase in the momentum of the single splashed particle, such that

$$\begin{equation} \overline {v_{{\rm spl}}} =\alpha v_{{\rm imp}} m_{{\rm imp}} /\overline {m_{{\rm spl}}}\tqs (N \sim 1). \end{equation} \tag{ 2.12 }$$

This linear dependence of the speed of splashed particles on the impacting momentum is supported by the numerical simulations of Anderson and Haff (1988, see figure 9). Conversely, for collisions with N ≫ 1, momentum conservation and equation (2.11) require that the average speed of ejected particles becomes independent of the impacting momentum, as confirmed by the experiments of Werner (1987, 1990), Rioual et al (2000), Beladjine et al (2007) and Oger et al (2008). That is,

$$\begin{equation} \frac{\overline {v_{{\rm spl}}}}{\sqrt {gD_{\rm p}}}=\frac{\alpha}{a},\tqs (N \gg 1), \end{equation} \tag{ 2.13 }$$

where $\sqrt{gD_{\rm p}}$ nondimensionalizes the mean splashed particle speed, and the dimensionless parameter a is a proportionality constant that scales the number of splashed particles (Kok and Renno 2009a):

$$\begin{equation} N=a\frac{m_{{\rm imp}}}{\overline {m_{{\rm spl}}}}\frac{v_{{\rm imp}}}{\sqrt {gD}}. \end{equation} \tag{ 2.14 }$$

Measurements (McEwan and Willetts 1991, Rice et al 1995, 1996) and results from the numerical saltation model COMSALT (Kok and Renno 2009a) indicate that a is ∼0.01–0.03. The simplest way to capture the two limits of equations (2.12) and (2.14) is then as follows:

$$\begin{equation} \frac{\overline {v_{{\rm spl}}}}{\sqrt {gD}}=\frac{\alpha}{a}\left[ {1-\exp \left( {-a\frac{m_{{\rm imp}}}{\overline {m_{{\rm spl}}}}\frac{v_{{\rm imp}}}{\sqrt {gD}}} \right)} \right], \end{equation} \tag{ 2.15 }$$

which produces good agreement with numerical and laboratory experiments (figure 9). The speed of splashed particles is thus generally about an order of magnitude less than that of the impacting particle (Rice et al 1995). Most splashed particles therefore move in low-energy reptation trajectories and quickly settle back to the soil bed. However, some splashed particles do gain enough momentum from the wind to participate in saltation and splash up more particles. This multiplication process produces a rapid increase in the particle concentration upon initiation of saltation (section 2.2).

Zoom In Zoom Out Reset image size

Measurements indicate that the distribution of splashed speeds at a given $\overline{v_{{\rm ej}}}$ follows either an exponential or a lognormal distribution (Mitha et al 1986; Beladjine et al 2007, Ho et al 2012; see figure 5 in Kok and Renno (2009a) for a compilation). Furthermore, laboratory and numerical experiments indicate that the mean angle at which particles are splashed is ∼40°–60° from horizontal (Willetts and Rice 1985, 1986, 1989, Anderson and Haff 1988, 1991, Werner 1990, McEwan and Willetts 1991, Rice et al 1995, 1996, Gordon and Neuman 2011).

The large length scale of the atmospheric boundary layer in which saltation occurs causes the Reynolds number of the flow to be correspondingly large, typically in excess of 10⁶ for both Earth and Mars, such that the flow in the boundary layer is turbulent. Since the horizontal fluid momentum higher up in the boundary layer exceeds that near the surface, eddies in the turbulent flow on average transport horizontal momentum downward through the fluid. Together with the much smaller contribution due to the viscous shearing of neighboring fluid layers, the resulting downward flux of horizontal momentum constitutes the fluid shear stress τ. Because the horizontal fluid momentum is transported downward through the fluid until it is dissipated at the surface, τ is approximately constant with height above the surface for flat and homogeneous surfaces (Kaimal and Finnigan 1994).

The downward transport of fluid momentum through both viscosity and turbulent mixing equals (e.g. Stull 1988)

$$\begin{equation} \tau =\left( {K\rho _{\rm a} +\mu} \right)\frac{\partial \overline {U_x} \left( z \right)}{\partial z}\approx K\rho _{\rm a} \frac{\partial \overline {U_x} \left( z \right)}{\partial z}, \end{equation} \tag{ 2.16 }$$

where we neglect the contribution from viscosity since the turbulent momentum flux exceeds the viscous momentum flux by several orders of magnitude for turbulent flows (e.g. White 2006). The eddy (or turbulent) viscosity K quantifies the transport of momentum in turbulent flows in analogy with the exchange of fluid momentum through viscosity, and $\overline{U_x}(z)$ is the mean horizontal fluid velocity at a height z above the surface. The size of turbulent eddies that transport the fluid momentum is limited by the distance from the surface, such that the eddy viscosity is usually assumed proportional to z (Prandtl 1935, Stull 1988) as well as the shear velocity u_* (defined in equation (2.4)),

$$\begin{equation} K=\kappa u_\ast z, \end{equation} \tag{ 2.17 }$$

where κ ≈ 0.40 is von Kármán's constant. (Note that recent results by Li et al (2010) have called into question the presumed constancy of κ during saltation.) Combining equations (2.16) and (2.17) then yields

$$\begin{equation} \overline {U_{x}} \left( z \right)=\frac{u_\ast}{\kappa}\ln \left( {\frac{z}{z_0}} \right), \end{equation} \tag{ 2.18 }$$

where z₀ is the aerodynamic surface roughness, which denotes the height at which the logarithmic profile of (2.18), when extrapolated to the surface, yields zero wind speed. The value of z₀ for different flow conditions is discussed in more detail below.

Equation (2.18) is referred to as the logarithmic law of the wall and is widely used to relate the shear velocity u_*, and thus the wind shear stress τ, to the vertical profile of the mean horizontal wind speed. However, there are several limitations to its use (Bauer et al 1992, Namikas et al 2003). First, equation (2.18) was derived assuming that the shear stress in the surface layer is constant with height. This is a realistic approximation for flat, homogeneous surfaces, but can be unrealistic for other conditions, such as for surface with non-uniform surface roughness or substantial elevation changes (e.g. Bauer et al 1992). Second, as discussed in the next section, the drag by saltating particles reduces the horizontal momentum flux carried by the wind. Therefore, determining u_* requires measurements of the wind speed above the saltation layer. Third, turbulence causes the instantaneous wind speed to substantially vary over time, such that equation (2.18) is only valid on timescales long enough for the fluid to access all relevant frequencies of the turbulence (Kaimal and Finnigan 1994, van Boxel et al 2004). Because of the large length scale of the atmospheric boundary layer (∼1000 m), this time scale is of the order of ∼10 min (Stull 1988, McEwan and Willetts 1993). And, finally, equation (2.18) is technically only valid when the stability of the atmosphere is neutral (that is, the vertical temperature profile follows the adiabatic lapse rate). Unstable conditions produce a lower wind speed at given values of z and u_*, whereas stable conditions produce a higher wind speed. However, the corrections required to account for the stability of the atmosphere in the wind speed profile of equation (2.18) are small close to the surface (Kaimal and Finnigan 1994), such that this equation can be considered sufficient for most studies of saltation.

An exception to the predominantly turbulent flow in the atmospheric boundary layer occurs near the surface, where surface friction can cause viscous forces to dominate over inertial forces, resulting in laminar flow. For a smooth surface, the thickness of this viscous (or laminar) sublayer is approximately (White 2006)

$$\begin{equation} \delta_{{\rm vis}} \cong \frac{5\mu}{\rho _{\rm a} u_\ast}, \end{equation} \tag{ 2.19 }$$

where μ is the dynamic viscosity. The thickness of the viscous sublayer at the fluid threshold is ∼0.4 mm on Earth and ∼2 mm on Mars. Consequently, sand particles protrude into the viscous sublayer, thereby enhancing turbulent transport by shedding wakes and partially dissipating the viscous sublayer. The thickness of the viscous sublayer thus depends on the size of these roughness elements, as denoted by the roughness Reynolds number (White 2006),

$$\begin{equation} {\it Re}_{\rm r} =\frac{\rho _{\rm a} k_{\rm s} u_\ast}{\mu}, \end{equation} \tag{ 2.20 }$$

where k_s is the Nikuradse roughness (Nikuradse 1933). For a homogeneously arranged bed of monodisperse spherical particles we have that k_s ≈ D_p (Bagnold 1938), whereas for a more realistic irregular surface of mixed sand k_s equals two to five times the median particle size (Thompson and Campbell 1979, Greeley and Iversen 1985, Church et al 1990). For Re_r > ∼60, the turbulent mixing generated by the roughness elements is sufficient to destroy the viscous sublayer and the flow is termed aerodynamically rough (Nikuradse 1933). Consequently, the aerodynamic surface roughness z₀ (defined by equation (2.18)) is related in a straightforward manner to the Nikuradse roughness (White 2006),

$$\begin{equation} z_0 =k_{\rm s} /30\tqs ({\it Re}_{\rm r} > 60). \end{equation} \tag{ 2.21 }$$

In contrast, for Re_r <∼ 4 the roughness elements are too small to substantially perturb the viscous sublayer, and the flow is termed aerodynamically smooth (Nikuradse 1933). In this case, the aerodynamic surface roughness is thus proportional to the thickness of the viscous sublayer (equation (2.19)) and equals (Greeley and Iversen 1985)

$$\begin{equation} z_0 =\frac{\mu}{9\rho _{\rm a} u_\ast}\tqs ({\it Re}_{\rm r} < 4). \end{equation} \tag{ 2.22 }$$

Aeolian saltation on Earth takes place for roughness Reynolds numbers of ∼1–100 and thus usually occurs in the transition zone between the smooth and rough aerodynamic regimes (Dong et al 2001). Since the roughness in the transition regime does not differ much from that in the aerodynamically rough regime (see figure 11 in Nikuradse 1933), most studies have used equation (2.21) to approximate the surface roughness (e.g. Anderson and Haff 1991, Sorensen 1991, 2004, Kok and Renno 2009a, Jenkins et al 2010). In contrast, the much lower air density on Mars causes saltation there to usually take place in the smooth aerodynamic regime, resulting in values of z₀ substantially larger than predicted by equation (2.21). However, most studies of martian saltation have still used this equation for z₀ (e.g. White 1979, Almeida et al 2008, Kok 2010a, 2010b) because, once saltation is initiated, the motion of the saltating particles through the viscous sublayer dramatically increases the vertical mixing of horizontal fluid momentum (e.g. Owen 1964), which likely eliminates the viscous sublayer (White 1979).

In the presence of saltation, the logarithmic wind profile of equation (2.18) is modified by the transfer of momentum between the fluid and saltating particles. Some of the downward horizontal momentum flux in the saltation layer is thus carried by saltating particles, with the total downward momentum flux remaining constant. That is (e.g. Raupach 1991),

$$\begin{equation} \tau =\tau _{\rm a} \left( z \right)+\tau _{\rm p} \left( z \right), \end{equation} \tag{ 2.23 }$$

with the particle momentum flux τ_p given by (Shao 2008)

$$\begin{equation} \tau _{\rm p} (z)=\phi (z)\Delta v_{\rm x} (z), \end{equation} \tag{ 2.24 }$$

and where the fluid momentum flux in the saltation layer τ_a(z) reduces to τ for z above the saltation layer, φ(z) is the mass flux of particles passing the height z in either the upward or downward direction, and Δv_x(z) is the difference between the average velocity of descending and ascending particles. Numerical models find that τ_p decays approximately exponentially from a maximum value near the surface (figure 10). Pähtz et al (2012) derive the dependence of the e-folding length of the particle shear stress on physical parameters such as particle size and the wind shear velocity.

Zoom In Zoom Out Reset image size

In analogy with equation (2.4), the shear velocity within the saltation layer equals

$$\begin{equation} u_{\ast_{\rm a}} (z)=\sqrt {\tau _{\rm a} (z)/\rho _{\rm a}}, \end{equation} \tag{ 2.25 }$$

which similarly reduces to u_* (the shear velocity above the saltation layer) for z above the saltation layer. Combining equation (2.25) with the analogs of equations (2.16) and (2.17) for the saltation layer then yields

$$\begin{equation} \frac{\partial \overline {U_{\rm x}} (z)}{\partial z}=\frac{1}{\kappa z}\sqrt {u_\ast ^2-\tau _{\rm p} \left( z \right)/\rho _{\rm a}}. \end{equation} \tag{ 2.26 }$$

Although equation (2.26) does not have a straightforward analytical solution (Sorensen 2004, Durán and Herrmann 2006a, Pähtz et al 2012), its implementation in numerical saltation models allows a straightforward calculation of the wind profile (see section 2.3.2.4 and, e.g., Werner 1990). Moreover, by assuming a plausible function for the vertical profile of the particle shear stress, such as an exponential profile (figure 10), analytical treatments have been able to derive approximate expressions of the wind profile (Sorensen 2004, Li et al 2004) and the aerodynamic roughness length in saltation (Raupach 1991, Durán and Herrmann 2006a, Pähtz et al 2012; see section 2.3.2.4).

As discussed above, the dominance of splash entrainment in steady-state saltation requires the mean particle speed at the surface to be independent of u_*. Since the mean speed of particles high in the saltation layer can be expected to increase with u_*, the vertical profiles of particle speed must converge towards a common value near the surface. A range of recent wind tunnel measurements (Rasmussen and Sorensen 2008, Creyssels et al 2009, Ho et al 2011), field measurements (Namikas 2003) and numerical saltation models (Kok and Renno 2009a, Durán et al 2011a) confirm both these results (figures 14 and 15). Note that the independence of surface particle speeds with u_* cannot be reproduced by numerical saltation models using Owen's hypothesis (figure 15).

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Although the mean speed of particles at the surface thus remains approximately constant with u_*, the probability distribution of particles speeds at the surface does not. Indeed, increases in u_* produce increases in the wind speed above the Bagnold focus but decreases in wind speed below the Bagnold focus (figure 12). Consequently, the speed of energetic particles moving mostly above the Bagnold focus increases with u_*, whereas the speed of less energetic particles moving mostly below the Bagnold focus decreases (figure 16). In other words, the probability distribution of particle speeds at the surface will broaden. In particular, the tail of very large impact speeds increases strongly with u_*, which is offset by a decrease in the probability of lower particle speeds (figure 16).

Zoom In Zoom Out Reset image size

Since particle trajectories vary widely (e.g. Anderson and Hallet 1986, Namikas 2003), there is no universal measure of the height of the saltation layer. However, since vertical profiles of the horizontal mass flux are commonly measured in wind tunnel and field experiments (e.g. Greeley et al 1996, Namikas 2003, Dong and Qian 2007), convenient ways to define the saltation layer height include the e-folding length in the vertical profile of the horizontal mass flux (e.g. Kawamura 1951, Farrell and Sherman 2006) and the height below which 50% of the mass flux occurs (Kok and Renno 2008). The height of the saltation layer is thus determined by (i) the distribution of speeds with which saltators leave the surface, which sets the vertical concentration profile of saltators, and (ii) the mean horizontal saltator speed at each height. As such, measurements of the saltation layer height can constrain theoretical and numerical predictions of important saltation properties such as the mass flux and particle impact speed.

Early theoretical predictions of the saltation layer height z_salt assumed that the speed with which particles leave the surface scales with u_* (Bagnold 1941, Owen 1964), resulting in

$$\begin{equation} z_{{\rm salt}} =c_{\rm z} u_\ast ^2/2g, \end{equation} \tag{ 2.29 }$$

with c_z of the order of 1. However, as discussed above, advances in theory, numerical saltation models, and measurements have shown that the average speed with which saltating particles are launched from the surface is approximately constant with u_* (figure 15), such that equation (2.29) is likely incorrect (Pähtz et al 2012). And indeed, field measurements of the saltation layer height (Greeley et al 1996, Namikas 2003) find no evidence of the sharp increase with u_* predicted by equation (2.29). Instead, the field measurements of Greeley et al (1996) and Namikas (2003) suggest that the height of the saltation layer remains approximately constant with u_* (figure 17), whereas the recent field measurements of Dong et al (2012) suggest a slight increase in the saltation layer with wind speed (see their figure 7). This latter result is also what is expected from theory and predicted by models: although the particle speed at the surface stays constant, numerical simulations and wind-tunnel measurements find that the particle speed does increase with u_* above the surface (figure 14). This means that the mass flux higher up in the saltation layer increases relative to that in lower layers, producing a slight increase in the saltation layer height with u_* (blue line in figure 17).

Zoom In Zoom Out Reset image size

The saltation mass flux Q, obtained by vertically integrating the horizontal flux of saltating particles, represents the total sand movement in saltation and is thus a critical measure for wind erosion and dune formation studies. Consequently, much effort has been devoted to formulating equations that effectively predict the saltation mass flux (e.g. Bagnold 1941, Kawamura 1951, Owen 1964, Lettau and Lettau 1978, Sorensen 2004).

The saltation mass flux can be derived from the momentum balance in the saltation layer, that is, that the sum of the horizontal momentum fluxes due to particles (τ_p) and the fluid (τ_a) equals the fluid momentum flux at the top of the saltation layer (τ) (equation (2.23)). Applying this relation at the surface yields

$$\begin{equation} \tau _{\rm p} =\rho _{\rm a}({u_\ast ^2-u_{\ast _{{\rm sfc}}}^2}). \end{equation} \tag{ 2.30 }$$

The saltation mass flux is thus generated by the absorption of the excess horizontal wind momentum flux of equation (2.30) by saltating particles. The contribution q of a typical saltating particle hop to the mass flux is the lengthwise mass transport per unit time,

$$\begin{equation} q=mL/t_{{\rm hop}}, \end{equation} \tag{ 2.31 }$$

where m is the particle mass, L is the particle's hop length, and t_hop is the duration of the hop. This hop extracts an amount of momentum p per unit time from the wind,

$$\begin{equation} p=m({v_{{\rm imp}}-v_{{\rm lo}}})/t_{{\rm hop}}, \end{equation} \tag{ 2.32 }$$

where v_imp and v_lo are respectively the horizontal particle speed upon impact and lift-off from the surface. The steady-state saltation mass flux Q is then the available horizontal wind momentum flux τ_p available to accelerate particles (equation (2.30)), multiplied by the mass flux q (equation (2.31)) generated by a unit momentum flux p (equation (2.32)). That is,

$$\begin{equation} Q=\tau_{\rm p} \frac{q}{p}=\rho _{\rm a} ({u_\ast^2-u_{\ast _{{\rm sfc}}}^2} )\frac{L}{\Delta v}, \end{equation} \tag{ 2.33 }$$

where Δv is the average difference between the particle's impact and lift-off speeds.

Different assumptions about the dependence of L, Δv and $u_{\ast_{\rm sfc}}$ on u_* have resulted in different equations relating Q to u_*, as summarized in table 1. In particular, $u_{\ast_{\rm sfc}}$ is usually approximated with $u_{\ast_{\rm it}}$, which, as we noted in section 2.3.1, is incorrect and produces errors when implemented in analytical or numerical saltation models (figures 12 and 15). However, this approximation is more reasonable in this instance because when u_* is close to $u_{\ast_{\rm it}}$, we have that $u_{\ast_{\rm it}} \approx u_{\ast_{\rm sfc}}$ (figure 13), whereas when $u_\ast \gg u_{\ast_{\rm it}}$, both $u_\ast^2-u_{\ast_{\rm it}} ^2$ and $u_\ast ^2-u_{\ast_{\rm sfc}}^2$ approximate $u_\ast ^2$.

Table 1. List of the most commonly used saltation mass flux relations.

Mass flux equation	Comments	Study
$\dsty Q_{{\rm Bagnold}} =C_{\rm B} \sqrt {\frac{D_{\rm p}}{D_{{\rm 250}}}} \frac{\rho _{\rm a}}{g}u_\ast ^3$	CB = 1.5, 1.8 or 2.8 for uniform, naturally graded and poorly sorted sand, respectively.	Bagnold (1941)
$\dsty Q_{{\rm Kawamura}} =C_{\rm K} \frac{\rho _{\rm a}}{g} u_\ast ^3\left( {1-\frac{u_{\ast _{{\rm it}}}^2}{u_\ast ^2}} \right)\left( {1+\frac{u_{\ast _{{\rm it}}}}{u_\ast}} \right)$	CK = 2.78 (Kawamura 1951) or 2.61 (White 1979). The origin of this relation is often confused to be White (1979); see Namikas and Sherman (1997).	Kawamura (1951)
$\dsty Q_{{\rm Owen}} =\frac{\rho _{\rm a}}{g}u_\ast ^3\left( {0.25+\frac{v_{\rm t}}{3u_\ast}} \right)\left( {1-\frac{u_{\ast _{{\rm it}}}^2}{u_\ast ^2}} \right)$	vt is a particle's terminal fall speed.	Owen (1964)
$\dsty Q_{{\rm Lettau}} =C_{\rm L} \sqrt {\frac{D_{\rm p}}{D_{{\rm 250}}}} \frac{\rho _{\rm a}}{{g}}u_\ast ^3\left( {1-u_{\ast _{{\rm it}}} /u_\ast} \right)$	CL = 6.7	Lettau and Lettau (1978)
$\dsty Q_{{\rm UH}} =C_{{\rm UH}} \rho _{\rm a} \sqrt {\frac{D_{\rm p}}{{g}}} u_\ast ^2\left( {1-\frac{u_{\ast _{{\rm sfc}}}^2}{u_\ast ^2}} \right)$	Ungar and Haff (1987) did not estimate a value of CUH.	Ungar and Haff (1987)
$\dsty Q_{{\rm Sorensen}} =\frac{\rho _{\rm a}}{{g}}u_\ast ^3( {1-u_{\ast _{{\rm it}}}^2 /u_\ast ^2} )( {\alpha +\gamma u_{\ast_{\rm it}} /u_\ast +\beta u_{\ast _{{\rm it}}}^2 /u_\ast ^2} )$	α, β and γ are parameters that characterize the dimensions of a typical saltation hop.	Sorensen (2004)
$\dsty Q_{{\rm DK}} =C_{{\rm DK}} \frac{\rho _{\rm a}}{g}u_{\ast_{\rm it}} u_\ast ^2\left( {1-\frac{u_{\ast_{\rm it}} ^2}{u_\ast ^2}} \right)$	CDK ≈ 5	Proposed here and in Durán et al (2011a)

D₂₅₀ is a reference diameter of 250 µm.

Another assumption commonly made to simplify equation (2.33) is that the particle speed in saltation scales with u_*, resulting in L/Δv ∝ u_* (Bagnold 1941). Using this assumption results in a scaling of Q with $u_\ast^3$, as for example proposed by Bagnold (1941), Kawamura (1951) and Owen (1964). However, more recent studies have questioned the linear scaling of particle speeds with u_* and thus the scaling of Q with $u_\ast^3$ (Ungar and Haff 1987, Andreotti 2004, Kok 2010a, 2010b, Durán et al 2011a, Ho et al 2011, Pähtz et al 2012). As discussed in sections 2.3.1 and 2.3.2.1, these studies are supported by recent experimental, numerical and theoretical results showing that particle speeds in steady state, transport limited saltation are not proportional to u_* (e.g. figures 14 and 15). Instead, the mean particle speed at the surface is independent of u_*, whereas the mean particle speed above the surface increases gradually with u_* (as well as with height). The speed with which particles are launched from the surface (and thus Δv) is thus independent of u_*, and the saltation hop length is only a weak function of u_* (see figure 15 and Namikas 2003, Rasmussen and Sorensen 2008, Kok 2010a, Ho et al 2011).

We can obtain a more physically realistic relation for the mass flux from equation (2.33) by neglecting the weak scaling of L with u_* (e.g. Namikas 2003) and using that particle speeds scale with wind speed in the saltation layer, which in turn scales with $u_{\ast_{\rm it}}$. This yields $L/\Delta v=C_{{\rm DK}} u_{\ast_{\rm it}} /g$, such that the mass flux can be approximated by (Kok and Renno 2007, Durán et al 2011a)

$$\begin{equation} Q_{{\rm DK}} =C_{{\rm DK}} \frac{\rho _{\rm a}}{g}u_{\ast_{\rm it}} ( {u_\ast ^2-u_{\ast _{{\rm it}}}^2} ). \end{equation} \tag{ 2.34 }$$

Experiments and numerical simulations suggest that, for ∼250 µm sand, we have that L ≈ 0.1 m (Namikas 2003), Δv ≈ 1 m s⁻¹ (figure 15), and $u_{*_{\rm it}} \approx 0.2$ (e.g. Bagnold 1937), such that the parameter C_DK ≈ 5. We now use equation (2.34) to non-dimensionalize the mass flux, which provides a convenient way to compare different mass flux relations. That is, we take the dimensionless mass flux as

$$\begin{equation} \hat {Q}=\frac{gQ}{\rho _{\rm a} u_{\ast_{\rm it}} ( {u_\ast ^2-u_{\ast_{\rm it}} ^2} )}. \end{equation} \tag{ 2.35 }$$

Figure 18 shows the dimensionless mass flux predicted by the relations in table 1. Although the scatter in the experimental data is substantial, most mass flux data sets seem to support the scaling of equations (2.34) and (2.35).

Zoom In Zoom Out Reset image size

The extraction of wind momentum by saltating particles produces a steady state wind profile that accelerates saltating particles to an impact speed that, on average, results in a single particle leaving the soil bed for each particle impacting it (section 2.3.1). However, as u_* increases, the wind speed higher up in the saltation layer will increase as well (figure 12), thereby increasing the saltator impact speed. But since the mean saltator impact speed must remain constant with u_* (figure 15), this increase of the wind speed higher in the saltation layer must be compensated by a decrease of the wind speed lower in the saltation layer. Consequently, when wind profiles at different values of u_* are plotted on the same graph, the wind profiles intersect in the Bagnold focus at a height of ∼1 cm above the surface for saltation over loose sand (Bagnold 1936). The wind speed below the focus, as well as the wind speed gradient and thus $u_{\ast_{\rm a}}$, actually decreases with u_* (figures 12, 13 and 19).

Zoom In Zoom Out Reset image size

The absorption of wind momentum by saltating particles acts as roughness elements that increase the aerodynamic roughness length sensed by the flow above the saltation layer (Owen 1964). The wind speed above the saltation layer is thus determined by this increase of the aerodynamic roughness length (Pähtz et al 2012). Several models have been proposed to relate z_0s, the aerodynamic roughness length in saltation, to u_*. These include the empirical Charnock (1955) and modified Charnock (Sherman 1992, Sherman and Farrell 2008) models, and the Raupach (1991) model, which is based on Owen's hypothesis (equation (2.28)). More recent physically based models of the aerodynamic roughness length include Durán and Herrmann (2006a) and the Pähtz et al (2012) model, which uses the balances of momentum and energy in the saltation layer to derive a physically based expression for z_0s. The various relations are listed in table 2, and comparisons of some these relations with measurements are given in Sherman and Farrell (2008).

Table 2. Relations for the aerodynamic roughness length in saltation. The relations of Durán and Herrmann (2006a, see their equations (20)–(24)) and Pähtz et al (2012), see their equations (67a)–(67e) are not included here due to their complexity.

Relation	Comments	Study
$z_{{\rm 0S}} =C_{\rm c} u_\ast ^2/g$	Cc ≈ 0.010 for wind tunnels and ≈0.085 for natural saltation (Sherman and Farrell 2008)	Charnock (1955)
$\dsty z_{{\rm 0S}} =\left( {\frac{Au_\ast ^2}{2g}} \right)^{1-r}z_0^r$	$r=u_{\ast_{\rm it}}/u_\ast$; A ≈ 0.2–0.4 (Raupach 1991, Gillette et al 1998); a similar relationship was derived by Owen (1964).	Raupach (1991)
$\dsty z_{{\rm 0S}} =z_0 +C_{\rm m} \frac{\left( {u_\ast -u_{\ast _{{\rm it}}}} \right)^2}{g}$	Cc ≈ 0.012 for wind tunnels and ≈0.13 for natural saltation (Sherman and Farrell 2008)	Sherman (1992)

As discussed in section 2.3.1, after saltation has been initiated by the aerodynamic entrainment of surface particles, subsequent lifting of surface particles occurs predominantly through splash. The fluid threshold therefore has little influence on the physical properties of steady-state saltation. In particular, the size distribution of saltating particles is not determined by the dependence of the fluid threshold on particle size, as many previous studies have assumed (e.g. Marticorena and Bergametti 1995, Shao et al 1996), but is rather determined by (i) the probability that a particle of a given size is splashed into the fluid stream, and (ii) the time the particle spends in saltation before settling back to the soil bed. The experiments of Rice et al (1995) indicate that particles with a larger cross-sectional area have a correspondingly higher change of being splashed by an impacting saltator. Conversely, larger particle have short lifetimes because their high inertia causes them to attain lower speeds when splashed (equation (2.15)), leading them to quickly settle back to the soil surface (Kok and Renno 2009a). These two competing effects appear to result in a size distribution of saltating particles that is similar to the soil size distribution (figure 20).

Zoom In Zoom Out Reset image size

The size distribution of saltating particles is thus primarily determined by the speed (or momentum) of impacting saltators. Since the impact speed of saltators remains constant with u_* (figure 15), the size distribution of saltating particles should also be relatively invariant to changes in u_*, as indeed indicated by field measurements (figure 20).

A foremost open question is how commonly saltation occurs on Mars. Sporadic lander measurements suggest that winds in the light martian atmosphere rarely exceed the saltation fluid threshold (Zurek et al 1992, Sullivan et al 2000; Holstein-Rathlou et al 2010), which is consistent with simulations of both mesoscale (Fenton et al 2005, Chojnacki et al 2011) and global circulation models (Haberle et al 2003). These findings are difficult to reconcile with recent reports of active sand transport at many locations on Mars (Fenton et al 2005, Fenton 2006, Greeley et al 2006, Sullivan et al 2008, Bourke et al 2008, Geissler et al 2010, Silvestro et al 2010, 2011, Chojnacki et al 2011, Hansen et al 2011, Bridges et al 2012b, Kok 2012), which in at least some cases results in sand fluxes that are similar to those on slow-moving dunes on Earth (Bridges et al 2012a).

Part of the explanation for these seemingly contradictory findings might be that small-scale topography and convection generate strong localized winds, which cannot be accurately simulated by the coarse resolution of martian atmospheric circulation models (Fenton and Michaels 2010). Moreover, once such strong localized winds initiate saltation, recent results indicate that it can be sustained by wind speeds reduced by up to a factor of ∼10 (Kok 2010a, 2010b). That is, the martian impact threshold is only ∼10% of the fluid threshold (see figure 21(b)), thus allowing saltation to occur at much lower wind speeds than previously thought possible. The ratio of impact to fluid thresholds on Mars is lower than on Earth (table 3) because the lower gravity and air density allows particles to travel higher and longer trajectories on Mars, causing them to be accelerated by wind for a longer duration during a single hop than on Earth (Kok 2010a). This effect combines with the lower atmospheric density on Mars to produce an impact threshold that is comparable to that on Earth (figure 21). However, the lower atmospheric density causes the martian fluid threshold to be an order of magnitude larger than on Earth (equation (2.5)). Consequently, the ratio of the impact to fluid thresholds is much smaller on Mars than it is on Earth (table 3).

The trajectories of saltating particles on Mars are largely determined by the typical speed with which they are launched from the surface. Interestingly, the mean saltator launch speed during steady-state saltation over a bed of loose sand is predicted to be very similar on Mars and Earth (Kok 2010a, 2010b). The reason for this is that the efficiency of the splashing process is essentially independent of the air density (Durán et al 2011a) and only weakly dependent on particle size and gravity. However, the lower gravity and vertical air drag cause martian saltators to undergo trajectories that are almost an order of magnitude longer and higher than those on Earth (see table 3 and figure 3(b) in Kok 2010b). These results contrast with the substantially larger saltation trajectory estimates of Almeida et al (2008), who did not explicitly account for splash. Since steady state on Mars is probably determined by splash entrainment, neglecting splash can yield unphysical results, such as an increase in the saltator impact speed with u_* (e.g. figure 5 in Almeida et al 2008). The saltation trajectory estimates of Almeida et al (2008) can nonetheless be considered plausible for supply limited saltation (Nickling and McKenna Neuman 2009), for which saltator impact speeds generally do increase with u_* (Houser and Nickling 2001, Ho et al 2011).

Rubin and Hunter (1987) designed an experiment to simulate aeolian ripple formation under bimodal winds. A sand-covered turntable subject to a nearly steady wind of constant direction was periodically rotated by a divergence angle θ_w. The table was held at one position for a time T_w1, then turned to the other position, where it remained a time T_w2, then turned back to the first position, where the process was repeated in a cyclic manner. Ripples always oriented transversely (longitudinally) to the net transport trend when θ_w was smaller (larger) than θ_L = 90°. A 'mixed state' was observed for θ_w within the range 90° < θ_w < 112°. 'Oblique bedforms' occurred for obtuse divergence angle when the transport ratio R = T_w1/T_w2 differed from unity (Rubin and Hunter 1987). The experiments showed that bedform alignment is such that it yields the maximum gross bedform-normal transport—which is defined as the amount of sand transported normal to the bedform axis (Rubin and Hunter 1987). Gross transport means that transport from opposite directions do not cancel, instead both contribute to build bedform surface. The prediction θ_L = 90° has been confirmed in the field (McKee and Tibbitts 1964, Tsoar 1983, 1989, Rubin and Hunter 1985, Bristow et al 2007, Rubin et al 2008, Dong et al 2010), in water-tank experiments (Rubin and Ikeda 1990, Reffet et al 2010), and in numerical simulations of aeolian dunes (Werner 1995, Reffet et al 2010, Parteli et al 2009, see sections 3.2.1.4 and 3.3.2.2).

Very large dunes commonly display superimposed dunes of a different morphology—such dune compounds are said to be complex dunes (McKee 1979). Examples of complex dunes are the giant longitudinal dunes in Namibia, which host smaller superimposed transverse bedforms migrating obliquely to the main dune trend (Lancaster 1988, Andreotti et al 2009). The origin of complex dunes is that bedforms of different sizes respond to wind regimes associated with different time scales (see section 3.1.6.1).

Some natural agents can serve as sand stabilizers, thus greatly affecting bedform morphology. For instance, in vegetated areas U-shaped parabolic dunes constitute the predominant dune morphology (Tsoar and Blumberg 2002; section 3.2.3.1). Straight-crested, downwind-elongating linear dunes (where 'linear' refers to dunes that are much longer than they are wide; see e.g. Rubin and Hesp 2009) can also form due to the action of vegetation, though the formative processes of these dunes are still poorly understood (Pye and Tsoar 1990). Some authors hypothesized that sand induration or cementation by salts or moisture in the dune could also lead to linear dunes elongating longitudinally to a single wind direction (Schatz et al 2006, Rubin and Hesp 2009; see section 3.3.2.3). Furthermore, the presence of an exposed water table may affect local rates of erosion and deposition and thus play a relevant role for dune morphology and mobility (Kocurek et al 1992, Luna et al 2012).

The best-known type of dune is the barchan, which is also of particular interest due to its high mobility: barchans can migrate between 30 and 100 m in a single year (Bagnold 1941, Pye and Tsoar 1990). Many authors reported linear relations between barchan's width (W), length (L) and height (H) and also a scaling of dune migration velocity (v) with the inverse of dune size (Bagnold 1941, Finkel 1959, Long and Sharp 1964, Hastenrath 1967, Embabi and Ashour 1993, Hesp and Hastings 1998, Sauermann et al 2000, Hersen et al 2004, Elbelrhiti et al 2005).

The migration velocity of a barchan dune can be approximately calculated for a dune in steady state by applying mass conservation at the dune's crest (Bagnold 1941). The dune migration velocity scales as v_m ∼ cQ₀/W, where c ≈ 50 (Durán et al 2010). Since the barchan's volume V scales with W³ with a constant of proportionality of about 0.017 (Durán et al 2010), one can write $v_{\rm m}\approx {\it kQ}_{0}/V^{1/3}$, with k ≈ 12.8. The turnover time (T_m ∼ W/v_m) of a dune can be understood as the time scale needed for a sand grain, once buried at the lee, to reappear at the windward dune's foot (Allen 1974, Lancaster 1988, Andreotti et al 2002a). T_m follows the scaling relation,

$$\begin{equation} T_{\rm m} \sim {\it aW}^{2}/Q_0, \end{equation} \tag{ 3.3 }$$

where Q₀ is the bulk sand flux and a ≈ 0.02 from field observations (Hersen et al 2004). So the smaller the dune, the shorter its turnover time, and thus the faster it readapts to a change in flow conditions.

Indeed, the coexistence of dunes of different types associated with distinct sizes in an environment of complex wind regime (e.g. small transverse dunes occuring superimposed on or nearby much larger star or longitudinal dunes; see Lancaster 1988) can be understood on the basis of equation (3.3): sufficiently small bedforms can orient transversely to average wind directions that prevail for a short time scale T_w comparable to the turnover time of these small dunes (10–100 years); in contrast, the shape of much larger dunes provides a proxy for changes in wind regime occurring within time scales similar to their much longer turnover time (1000–100 000 years). Consequently, the aforementioned scaling relations have crucial implication for understanding the different shapes of dunes occurring in Earth and Mars deserts and their connection to local variations in wind regimes occurring at a range of time scales.

As explained previously, dune dynamics respond to variations in wind regime, availability of mobile sediments and different factors that are influenced by local climate, such as vegetation cover. The shape of mobile dunes prevailing in a given field may change over time in response to changes in the local environmental conditions, with new dunes emerging superimposed on more ancient bedforms (Levin et al 2007). Indeed, it is possible to infer information about past climates from the shape of fossil dunes (Lancaster et al 2002). Also, knowledge of how dune dynamics react to climatic changes has potential application to predicting the remobilization of currently inactive sand systems due to global warming (Thomas et al 2005), or to infer information about the history of the martian climate (Fenton and Hayward 2010, Gardin et al 2012).

Changes in wind regimes not only affect the shape of dunes but can also serve as a control mechanism for dune size in a dune field. Large barchans display in general small superimposed dunes ('surface-wave instabilities') induced by variations in wind trend (Elbelrhiti et al 2005). When subject to a storm wind blowing from a direction that makes a small angle with the primary wind, the dune's surface is unstable. On the windward side, small transverse bedforms with size of the order of the minimal dune occur, which then migrate along the dune surface finally escaping through the limbs. This mechanism reduces the mass of large dunes, thus hindering the formation of a single giant dune through merging of all dunes in the field (Elbelrhiti et al 2005).

A further factor relevant for dune size selection is the occurrence of collisions. Due to the scaling of dune velocity with the inverse of dune width, if a small barchan is migrating behind a larger one, a collision will occur. Collisions between barchans in a corridor may generate small 'baby'-barchans escaping through the limbs of large dunes (Schwämmle and Herrmann 2003, Endo et al 2004, Hersen 2005, Durán et al 2005) and could therefore also contribute to regulating the dune size. The respective roles of dune collisions and wind trend fluctuations for the stability of dune fields are currently a matter of substantial debate (Hersen and Douady 2005, Elbelrhiti et al 2008, Durán et al 2009, Diniega et al 2010, Elbelrhiti and Douady 2011, Durán et al 2011b).

Excellent reviews on field investigations of dunes can also be found elsewhere (Pye and Tsoar 1990, Lancaster 1995, Livingstone et al 2007, among others). Field observations provide the basis for conceptual models, as for example the model by Tsoar and Blumberg (2002) for parabolic dunes (section 3.2.3.1) or the model by Lancaster (1989a) for star dunes. The latter author suggested that a dune shaped by two opposite wind directions (i.e. a transverse dune of reversing type) may transform into a star dune with four arms owing to leeside 'secondary flow patterns', which lead to along-axis sand transport from the sides towards the center of the dune. First, the secondary winds build two longitudinal arms, one at the windward side and the other at the lee, whereupon the growth of the emerging arms is accentuated by a third oblique wind direction. Conceptual models have also been proposed for the common occurrence of asymmetric barchan shapes, where one of the barchan limbs extends downwind, eventually developing into a seif dune (Bagnold 1941, Tsoar 1983, Bourke 2010). Experiments or numerical simulations, discussed next, might be employed to test such models.

Water tank experiments provide an alternative approach to study the long time scale processes of dune dynamics (Rubin and Ikeda 1990, Hersen et al 2002, Endo et al 2004). The major difference between subaqueous and aeolian dunes lies in their scale. Owing to the difference between the fluid densities of water and air (see equation (3.1)), underwater dunes are almost 1000 times smaller than their aeolian counterparts. Hersen et al (2002) designed an experimental set-up to produce unimodal water streams strong enough to put glass beads, with size and density similar to those of desert sand grains, into saltation. The experiments produced subaqueous barchan-like bedforms displaying nearly the same morphological relations as those of aeolian dunes, thus suggesting that dune formation should not be strictly dependent upon the details of the transport mechanism (Hersen et al 2002). By using water tank experiments, it has been possible to reproduce in the laboratory different dune morphologies appearing under varying flow conditions (Rubin and Ikeda 1990, Hersen 2005, Reffet et al 2010) as well as to investigate in a systematic manner the collision between barchans (Endo et al 2004) and the formation of a barchan belt emerging from a sand source (Katsuki et al 2005a, 2005b). Water tank experiments have been also employed as a means to test hypotheses on the origin of some extraterrestrial dune shapes (Taniguchi and Endo 2007; see also section 3.3.2.2).

A surface with dunes introduces a perturbation in the shear stress profile over a flat ground. Computational fluid dynamics (CFD) simulations have been employed to solve the Navier–Stokes equations of the turbulent wind over barchans and transverse dunes (Wiggs 2001, Parsons et al 2004, Herrmann et al 2005, Schatz and Herrmann 2006, Livingstone et al 2007). Such an approach has been helpful in elucidating the patterns of recirculating flow at the dune lee (figure 24) as well as the modification of the flow due to the presence of an array of many dunes. However, such calculations are too computationally expensive to be utilized in the framework of a continuum three-dimensional dune model.

A more suitable, less computationally expensive approach consists of using an analytical model to calculate the average turbulent wind shear stress field over the sand terrain. Such a model has been developed in the course of several years (Jackson and Hunt 1975, Weng et al 1991), and has been widely used in dune models (Stam 1997, Kroy et al 2002, Andreotti et al 2002b). It computes the Fourier-transformed components of the shear stress perturbation $(\vec{\hat{\tau}})$ in the directions longitudinal and transverse to the wind, whereupon the (perturbed) shear stress τ can be calculated,

$$\begin{equation} \vec{\tau}=|\vec{\tau}_0|(\vec{\tau_0}/|\vec{\tau}_0|+\vec{\hat{\tau}}), \end{equation} \tag{ 3.5 }$$

where $\vec{\tau}_0$ is the undisturbed shear stress over a flat ground. Note that the shear stress formulation is a linear solution of the Reynolds averaged Navier–Stokes equations. A consequence of this is that the shear stress perturbation is scale invariant: it only depends on the shape of the dune and not its size.

This analytical model works only for smooth heaps and gentle slopes. It does not apply for a dune with a slip face, since it does not account for flow separation occurring at sharp edges and steep slopes. Zeman and Jensen (1988) therefore suggested a heuristic approach by introducing a separation bubble at the dune lee, which comprises the area of recirculating flow. For the purpose of calculating the flow at the windward side of the dune, the wind model can be solved for the envelope,

$$\begin{equation} h_{s}(x,y) = \max\{h(x,y),s(x,y)\}, \end{equation} \tag{ 3.6 }$$

comprising the separation streamlines (s(x, y)) and the dune surface (h(x, y)) (figure 27). Then, at the lee, the shear stress within the separation bubble, where the resulting net transport essentially vanishes, can be set as zero.

Zoom In Zoom Out Reset image size

In addition to the shear stress distribution, a continuum dune model requires an equation for the mass flux, q(x, y), which is then combined with equation (3.4) to give the bed evolution. At this stage, many authors (Wippermann and Gross 1986, Zeman and Jensen 1988, Stam 1997) attempted to use an equation for the saturated flux (q_s), such as the one by Lettau and Lettau (1978), see section 2.3.2.3. However, the flux upwind of a barchan moving on the bedrock is well below the saturated value (Fryberger et al 1984, Wiggs et al 1996). The simulation of a barchan using an equation for saturated flux therefore leads to an unphysical result: the decrease in wind strength close to the dune at the upwind causes deposition at the windward foot (which is not observed in the field; Wiggs et al 1996). The dune's foot then stretches and the dune flattens. This problem was illustrated by van Dijk et al (1999), who added to the saturated flux equation of Lettau and Lettau (1969) a phenomenological 'adaptation length' to account for the distance over which sediment transport adapts to a new equilibrium condition due to a spatial change in transport conditions.

The continuum model by Sauermann et al (2001) accounts for a physical modeling of the saturation transients of the flux. The cloud of saltating particles is regarded as a thin fluid-like layer moving on top of the immobile bed. The flux is described by calculating the average density and velocity of grains jumping with a 'mean saltation length', whereas it is not distinguished between reptating and saltating particles. The spatial evolution of the average height-integrated mass flux over the sand bed is computed by solving the following equation, i.e. the well-known logistic equation for population growth,

$$\begin{eqnarray} &&\vec {\nabla }\cdot \vec {q}=\left[ {1-{\left| {\vec {q}} \right|}/{q_{\rm s} }} \right]\cdot {\left| {\vec {q}} \right|}/{L_{{\rm sat}}}, \end{eqnarray} \tag{ 3.7 }$$

where the saturated flux q_s and the saturation length L_sat, which respectively scale with [τ − τ_t] and [τ/τ_t − 1], encode the information on grain diameter, gravity, fluid viscosity and densities of fluid and particles (Sauermann et al 2001, Durán and Herrmann 2006a, Parteli and Herrmann 2007b). It was shown that the minimal dune width W_min in the simulations is around 12L_sat (Parteli et al 2007a).

So the scaling of the saturation length in the model of Sauermann et al (2001) includes a dependence on the wind shear velocity which is not included in the scaling L_sat ≈ 2l_drag (see section 3.1.1). Experiments show that the distance measured from the edge of a flat sand patch up to the horizontal position downwind where the flux is nearly saturated scales with the shear stress as predicted in the model of Sauermann et al (2001) for L_sat (Andreotti et al 2010). However, L_sat appears to be nearly independent of τ when q is close to q_s. Hersen et al (2002) proposed a simplified version of equation (3.7), where the RHS of this equation is replaced by (q_s − q)/L_sat, with L_sat independent of the shear stress. There is an ongoing debate about the modeling of the saturation length, and so we refer to ongoing research papers for a more detailed discussion (Andreotti and Claudin 2007, Parteli et al 2007b, Fourrière et al 2010).

The continuum model then computes the bed evolution by inserting equation (3.7) into the mass conservation equation (equation (3.4)) (Kroy et al 2002). Wherever the local slope exceeds the angle of repose of the sand (θ_r), the surface must relax through avalanches in the direction of the steepest descent. Avalanches can be considered to occur instantaneously, since their time scale is much smaller than the time scale of dune evolution. Durán et al (2010) presented a model of relatively fast convergence for computing avalanches at the dune lee.

The continuum dune model by Kroy et al (2002) was applied successfully by Sauermann et al (2003) to predict the average shear stress and flux profiles over the symmetry axis of a real barchan dune in northeastern Brazil. It also reproduced quantitatively the three-dimensional shape of a barchan in Morocco (figure 28).

Zoom In Zoom Out Reset image size

Also, important insights into the formation and dynamics of dunes and dune fields were gained using the continuum model.

In particular, the simulations showed that collisions between dunes lead to different dynamics depending on the volume ratio of the small and big dunes, r = V_small/V_big. The dunes merge if r < 7%. If r > 25%, then sand from the larger dune downwind is caught in the wake of the smaller dune behind, which then increases in size, eventually becoming larger (and therefore slower) than its downwind companion. Because the dynamics resembles a small dune crossing over the larger one (Besler 1997), it has been ill-named 'solitary wave-like' behavior of dunes (Schwämmle and Herrmann 2003, Andreotti et al 2002b, Durán et al 2005), whereas the simulation showed that in reality what happens is that large and small dunes simply interchange their roles through the net release of sand from the large dune to the small one (Schwämmle and Herrmann 2003). For r between 14% and 25%, the collision results in the ejection of small 'baby barchans' from the limbs of the larger dune downwind, whereas a further dune is realeased from the central part for 7% < r < 14% (Schwämmle and Herrmann 2003, Durán et al 2005). Indeed, the ejection of small dunes resulting from collisions between dunes has implication for maintaining an equilibrium size distribution in barchan corridors (Durán et al 2009, 2011b, see discussion in sections 3.1.4 and 3.1.6.3).

Furthermore, by using an extension of the model to study dune formation under bimodal winds (Parteli et al 2009) it was possible to reproduce many exotic dune shapes that occur on Mars (section 3.3.2.2). The simulations could also help to understand the meandering of longitudinal seif dunes—which depends, as shown in figure 29, on the duration (T_w) of each of both prevailing winds relative to the dune turnover time (T_m).

Zoom In Zoom Out Reset image size

Moreover, application of the model to investigate the genesis and evolution of barchan dune fields from a sand source has yielded important insights into the origin of dunes under realistic conditions (Durán et al 2010). It was shown that, when a flat sand sheet is subjected to a saturated sand influx, it becomes a source for transverse dunes which emerge from sand-wave instabilities (see section 3.1.1) on the sand surface. The transverse dunes become unstable after reaching the bedrock (Reffet et al 2010, Parteli et al 2011, Niiya et al 2012, Melo et al 2012) and decay into chains of barchans (figure 30). The genesis of barchan dune fields has also been studied with inclusion of vegetation (Luna et al 2009, 2011, see also section 3.2.3) and an exposed water table (Luna et al 2012).

Zoom In Zoom Out Reset image size

While the continuum model can potentially serve as a helpful tool in the quantitative assessment of planetary dunes, some important limitations must be emphasized. One major issue concerns the modeling of the saturation length, L_sat. The saturation length of the dune model does account for inertia of the transported sediments, which limits its application to small values of wind shear stress (Parteli et al 2007b). Indeed, due to the scaling of L_sat with (τ − τ_t)⁻¹, a rapid increase of the minimal dune size with decreasing values of wind shear stress close to τ_t is predicted. While measurements confirming this prediction are scarce (Andreotti et al 2010), the physical processes leading to flux saturation are currently the matter of substantial debate (Andreotti and Claudin 2007, Parteli et al 2007b, Durán et al 2011b). Moreover, the dependence of the threshold friction speed of the wind on the local slope is not accounted for in the model (Andreotti and Claudin 2007, Parteli et al 2007b). One further issue concerns the model for the separation bubble (figure 27). For isolated dunes, setting the wind shear at the lee to zero might have minor consequences for the barchan shape (Kroy et al 2002). However, complex three-dimensional flow patterns at the lee could be critical for the dynamics of closely spaced dunes (Parteli et al 2006) or the elongation of seif dunes (Bristow et al 2000) and should be considered in future modeling.

In comparison to morphodynamic modeling CA algorithms may offer simpler and computationally less expensive numerical tools for simulating dune dynamics. One drawback of CA models is the lack of a physical modeling of the sand flux and the wind profile, which is essential for the quantitative assessment of dunes. Length and time scales are not explicitly established in CA transport rules (Werner 1995), and it is also not possible to perform a stability analysis of dunes using this type of modeling (Narteau et al 2009). Recently, attempts were made to develop more realistic CA models through including quantitative descriptions of the average turbulent wind profile or saltation flux. For example, the model by Zheng et al (2009) considered a wind friction velocity (u_*) that increased linearly with the surface height, while in the model by Pelletier (2009), Werner's transport algorithm was coupled to a generalized version of the boundary layer flow model by Jackson and Hunt (1975). Narteau et al (2009) introduced a three-dimensional model in which the Navier–Stokes equations for the turbulent fluid were solved using a lattice gas cellular automaton model (LGCA). This model was coupled to a probabilistic sand transport model that considered the erosion rate proportional to the excess shear stress, τ–τ_t, thus accounting for the exchange of momentum between particles and fluid that lead to flux saturation. The simulations produced secondary (superimposed) bedforms on the surface of large dunes (Narteau et al 2009, Zhang et al 2010), thus indicating the minimal dune size in units of elementary cells. Therefore, it was possible to set real length associated with the size of a unit cell from the (real) size of minimal dunes on Earth deserts. As explained by Narteau et al (2009), a similar procedure could be used in order to simulate dunes in other environments using the model.

Besides Werner's algorithm, different phenomenological models with deterministic transport rules were proposed (Nishimori et al 1998, Katsuki et al 2005a). The model by Nishimori et al (1998), which included gravity-driven creep flux and lateral diffusion, also produced different dune geometries consistent with the imposed wind regime (see section 3.1.3). Katsuki et al (2005a) introduced a simpler model including empirical deterministic rules for saltation flux on the windward side of dunes and avalanches at the lee. This model was applied to study collisions between barchans (Katsuki et al 2011) and the genesis of a barchan belt from a sand source (Katsuki et al 2005b). Katsuki and Kikuchi (2011) found that a field of barchans forms only if the sand source is subject to a large influx, which is in qualitative agreement with the findings of Durán et al (2010), see section 3.2.1.4.

Insights gained from continuum and CA models for dunes can provide a basis for a different class of models where dunes are modeled as 'particles' interacting according to phenomenological rules (Diniega et al 2010). Such models can be potentially useful to assess the long time scale processes involved in the evolution of a field of barchans (Lima et al 2002, Durán et al 2011b) or transverse dunes (Parteli and Herrmann 2003, Lee et al 2005).

A conceptual model for the barchan–parabolic transition—Plants locally reduce the velocity of the wind blowing over dunes, thus inhibiting sand erosion and enhancing deposition (see section 4.1.1.2). The first places colonized by vegetation are those where erosion or deposition is small, i.e. the dune's arms, crest and the surrounding terrain. Indeed, erosion is an important factor for the dying of vegetation, since sand removal denudes the plants' roots and exposes deep layers to evaporation (Tsoar and Blumberg 2002). Since erosion prevents vegetation from growing on the windward side, the central part of a barchan or transverse dune can still migrate while two marginal ridges are left behind. A barchan that enters a stabilizing vegetated area thus undergoes a transition to a U-shaped parabolic dune (figure 32): its shape 'inverts', after which the vegetated limbs point opposite to the migration direction. When the rates of erosion or deposition at the crest become sufficiently small, vegetation finally stabilizes also the central part and the dune is rendered inactive (Tsoar and Blumberg 2002).

Zoom In Zoom Out Reset image size

Numerical simulations—The first successful modeling of the barchan–parabolic transition (Tsoar and Blumberg 2002) was achieved by Durán and Herrmann (2006b), who extended the model by Kroy et al (2002) to account for vegetation growth. The simulations indicated that the relevant parameter for the dune shape is the so-called fixation index (Durán and Herrmann 2006b),

$$\begin{equation} \theta\equiv Q_0/(V^{1/3}V_{\rm v}), \end{equation} \tag{ 3.8 }$$

where Q₀ is the saturated bulk sand flux, V is the dune volume and V_v is the characteristic growth rate of the vegetation cover. The fixation index gives the ratio between the erosion rate, which is proportional to Q₀/V^1/3, and the vegetation growth rate, V_v. Durán and Herrmann (2006b) found that parabolic dunes form when θ < θ_c (≈0.5), while barchans appear for larger θ (Durán and Herrmann 2006b). This prediction roughly matches field observations (Reitz et al 2010), and good quantitative agreement is found between the shape of real and simulated parabolic dunes (Durán et al 2008). Using the same model, Luna et al (2011) found that a sand beach evolves into a vegetated sand barrier (or 'foredune', see Hesp 1989), rather than into a field of barchans, when θ < 0.1.

In an earlier (and simpler) model, Nishimori and Tanaka (2001) extended the algorithm of Nishimori et al (1998) to include a phenomenological vegetation cover density with growth rate depending on local erosion and deposition rates. Geometries resembling barchans (parabolic dunes) were obtained for strong (weak) winds, i.e. for long (short) hop length of the slabs (Nishimori and Tanaka 2001), essentially the information encoded in equation (3.8). Also, Baas (2002) extended Werner's model (section 3.2.2) by adding a vegetation cover density, which linearly altered the probability of erosion at a given cell. Again, the vegetation growth rate depended on the local rates of erosion and deposition (Durán and Herrmann 2006b, Nishimori and Tanaka 2001). The barchan–parabolic transition was achieved by allowing different vegetation species (each associated with a maximum growth velocity and response rate to changes in erosion or deposition) to grow concurrently: pioneer species slowed down the flanks (first) for the stabilizer species to subsequently immobilize the edges into depositional ridges (Baas and Nield 2007).

Modeling has been used to investigate the evolution of large sand systems in response to variations in vegetation growth due to climatic changes (Nield and Baas 2008b) or altitude-induced changes in soil fertility (Pelletier et al 2009). In some places, barchans and parabolic dunes can even coexist under the same climatic conditions. This can be understood by noting that dunes with different sizes have different fixation indices (equation (3.8)). Yizhaq et al (2007) proposed a simple mathematical model that elucidates this bistability phenomenon. For moderate wind power condition, both barchans and parabolic dunes can form, while the former (latter) prevail for high (low) wind power (Yizhaq et al 2007).

Forrest and Haff (1992) introduced a discrete CA simulation following the trajectory of each moving particle. The initial ripples grew from the spatial asymmetry in the flux of impacting grains (screening instability) in agreement with Anderson (1987a). In the nonlinear regime of the bed evolution, merging of ripples automatically filtered out ripples of significantly different sizes (e.g. Anderson 1990). Forrest and Haff (1992) suggested that bedform merging should play a major role for attaining a fully developed state of ripples with constant size and spacing after a long time.

This hypothesis was explored in a toy-model introduced by Werner and Gillespie (1993), in which ripples were represented as one-dimensional entities ('worms') of different lengths moving around a discretized ring of N segments. The head of each 'ripple' moved one segment clockwise with probability inversely proportional to its length, thereby absorbing one segment of the 'downwind ripple'. Without any account for the detailed physics of ripples, the 'worm' model (as well as its two-dimensional version) produced a logarithmic increase of 'ripple' size with time (Werner and Gillespie 1993). Landry and Werner (1994) then introduced a three-dimensional model including a slope-dependent reptation length (Werner and Haff 1988) as well as rolling: after landing, the particle moved down the steepest descent (with a prescribed 'pseudomomentum') until halting on a stable pocket on the granular bed. Bedform spacing grew at a decreasing rate with time, however without indication of a steady-state value for the mean height and spacing of ripples (Landry and Werner 1994).

In a different approach, Pelletier (2009) applied an improved version of Werner's dune model (section 3.2.2) to ripple dynamics. Pelletier (2009) interpreted the (smaller) bedforms in the simulations as ripples, which then provided the roughness elements for simulations of dunes. Both ripples and dunes achieved a constant size and spacing, both increasing functions of grain size and excess shear stress (Andreotti et al 2006). As explained by the author, this steady state occurred because of the correction of the migration velocity of taller bedforms: taller dunes, subjected to larger values of bed shear stress (figure 27), migrated a little faster than in the original model by Werner (1995). In this manner, the model produced a smaller rate with which small and tall dunes merge, compared to the original Werner's model, eventually achieving a nearly stationary state where ripple spacing remained constant (Pelletier 2009). The author's findings suggested an intriguing genetic linkage between ripples and dunes via the scaling between bedform size and aerodynamic roughness. Such a model however differs from Anderson's picture, since it assumes that saltation (rather than reptation) is the transport mode for ripple migration (Werner 1995).

Prigozhin (1999), following Anderson (1990), improved the continuum model of Anderson (1987a) by allowing reptating grains to roll a certain distance upon landing, due to the action of gravity. The flux of rolling particles was described using the so-called BCRE theory of surface flow of granular materials over a static sand bed, which was introduced by Bouchaud et al (1994) for avalanche flows in sandpiles and adapted by Terzidis et al (1998) to the problem of ripple instability. An alternative model for rolling was also proposed by Hoyle and Woods (1997) and Hoyle and Mehta (1999). These authors and Prigozhin (1999) studied numerically the nonlinear regime of ripples growth (see also the theoretical works by Valance and Rioual (1999) and Misbah and Valance (2003)).

Towards a fully developed state—Following Prigozhin (1999), Yizhaq et al (2004) introduced a phenomenological correction that made the reptation flux slightly larger on the lee side. This correction artificially shifted the largest flux value closer to the ripple's crest, thus acting to limit the growth of the ripples. After a sufficiently long simulation time, the system apparently approached a saturated (fully developed) state with nearly constant ripple size (Andreotti et al 2006). In an improved version of Prigozhin's model, rolling occurred not only due to gravity but also due to aeolian drag (Manukyan and Prigozhin 2009). This model also accounted for mass exchange between the bed and the saltation cloud (Sauermann et al 2001). Erosion of ripple crests and deposition in throughs produced much more flat ripples than those obtained previously (Manukyan and Prigozhin 1999). These numerical simulations provided convincing evidence that the saturation of ripple growth (Andreotti et al 2006) occurs due to removal of sand from crests into throughs due to the action of the wind, as suggested by Bagnold (1941). Manukyan and Prigozhin (1999) also suggested that ripples produced in a wind tunnel longer than the saturation length should display both size and spacing independent of wind strength. This prediction remains to be confirmed experimentally (Andreotti et al 2006).

In sand beds with two grain sizes, ripples in general display coarse-grained crests, while in their inner parts and in the troughs smaller grains dominate (Sharp 1963, Hunter 1977). Ripples displaying this so-called 'inverse grading' have been named granule ripples or megaripples as they may achieve much larger wavelengths than ripples of uniform sand (Sharp 1963, Ellwood et al 1975, Tsoar 1990, Fryberger et al 1992, Sullivan et al 2005, Jerolmack et al 2006, Isenberg et al 2011).

Anderson and Bunas (1993) introduced the first model for aeolian ripples incorporating two grain sizes. Their model included a phenomenological description of the wind flow around the evolving topography, such as to produce a horizontally uniform wind speed above a 'ceiling height', roughly one ripple wavelength above the mean surface altitude. The trajectory of the reptating particles was different for coarse and fine grains: since smaller grains can jump higher, they are subjected to larger wind velocities, thus experiencing larger acceleration and hopping further. Small particles, once ejected from the ripples' crest, can jump trajectories long enough to land in the throughs. In contrast, coarse grains do not hop far enough to escape the lee, thus accumulating around the crest; after being buried during migration of the ripple, they reappear on the stoss surface, thereafter skittering up the crest to repeat the process. The difference in jump lengths for coarse and fine grains thus provides an explanation for sorting and 'inverse grading' of granule ripples (Werner et al 1986).

Following these findings, Makse (2000) extended previous continuum modeling (Terzidis et al 1998, Prigozhin 1999) to include a two-species model where coarse grains jumped shorter trajectories than smaller ones; and Manukyan and Prigozhin (2009) also adapted their model (section 3.2.4.2) to account for two grain sizes. The latter study showed that megaripples with larger grain-size ratios (of the order of 5) than those considered in previous models (Anderson and Bunas 1993, Makse 2000) can be obtained by accounting for wind-induced motion of the coarse particles armoring the ripple crests (Manukyan and Prigozhin 2009). Yizhaq (2004) also extended the model by Yizhaq et al (2004) to describe two types of sand flux, namely reptation of coarse grains and saltation of fines. This model links the wavelength of megaripples not to reptation, as in Anderson's model, but to saltation, as suggested by Bagnold (1941).

Other recent attempts to model granule ripples considered directly solving the Navier–Stokes equations for the wind flow over the topography (Wu et al 2008, Zheng et al 2008). One particular model (Zheng et al 2008), which considered three grain sizes, produced coarse-crested ripples as observed in nature and simulations. The model supported the conclusion that inverse grading of granule ripples results from the dependence of reptation length on grain size (Anderson and Bunas 1993).

Parteli et al (2005, 2007a) presented the first calculation of the shape of barchan dunes at Arkhangelsky crater on Mars (figure 33(a)). The simulations used the model by Kroy et al (2002) considering a grain diameter d = 500 µm (Edgett and Christensen 1991, see also section 2.4) and wind speeds larger than the corresponding threshold for direct entrainment under martian conditions, i.e. $u_{\ast_{\rm ft}}=2.0{\hbox{--}}.5\,{\rm m}\,{\rm s}^{-1}$ (Parteli et al 2007a). A value of impact threshold $u_{\ast_{\rm t}}=0.8u_{\ast_{\rm ft}}$ (≈2.0 m s⁻¹) was used in analogy with saltation on Earth (Bagnold 1941) (although more recent simulations and theoretical models have suggested a substantially lower ratio $u_{\ast_{\rm it}}/u_{\ast_{\rm ft}}$ for Mars; Almeida et al 2008, Kok 2010a, 2010b, Pähtz et al 2012; see section 2.4). From the minimal dune width (W_min), which scales with 1/[τ/τ_t − 1] (section 3.2.1.2), the authors obtained an estimate for the average shear velocity of sand-moving winds at Arkhangelsky crater, namely u_* = 3.0 m s⁻¹. Excellent quantitative agreement was found between the morphological relations of the barchans obtained in the simulation and those of the Arkhangelsky barchans (figure 33(a)).

Zoom In Zoom Out Reset image size

The calculations were then adapted to reproduce the shape of barchans at other locations on Mars (Parteli and Herrmann 2007b). In spite of the different minimal threshold velocity values predicted for each field owing to variations in atmospheric pressure, the values of u_* estimated from the simulations in all cases were around 3.0 m s⁻¹ (Parteli and Herrmann 2007b). This value is well within the maximum range of values of u_* on Mars (between 2.0 and 4.0 m s⁻¹) as it has been estimated indirectly from observations of surface erosion and ripple formation by the rovers (where estimations also assumed $u_{\ast_{\rm t}}\approx 2.0\,{\rm m}\,{\rm s}^{-1}$; Arvidson et al 1983, Moore 1985, Sullivan et al 2005). Observations from orbiters and landers suggested that such extreme values of wind shear velocity occur occasionally on Mars, during global dust storms (Arvidson et al 1983, Moore 1985, Sullivan et al 2005).

As mentioned above, there is strong theoretical evidence provided by recent works (Almeida et al 2006, Kok 2010a, 2010b, Pähtz et al 2012) for much lower values of impact threshold on Mars than the one (2.0 m s⁻¹) previously estimated (figure 21(b)). Future work should thus focus on reproducing the observed morphology of martian dunes using these revised estimates of the martian impact thresholds.

Experiments in a water tank and computer simulations have shown that some exotic dune shapes occurring on Mars could form under directionally varying flows (figures 33(b)–(d)). Taniguchi and Endo (2007) performed experiments of barchan ripples under reversing flow conditions, and found barchans with deformed limbs that resembled dunes at Proctor crater (figure 33(d)). Further, Parteli et al (2009) extended the simulations of dune formation under bimodal winds (see section 3.2.1.4) to low sand availability condition. The authors found that the resulting dune shape depends on θ_w (which determines dune alignment, see section 3.1.5.1) as well as on the duration of the bimodal wind relative to the dune migration time, i.e. t_w = T_w/T_m: a value of t_w smaller than 10⁻³ forms rounded barchans with reduced slip face ('occluded barchans'), if θ_w = 40°–80°; non-elongating 'wedge' dunes (figure 33(b)), if θ_w = 100° (also produced in water tank experiments by Hersen 2005), or elongating longitudinal bedforms (figure 33(c)), if θ_w > 110° (see Parteli et al 2009 for a phase diagram of the morphology.

Increased soil moisture is an important contributor to the immobilization of dunes. Crusting or induration of sand has been observed on Mars by Viking 1, Viking 2, Mars Pathfinder, Spirit and Opportunity landers (Arvidson et al 2004, Moore et al 1999, Thomas et al 2005, Sullivan et al 2008, 2011). The classical experiments by Kerr and Nigra (1952) illustrated the dramatic consequences sand induration can have for dune morphodynamics. The experiments consisted of spraying crude oil onto advancing barchans, in order to halt their movement. The sand arriving at the fixed dunes from the upwind was deposited at the lee, whereupon the dune was soaked with oil again. As this process continued, the slip face of the dunes became progressively smaller and each dune more elongated. The oil was sprayed on the dunes in successive stages until all motion was stopped, and the dunes achieved a rounded, elliptical dome-shaped form.

Rounded barchans occur in the north polar region of Mars known as Chasma Boreale (Schatz et al 2006). Their discovery raised the question whether sand induration due to frozen carbon dioxide or the presence of salts as an intergranular cement (Clarck et al 1982) could be a relevant agent for dune processes on Mars. Through using numerical modeling, Schatz et al (2006) succeeded in reproducing the transition from a crescent to a rounded dune shape (Kerr and Nigra 1952). Induration (non-mobility) was simulated for each successive stage of dune development, while at the same time new sand was added from upwind (figure 34). A further argument in favor of induration of the Chasma Boreale dunes is the surprising co-existence of the rounded barchans with straight-crested linear dunes aligned parallel to the orientation of the barchans—the linear dunes should simply decay in barchans if the wind was unidirectional (section 3.2.1.4). Possibly, sand saltating downwind along both sides of the indurated linear dune could deposit at the lee due to secondary flow effects (not included in the model by Schatz et al 2006) at the dune's downwind extremity (Tsoar 2001), thus yielding dune elongation.

Zoom In Zoom Out Reset image size

The presence of soil moisture can create substantial interparticle forces that inhibit the initiation of saltation, especially for sandy soils (Chepil 1956, Belly 1964, McKenna Neuman and Nickling 1989, Sherman et al 1998). For low relative humidities (below ∼65%), these interparticle forces are produced primatirly by bonding of adjacent adsorbed water layers (hygroscopic forces), whereas for high relative humidities (above ∼65%) this occurs primarily through the formation of water wedges around points of contact (capillary forces) (Hillel 1980, Ravi et al 2006, Nickling and McKenna Neuman 2009).

Water adsorption is governed by electrostatic interactions of the mineral surface with the polar water molecules. Since sandy soils generally contain a lower density of net electric charges, substantially less water can be adsorbed onto sandy soils than onto clayey soils (Hillel 1980). Consequently, water bridges form in sandy soils at a relatively low soil moisture content, thereby producing substantial capillary forces (e.g. Belly 1964). McKenna Neuman and Nickling (1989) used this observation to derive an expression for the increase of the saltation fluid threshold for sand due to the presence of soil moisture. Fécan et al (1999) then generalized the McKenna Neuman and Nickling expression to all soil types, obtaining the following empirical expressions:

$$\begin{eqnarray} &&\frac{u_{\ast_{\rm wt}}}{u_{\ast_{\rm ft}}}=1\tqs (w < w'),\\ &&\frac{u_{\ast_{\rm wt}}}{u_{\ast_{\rm ft}}}=\sqrt {1+1.21( {w-w'} )^{0.68}}\tqs (w \ge w'), \end{eqnarray} \tag{ 4.1 }$$

where $u_{\ast_{\rm wt}}$ and $u_{\ast_{\rm ft}}$ are respectively the fluid saltation thresholds in the presence and absence of soil moisture, and w is the soil's volumetric water content in percent. Since water adsorbs more readily onto clay soils (Hillel 1980), the volumetric fraction of water that the soil can absorb before capillary forces become substantial, w', increases with the soil's clay content. Using curve fitting of equation (4.1) to experimental results, Fécan et al (1999) obtained

$$\begin{equation} w'=0.17c_{\rm s} +0.0014c_{\rm s}^{\rm 2} , \end{equation} \tag{ 4.2 }$$

where c_s is the soil's clay content in percent.

Although the empirical parametrization of equations (4.1) and (4.2) ignores interparticle forces due to bonding of adsorbed water layers, which is likely substantial for soils with either a high clay content or a low soil moisture content (Hillel 1980, Ravi et al 2004, 2006), this fairly straightforward parametrization seems to reasonably match measurements and is widely used to parametrize the effect of soil moisture on the dust emission threshold in atmospheric circulation models (e.g. Zender et al 2003a).

In addition to soil moisture, the saltation fluid threshold is also affected by the presence of non-erodible roughness elements such as pebbles, rocks and vegetation. Since the flow resistance of these objects is substantially larger than that of bare soil, they extract momentum from the wind. The presence of roughness elements therefore reduces the wind shear stress on the intervening bare soil and increases the total threshold wind stress required to initiate saltation and dust emission (Raupach et al 1993).

The effect of roughness elements on the saltation fluid threshold depends on the partitioning of the total wind stress between the fraction absorbed by the bare soil (τ_S) and the fraction absorbed by the roughness elements (τ_R) (Schlichting 1936),

$$\begin{equation} \tau =\tau _{\rm S} +\tau _{\rm R}. \end{equation} \tag{ 4.3 }$$

The difficulty now lies in relating the partitioning of the fluid drag to the properties of the roughness elements. The absorption of fluid drag by roughness elements is largely determined by the frontal area presented to the flow (Marshall 1971), which is captured in the roughness density λ (e.g. Raupach 1992),

$$\begin{equation} \lambda =nbh, \end{equation} \tag{ 4.4 }$$

where n is the number of roughness elements per unit area, and b and h are the average width and height of the roughness elements. By first relating the fluid drag absorbed by roughness elements to the roughness density λ and then assuming that the area and volume of the wakes of different roughness elements are randomly superimposed, Raupach (1992) derived an approximate expression for the fraction of the fluid drag absorbed by the bare soil,

$$\begin{equation} \frac{\tau _{\rm S}}{\tau}=\frac{1}{1+\beta \lambda}, \end{equation} \tag{ 4.5 }$$

where β is the ratio of the drag coefficients of a typical roughness element and the bare soil, which is on the order of ∼100 for typical conditions. Raupach et al (1993) then expanded the Raupach (1992) drag partitioning theory to estimate the effect of roughness elements on the saltation fluid threshold and obtained an expression for the ratio R of the fluid threshold without and with roughness elements (Gillette and Stockton 1989):

$$\begin{equation} R=\frac{u_{{\ast_{\rm ft,bare}}}}{u_{{\ast_{\rm ft,rough}}}}=\left[ {\frac{1}{\left( {1-m\sigma \lambda} \right)\left( {1+m\beta \lambda} \right)}} \right]^{1/2}, \end{equation} \tag{ 4.6 }$$

where σ is the ratio of the basal area to the frontal area for the roughness elements, and m (0 < m ⩽ 1) is a factor meant to account for the fact that the fluid saltation threshold of a homogeneous soil is not determined by the average shear stress required for particle mobilization, but rather by the maximum shear acting on any erodible point on the surface. Raupach et al (1993) recommended m ≈ 0.5 for flat, erodible surfaces based on comparisons with the experimental data of Gillette and Stockton (1989). Subsequent wind tunnel and field measurements of the drag partition and the ratio R have generally confirmed the validity of equations (4.5) and (4.6), although most experiments had to adapt the values of σ, m and β from the values proposed by Raupach (Wyatt and Nickling 1997, Crawley and Nickling 2003, King et al 2005, Gillies et al 2007).

Despite the reasonable agreement of equations (4.5) and (4.6) with measurements, the roughness density λ does not fully capture the dependence of the drag partition on the physical properties, dimensions, spatial distribution, and orientation of the roughness elements (Gillies et al 2006, Nickling and McKenna Neuman 2009). For instance, the drag partition is affected by the porosity (Shao 2008) and Reynolds number dependency of the roughness element drag coefficients (Gillies et al 2002), which is highly relevant for determining the effects of vegetation on the saltation fluid threshold (Okin 2005, 2008). Moreover, the spatial distribution of the roughness elements can play an important role in determining their effect on the fluid threshold. For example, Okin and Gillette (2001) found that the distribution of vegetation in a semi-arid mesquite landscape adapts in such a manner as to form 'streets' of readily erodible bare soil, for which the saltation fluid threshold is substantially lower than would be predicted with equation (4.6). The non-uniformity of the spatial distribution of vegetation thus likely results in higher saltation and dust emission fluxes than would be predicted using the Raupach model (Okin 2005, 2008), which assumes a random distribution (Raupach 1992).

Although the proportionality constants C_S, C_F, C_K and C_GP in equations (4.12)–(4.16) are independent of u_*, these parameters depend on how efficiently impacting kinetic energy is converted into emitted dust aerosols. Specifically, soil properties such as the clay content (Marticorena and Bergametti 1995), the bonding strength of soil dust aggregates (Shao et al 1993a, Rice et al 1996), the dry aggregate size distribution (Chepil 1950, Zobeck 1991b), the presence of soil crusts (Belnap and Gillette 1998, Gomes et al 2003), and the soil plastic pressure (Lu and Shao 1999) are thought to determine these proportionality constants. Unfortunately, the exact dependence of the vertical dust flux on these physical parameters is both poorly understood and requires detailed knowledge of soil properties that are not normally available on regional or global scales (e.g. Shao 2001, Laurent et al 2008). Extensive further study is thus required to better relate the vertical dust flux to soil properties.

In order to parameterize observed regions of high dust emission potential despite these problems, atmospheric circulation models generally use empirical parametrizations (e.g. Marticorena and Bergametti 1995, Ginoux et al 2001, Tegen et al 2002, Zender et al 2003a). In particular, the soil's ability to produce dust aerosols is often captured in an empirical soil erodibility function. Prospero et al (2002) noted from multi-decadal observations by the Total Ozone Mapping Spectrometer (TOMS) onboard the Nimbus 7 satellite that areas with high dust loading coincided with topographic lows, which they hypothesized to be due to the fluvial deposition of fine-grained material in topographic lows. This finding was subsequently used to propose empirical soil erodibility parametrizations (Ginoux et al 2001, Tegen et al 2002, Zender et al 2003b); for instance, Ginoux et al (2001) proposed an empirical formulation of the soil erodibility based on the relative height of a model grid box compared to its surroundings. However, advection of dust plumes after emission causes the high dust loading observed by TOMS to be shifted downwind from dust source regions (Schepanski et al 2007), which later studies indicate are actually concentrated in mountain foothill regions and dry lake beds (Schepanski et al 2009). Although the use of soil erodibility parametrizations improves model agreement with measurements (Zender et al 2003b, Cakmur et al 2006), the empirical nature of these parametrizations suggests that improvements in understanding the dependence of the vertical dust flux on soil properties could bring corresponding improvements in simulations of the dust cycle. Indeed, current model simulations of the global dust cycle still show substantial discrepancies with measurements (e.g. Cakmur et al 2006).

The dust production model (DPM) assumes that sandblasting results in the emission of three separate lognormal modes of dust aerosols, with the relative contribution of each mode determined by its bonding energy and the kinetic energy of impacting saltators (Alfaro et al 1997, Alfaro and Gomes 2001). Using wind tunnel experiments with two separate soils, Alfaro et al (1998) determined the median diameters of the three modes to be 1.5, 6.7 and 14.2 µm. However, recent field measurements suggest that the bonding energies and median diameter of these three lognormal modes might need to be adapted for specific soils (Sow et al 2011).

Since the rupturing of interparticle bonds requires energy, the DPM predicts that larger saltating particle impact energies produce more disaggregated and thus smaller dust aerosols, as also suggested by experiments showing that glass broken at larger input energies produces a larger fraction of small fragments (Scheibel et al 1990, Weichert 1991). Because the DPM further assumes that the saltator impact speed is proportional to wind speed, it predicts a shift to smaller aerosol sizes with increasing wind speed. Although the assumption that saltator impact speeds scale with u_* conflicts with measurements for transport limited saltation (figure 15), several wind tunnel studies indeed observed a shift to smaller aerosol sizes with increasing wind speed (Alfaro et al 1997, 1998, Alfaro 2008). Consequently, the theory of Alfaro and Gomes (2001) is in agreement with these wind tunnel measurements (figure 38(a)).

Zoom In Zoom Out Reset image size

The possible dependence of the dust size distribution on the wind speed is discussed in more detail in section 4.1.3.4.

A second theory for the size distribution of emitted dust aerosols was formulated by Shao (2001, 2004) and is based on the insight that the emitted dust size distribution must be bound by the two extreme states of the soil: the undisturbed, minimally disaggregated state, and the disturbed, fully disaggregated state. That is,

$$\begin{equation} p_{\rm d} ( {D_{\rm p}} )=\gamma p_{\rm m} ( {D_{\rm p}} )+( {1-\gamma} )p_{\rm f} ( {D_{\rm p}} ), \end{equation} \tag{ 4.17 }$$

where p_d, p_m and p_f are respectively the particle size distributions of the emitted dust aerosols, the minimally disturbed soil and the fully disaggregated soil. Using the assumption that the saltator impact speed scales with u_*, Shao postulated a weighting factor γ that increases monotonically with the energy of impacting saltators:

$$\begin{equation} \gamma =\exp [ {-k( {u_\ast -u_{\ast {\rm t}}} )^n} ], \end{equation} \tag{ 4.18 }$$

where k and n are empirical coefficients that are determined from fitting of equations (4.17) and (4.18) to experimental data. An important similarity of the Shao model to the DPM is that it also predicts that an increase in u_* results in more disaggregated, and hence smaller, dust aerosols. Like the DPM, the Shao (2001, 2004) model is in agreement with the Alfaro wind tunnel measurements (figure 38(a)), although it requires the n and k parameters to be tuned accordingly. However, Shao et al (2011b) recently found that a constant value of γ (i.e. independent of u_*) produced the best agreement with field measurements.

The theories of Alfaro and Gomes (2001) and Shao (2001, 2004) provide a framework for predicting the dust size distribution of emitted dust aerosols. However, the use of these theories for regional and global circulation model parametrizations of dust emission is limited by the poor understanding of the interparticle forces that determine both the energy required to rupture the cohesive bonds between soil particles (in the case of the DPM; see Sow et al 2011) and the soil's aggregation state (in the case of Shao's theory). Moreover, the DPM and especially Shao's theory require knowledge of the soil size distribution, which is not normally available on regional and global scales.

As a way to circumvent the problem of the poorly understood interparticle forces and bonding energies and the resulting use of empirical coefficients, Kok (2011a) developed an alternative theory for the size distribution of emitted dust aerosols, assuming that most dust emission is due to brittle material fragmentation. This idea was previously noted by Gill et al (2006) based on the observation from soil science that stressed dry soil aggregates fail as brittle materials (Lee and Ingles 1968, Braunack et al 1979, Perfect and Kay 1995, Zobeck et al 1999). Since the size distribution of fragments resulting from brittle material fragmentation is determined by the patterns in which cracks nucleate, propagate, and merge in the material (Astrom 2006), the analogy of dust emission with brittle material fragmentation mostly eliminates the need for a detailed understanding of interparticle forces to predict the emitted dust size distribution. Moreover, since the pattern with which cracks propagate and merge is largely scale-invariant (Astrom 2006), the size distribution of fragments produced by brittle material fragmentation follows a simple power law, with the exponent of the power law only weakly dependent on the brittle object's shape (Oddershede et al 1993) and seemingly invariant to the material type (see figure 1 in Kok 2011a). (Note that scale invariance occurs for many natural phenomena, including avalanches, the fragmentation of rocks and atomic nuclei, and earthquakes (Gutenberg and Richter 1954, Turcotte 1986, Bak et al 1987, Bondorf et al 1995, Astrom et al 2004).)

Evidence that the fragmentation of dust aggregates is indeed a form of brittle material fragmentation is that the power law that describes the size distribution of fragments produced by brittle fragmentation also describes the size distribution of emitted dust aerosols in a limited size range (∼2–10 µm; see figure 2 in Kok 2011a). The emission of aerosols >10 µm is depleted relative to this power law because of the finite propagation distance λ of the side branches of cracks created in the dust aggregate by a fragmenting impact (Astrom 2006). Conversely, the emission of aerosols <2 µm is depleted relative to the power law because the discrete particles composing the dust aggregates have a typical size of ∼1 µm, which inhibits the creation of smaller fragments. Since many dust aerosols are aggregates themselves (Okada et al 2001, E A Reid et al 2003, Chou et al 2008), Kok assumed that the production of dust aerosols with size D_d is proportional to the volume fraction of soil particles with size D_s ⩽ D_d that can contribute to the formation of the aerosol. By combining this assumption with the size distribution resulting from brittle material fragmentation (reviewed in Astrom 2006), Kok (2011a) derived relatively straightforward analytical expressions for the number and volume size distribution of emitted dust aerosols:

$$\begin{eqnarray} &&\fl\frac{\rmd N_{\rm d}}{\rmd\ln D_{\rm d}}=\frac{1}{c_{\rm N} D_{\rm d}^2}\left[ \!{1+{\rm erf}\left( {\frac{\ln ( {D_{\rm d} /\overline {D_{\rm s}}} )}{\sqrt 2 \ln \sigma _{\rm s}}} \right)} \right]\exp \left[\! {-\left( {\frac{D_{\rm d}}{\lambda}} \right)^3} \!\right],\nonumber\\ \end{eqnarray} \tag{ 4.19 }$$

and

$$\begin{eqnarray} &&\fl\frac{\rmd V_{\rm d}}{\rmd\ln D_{\rm d}}=\frac{D_{\rm d}}{c_{\rm V}}\left[ {1+{\rm erf}\left( {\frac{\ln \left( {D_{\rm d} /\overline {D_{\rm s}}} \right)}{\sqrt 2 \ln \sigma _{\rm s}}} \right)} \right]\exp \left[ {-\left( {\frac{D_{\rm d}}{\lambda}} \right)^3} \right],\nonumber\\ \end{eqnarray} \tag{ 4.20 }$$

where N_d and V_d are respectively the normalized number and volume of dust aerosols with size D_d, c_N = 0.9539 µm⁻² and c_V = 12.62 µm are normalization constants, $\overline{D_{\rm s}} \approx 3.4\,\mu \,{\rm m}$ and σ_s ≈ 3.0 are the median diameter by volume and the geometric standard deviation of the log-normal distribution that describes a 'typical' arid soil size distribution in the PM20 size range, erf is the error function, and λ ≈ 12 µm denotes the propagation distance of the side branches of cracks created in the dust aggregate by a fragmenting impact.

The brittle fragmentation theory applies only to dust emission events that are predominantly due to the fragmentation of soil aggregates. Therefore, equations (4.19) and (4.20) are for instance not valid for (i) aerodynamically lifted dust, (ii) dust emitted mainly by impacts in the damage regime (Kun and Herrmann 1999), which could occur either for very cohesive soils or for sandy soils where most of the PM20 dust exists as coatings on larger sand grains (Bullard et al 2004), and (iii) dust with diameters larger than ∼20 µm, which are more likely to occur as loose particles in the soil (Alfaro et al 1997, Shao 2001), such that their emission is not predominantly due to fragmenting impacts. Despite these limitations, the good agreement of equations (4.19) and (4.20) with seven different data sets of field measurements (figure 38(b)) suggests that many natural dust emission events are due to aggregate fragmentation and can thus be reasonably described by brittle fragmentation theory. Further research is required to investigate this issue.

Measurements and theories have yielded conflicting results on whether the dust size distribution depends on the wind speed at emission. Whereas the Alfaro and Gomes (2001) and Shao (2001, 2004) theories predict a decrease in the mean dust aerosol size with u_*, the brittle fragmentation theory of dust emission (equations (4.19) and (4.20)) predicts that the size distribution is independent of u_*. Similarly, the wind tunnel measurements of Alfaro and colleagues (Alfaro et al 1997, 1998, Alfaro 2008) reported a strong dependence of the emitted dust size distribution on u_*, whereas some field measurements have not found such a dependence (Gillette et al 1974, Shao et al 2011b). The study of Kok (2011b) attempted to settle this issue using a statistical analysis of the mean aerosol diameter of published field measurements (Gillette 1974, Gillette et al 1974, Sow et al 2009, Shao et al 2011b). This study showed that individual data sets had opposing and statistically insignificant trends with u_*. Moreover, a compilation of all published field measurements showed no statistically significant trend with u_*, indicating that the size distribution of emitted dust aerosols does not depend on the wind speed at emission. A similar conclusion was reached by Reid et al (2008), based on the similarity of measured size distributions of dust advected with different wind speeds at emission.

Kok (2011b) also offered a possible explanation for the discrepancy between field measurements on the one hand and wind tunnel measurements and the theories of Alfaro and Gomes (2001) and Shao (2001, 2004) derived from them on the other hand. The measurements of Alfaro et al (1997, 1998) were probably not in steady state since the wind tunnel was only 3.1 meters in length (Alfaro et al 1997, p 11243). Consequently, the saltator impact speed could have increased with u_*, which in turn could have produced a shift to smaller aerosols as proposed by Alfaro and colleagues. Similarly, the dependence of the emitted dust size distribution on u_* in the DPM and Shao theories is due to the assumption in these models that the saltator impact speed scales with u_*, which measurements indicate is incorrect for transport limited steady-state saltation (figure 15).

All three theories reviewed above also include a dependence on the soil size distribution, and many regional and global circulation models also account for a soil dependence of the emitted dust size distribution (e.g. Ginoux et al 2001, Laurent et al 2008). However, the small amount of scatter between the dust flux measurements compiled in figure 38(b), despite the widely varying soil types for which these measurements were made, suggests that changes in soil conditions have only a limited effect on the emitted dust size distribution. A similar conclusion is suggested by the insensitivity of dust aerosol size distributions to changes in the source region (Reid et al 2008). More research is clearly needed to better quantify the influence of the soil size distribution on the emitted dust size distribution, but the apparently limited dependence of the emitted dust size distribution to the soil size distribution is highly fortuitous for regional and global dust modeling of the dust cycle.

LLJs are horizontal winds characterized by a maximum of about 15 m s⁻¹ in the frictionally decoupled layer immediately above the surface layer (e.g. Blackadar 1957, Holton 1967). They are most commonly observed at nighttime and can occur over all continents above both flat and complex terrain and may extend over tens to hundreds of kilometers (May 1995, Davis 2000, Schepanski et al 2009). During nights with low surface wind speeds, near-surface air layers are well stratified and turbulence is suppressed. In such conditions, air layers above the near-surface are frictionally decoupled from the surface and hence are associated with high wind speeds (e.g. Hoxit 1975, Garratt 1992, Mauritsen and Svensson 2007, Mahrt 1999).

After sunrise, convective turbulence arises with the onset of solar heating at the surface. The decoupled air layer aloft becomes frictionally coupled to the surface, and LLJ momentum is mixed down (Blackadar 1957, Lenschow et al 1979). As a consequence, high surface wind speeds occur until the LLJ is eroded by the convection in the planetary boundary layer (PBL) set on by the solar heating at the surface. The LLJ peaks at nighttime while surface wind peaks in mid-morning when LLJ momentum is mixed down to the surface (figure 41, Washington et al 2006, Todd et al 2008). Over the Sahara, LLJs occur under clear skies and low surface wind speed conditions (e.g. Thorpe and Guymer 1977) and have been observed within the northeasterly 'Harmattan' flow (Washington and Todd 2005). About 65% of dust activity over the Sahara is due to the break-down of the nocturnal LLJs (e.g. Schepanski et al 2009). In particular, dust emission throughout the year in the Bodélé region in Chad is mainly related to LLJs dynamics (Washington et al 2006).

Zoom In Zoom Out Reset image size

Deep moist convection and the associated downdrafts of cold, humid air is another meteorological phenomenon capable of generating the high surface wind speeds required for dust emission. Such outflows take the form of density currents and can propagate many hundreds of kilometers away from the parent convective system (Simpson 1997, Williams et al 2009). Over arid and semiarid areas, these density currents cause strong dust emission (of the order of 1 Tg per event (Bou Karam et al 2011)). Over North and West Africa, such events are called 'haboobs' (figure 42) and are observed during the summer monsoon season. Haboobs are frequent during afternoons and evenings, when moist convection is at maximum (Sutton 1925, Idso et al 1972, Droegemeier and Wilhelmso 1987, Flamant et al 2007, Bou Karam et al 2008, 2011).

Zoom In Zoom Out Reset image size

Deep moist convection can also occur in mountain areas due to blocking of the atmospheric flow by the orography. Large scale density currents can form resulting from the subsequent evaporative cooling of cloud particles or precipitation. The high surface wind speeds associated with these density currents can lead to dust emission over arid and semiarid areas (Droegemeier and Wilhelmson 1987). Over the Sahara, such dust events can be observed for example in the Atlas Mountains (Knippertz et al 2007) and in the Hoggar Mountains domain (Cuesta et al 2010).

Cyclones are considered as major dynamic features for the mobilization and the upward mixing of dust by virtue of the associated strong near-surface cyclonic winds and convective turbulence (e.g. Liu et al 2003, Bou Karam et al 2010). This mechanism operates in two steps; first the dust is mobilized by strong momentum from the cyclonic surface wind and then it is mixed upward to high altitudes by the strong turbulence and the systemic ascending motion associated with the cyclone dynamics (Liu et al 2003, Bou Karam et al 2010). For example, the Mongolian Cyclone (figure 43(a)), affects eastern Asia in March, April and May. This intense system is associated with the East Asian trough (a trough is an elongated region of low atmospheric pressure, often associated with fronts) and leads to strong northwesterly near-surface winds up to 18 m s⁻¹, thereby generating dust storms (Shao and Wang 2003). In general, Asian dust storms are mainly associated with cyclonic cold fronts during the surges of cold continental air masses in late winter and early spring when most of the area is under the influence of the powerful Siberian–Mongolian anticyclone (e.g. Watts 1969, Littmann 1991).

Zoom In Zoom Out Reset image size

Cyclones are also observed over the Sahara desert in spring months when temperature gradients between the North African coast and the Mediterranean Sea are strongest (Pedgley 1972, Alpert and Ziv 1989, Trigo et al 2002). Observations indicate a frequent initiation of Saharan cyclones on the leeward side east and south of the Atlas Mountains (Barkan et al 2005, Prezerakos et al 1990, Alpert et al 1990, Alpert and Ziv 1989). Hence, both the lee-effect of the mountains and the coastal thermal gradient effect can explain the spring cyclogenesis, initiated by the presence of an upper-level trough to the west (e.g. Horvath et al 2006, Egger et al 1995, Schepanski et al 2009, Bou Karam et al 2010). Saharan cyclones are characterized by an active warm front, heavy dust storms and low visibilities, and a sharp cold front that is well defined at the surface by changes in temperature of 10–20 K (figure 43(b)). They move quickly eastward (faster than 10 m s⁻¹) mostly following the North African coast (Alpert et al 1990, Alpert and Ziv 1989, Bou Karam et al 2010). The dust load associated with Saharan cyclones is estimated to be of the order of 8 Tg per event (Bou Karam et al 2010).

Large and persistent outbreaks of dust over arid and semiarid regions can also be caused by cold fronts of synoptic scale (e.g. figure 44). These cold fronts mainly originate from the penetration of upper-level troughs from high latitudes into low latitudes (e.g. Jankowiak and Tanré 1992, Knippertz and Fink 2006, Tulet et al 2008, Liu et al 2003). The dust fronts generated during these events are related to density currents caused by strong evaporative cooling along the precipitating cloud-band that accompany the penetration of the upper-level cold front into arid regions (e.g. Knippertz and Fink 2006). During such storms, the dust is pushed ahead of the cold front into the rising warm air over the desert. Over northern China, most of the dust storms during spring are caused by cold fronts that form in connection with the dynamic of the Siberian High.

Zoom In Zoom Out Reset image size

Furthermore, cold-front dust emission can occur together with the surge of monsoon flows from the ocean toward the continent. Bou Karam et al (2008) have identified this new mechanism for dust emission over West Africa during the summer monsoon season. Highly turbulent winds at the leading edge of the monsoon nocturnal flow in the inter tropical front (ITF), the ITF marks the interface in the lower atmosphere between the moist southwesterly monsoon flow and the hot and dry northeasterly harmattan flow over Africa) region have been identified to generate dust uplifting. Also, Marsham et al (2008) have described, via in situ measurements, the characteristics of the dusty layer at the leading edge of the monsoon flow over West Africa. Bou Karam et al (2009) have estimated dust emission associated with the nocturnal monsoon flow over Niger to be of the order of 0.7 Tg per day. This mechanism is most active during the night when the monsoonal surge usually occurs.