An improved dust emission model – Part 1 : Model description and comparison against measurements

Simulations of the dust cycle and its interactions with the changing Earth system are hindered by the empirical nature of dust emission parameterizations in weather and climate models. Here we take a step towards improving dust cycle simulations by using a combination of theory and numerical simulations to derive a physically based dust emission parameterization. Our parameterization is straightforward to implement into large-scale models, as it depends only on the wind friction velocity and the soil’s threshold friction velocity. Moreover, it accounts for two processes missing from most existing parameterizations: a soil’s increased ability to produce dust under saltation bombardment as it becomes more erodible, and the increased scaling of the dust flux with wind speed as a soil becomes less erodible. Our treatment of both these processes is supported by a compilation of qualitycontrolled vertical dust flux measurements. Furthermore, our scheme reproduces this measurement compilation with substantially less error than the existing dust flux parameterizations we were able to compare against. A critical insight from both our theory and the measurement compilation is that dust fluxes are substantially more sensitive to the soil’s threshold friction velocity than most current schemes account for.


Introduction
The emission of mineral dust aerosols produces important impacts on the Earth system, for instance through interactions with radiation, clouds, the biosphere, and atmospheric chemistry (e.g., Miller and Tegen, 1998;Jickells et al., 2005;Cwiertny et al., 2008;Creamean et al., 2013).The inclusion of an accurate dust cycle in climate and weather models is thus critical.Yet, the current generation of dust modules shows substantial disagreements with measurements (Cak-Published by Copernicus Publications on behalf of the European Geosciences Union.mur et al., 2006;Huneeus et al., 2011;Evan et al., 2014), and commonly uses semiempirical "dust source functions" to help parameterize dust emission processes (e.g., Ginoux et al., 2001;Tegen et al., 2002;Zender et al., 2003b).
Here we aim to improve the dust cycle's representation in weather and climate models, in particular for climate regimes other than the current climate to which most models are tuned (Cakmur et al., 2006).We do so by presenting a physically based theory for the vertical dust flux emitted by an eroding soil.The functional form of the resulting dust flux parameterization is supported by a compilation of quality-controlled dust flux measurements, and our new parameterization reproduces these measurements with substantially less error than the existing parameterizations we are able to test against.Moreover, our new parameterization is relatively straightforward to implement since it uses only variables that are readily available in large-scale models.A critical insight from the theory is that the dust flux is substantially more sensitive to changes in the soil state than most climate models account for.
We derive our new dust emission parameterization in Sect.2, after which we compare our parameterization's predictions against a compilation of quality-controlled vertical dust flux measurements in Sect.3. We discuss the implications of the new parameterization and conclude the article in Sect. 4.

Derivation of physically based dust flux parameterization
Because of their small size, dust particles in soils ( <62.5 µm diameter; Shao, 2008) experience cohesive forces that are large compared to aerodynamic and gravitational forces.Consequently, dust aerosols are usually not lifted directly by wind (Gillette et al., 1974;Shao et al., 1993;Sow et al., 2009) and instead are emitted through saltation, in which larger sand-sized particles (∼70-500 µm) move in ballistic trajectories (Bagnold, 1941;Shao, 2008;Kok et al., 2012).Upon impact, these saltating particles can eject dust particles from the soil, a process known as sandblasting.Moreover, some saltating particles are actually aggregates containing dust particles.Upon impact, these aggregates can also emit dust aerosols (Shao et al., 1996).We aim to obtain an analytical expression that captures the main dependencies of the emitted flux of dust aerosols on wind speed and soil properties.An important limitation is that, to allow its implementation into climate models, this expression can only use parameters that are globally available.Our approach to achieve this objective combines a theoretical derivation with numerical simulations of dust emission.We start in the next section by providing a basic theoretical expression for the vertical dust flux, after which we derive the three main variables in this expression in the three subsequent sections.We then combine all these components together to give the full dust emission parameterization in Sect.2.5.

Basic theoretical expression of the vertical dust flux
The starting point of our theory is the insight that a saltator impact will produce dust emission only if a threshold impact energy is exceeded (Rice et al., 1999), with the nature and value of this threshold depending on the soil type and state.For instance, for a soil with only a small fraction of suspendable particles, much of the dust is present as coatings on larger sand particles (Bullard et al., 2004), such that the relevant threshold is likely the energy required to rupture these coatings (Crouvi et al., 2012).Conversely, for a soil containing a large fraction of suspendable dust particles, the threshold for fragmentation of brittle dust aggregates could be most important (Kok, 2011b).Since the theoretical size distribution predicted by brittle fragmentation theory is in good agreement with dust size distribution measurements (Albani et al., 2014;Mahowald et al., 2014;Rosenberg et al., 2014), and its implementation into large-scale models improves agreement with other measurements of the dust cycle (Johnson et al., 2012;Nabat et al., 2012;Li et al., 2013;Evan et al., 2014), the threshold for fragmentation of soil dust aggregates might be the most relevant threshold for dust emission under many conditions.For simplicity, we thus assume that the energy required for dust aggregate fragmentation is globally the most relevant dust emission threshold, but we note that the functional form of the dust flux parameterization derived below is likely relatively insensitive to the chosen threshold process (see further discussion in Sect. 3.6).Following the discussion above, the vertical dust flux F d (kg m −2 s −1 ) generated by a soil during saltation can be written as where f bare is the fraction of the surface that consists of bare soil; n s is the number of saltator impacts on the soil surface per unit area and time; f frag is the average fraction of saltator impacts resulting in fragmentation of either the impacted soil dust aggregate, or the saltator if that is an aggregate itself; m frag is the mean mass of emitted dust produced per fragmenting impact; and ε is the mass fraction of emitted dust that does not reattach to the surface and is transported out of the near-surface layer where it can be measured (Gordon and McKenna Neuman, 2009).Since ε likely depends predominantly on the flow immediately above the surface, which remains relatively constant with wind speed (Ungar and Haff, 1987;Shao, 2008;Kok et al., 2012), we expect ε to be approximately constant for different wind conditions for a given soil.Finally, we obtain n s from the balance of horizontal momentum in the saltation layer (Shao et al., 1996;Kok et al., 2012): where τ s denotes the wind stress exerted on the bare soil, and τ st denotes the threshold value of τ s above which saltation occurs.Furthermore, m s and v imp are the mean saltator mass and impact speed, and the constant C ns ≈ 2 (Kok et al., 2012).Substituting Eq. (2) into Eq.( 1) yields where we assumed that m frag / m s = γ f clay .That is, we assumed that m frag / m s scales with the volume fraction of the soil that contributes to the creation of dust aerosols (Sweeney and Mason, 2013).The size limit of dust relevant for climate is usually taken as ∼10 µm (Mahowald et al., 2006(Mahowald et al., , 2010)), but since the mass fraction of soil particles ≤10 µm is not available on a global scale, we instead use the soil clay fraction (f clay ; ≤2 µm diameter), which is globally available (FAO, 2012).The dimensionless coefficient γ likely depends on the relative sizes of soil dust aggregates and saltators.Because many saltators are aggregates (Shao, 2008), we expect only modest variations in γ between soils and take it as a constant.
Since we thus expect variations of γ and ε with wind and soil conditions to be less important (see above), we seek to understand the dependence of τ s , v imp , and f frag on wind and soil conditions in order to complete our theoretical expression for F d .In the next three sections, we derive these dependencies through a combination of insights from previous studies, new theoretical work, and simulations with the numerical saltation model COMSALT (Kok and Renno, 2009).

Friction velocity and the wind stress τ s on the bare soil surface
The dust flux emitted by an eroding soil depends on both the soil's properties and on the wind shear stress τ exerted on the surface (Marticorena and Bergametti, 1995;Shao et al., 1996;Alfaro and Gomes, 2001;Shao, 2001;Klose and Shao, 2012;Kok et al., 2012).This shear stress is characterized by the friction velocity, which is defined as (e.g., Bagnold, 1941;Shao, 2008;Kok et al., 2012) where ρ a is the air density.Dust emission often occurs in the presence of nonerodible elements such as rocks and vegetation.Thus, τ can be partitioned between the stress τ R exerted on nonerodible roughness elements and the stress τ s exerted on the bare erodible soil; only τ s produces dust emission (Raupach et al., 1993;Shao et al., 1996).In analogy with Eq. ( 4), we define the soil friction velocity corresponding to τ s as where f bare is the fraction of the surface that consists of bare, erodible soil (note that f bare corresponds to the quantity S / S in the terminology of Raupach, 1992).The soil friction velocity u * can be derived from u * using knowledge of the soil's roughness elementsf bare , the aerodynamic roughness length, and/or the spatial distribution and size of roughness elements -through the use of a drag partitioning model (e.g., Raupach et al., 1993;Marticorena and Bergametti, 1995;Okin, 2008) that yields the stress exerted on the bare erodible soil.Equation ( 5) thus accounts for the effect of wind momentum absorption by nonerodible roughness elements on aeolian transport through the wind stress on the bare soil, as captured by the soil friction velocity u * .However, with the exception of Okin (2008), most previous studies have accounted for the effects of roughness elements by using the ratio of τ s / τ to scale the value of the threshold friction velocity u * t at which transport is initiated (Raupach et al., 1993;Marticorena and Bergametti, 1995).Although phenomenologically correct, the result of this approach is that, in the presence of nonerodible roughness elements, the quantity ρ a u 2 * overestimates the wind shear stress exerted on the bare soil.For instance, Marticorena and Bergametti (1995) equate the wind stress driving saltation to τ sand = ρ a u 2 * − u 2 * t (in their Eq.24), rather than (Owen, 1964), where the soil threshold friction velocity u * t is defined in more detail in the next paragraph.Therefore, using u * and u * t to parameterize saltation properties likely results in an overestimation of aeolian transport in the presence of nonerodible roughness elements (Webb et al., 2014), which our approach avoids.
In analogy to the threshold friction velocity u * t , the soil threshold friction velocity u * t is the minimum value of u * for which the bare soil experiences erosion.u * t depends on both the properties of the fluid and on the gravitational and interparticle cohesion forces that oppose the fluid lifting of sand particles that initiates saltation (Shao and Lu, 2000;Kok et al., 2012).In principle, u * t can be estimated from dust or sand flux measurements, as long as a correction is made for the presence of nonerodible elements, as discussed above and in the Supplement.However, the theoretical interpretation of this threshold is complicated by several factors.For instance, the threshold friction velocities at which saltation is initiated (the fluid or static threshold u * ft ) and terminated (the impact or dynamic threshold u * it ) are not equal.For most conditions, the impact threshold is thought to be smaller than the fluid threshold, of the order of ∼85 % (Bagnold, 1941;Kok, 2010).Moreover, spatial and temporal variations in soil conditions (Wiggs et al., 2004;Barchyn and Hugenholtz, 2011), as well as large variations in instantaneous wind speed for a given friction velocity (Rasmussen and Sorensen, 1999), make it such that there is generally not a clear value of u * above which saltation does occur and below which it does not (Wiggs et al., 2004).Despite these problems, we neglect here for simplicity the temporal and spatial variability of u * t and also assume that u * t = u * ft = u * it , as previous dust emission parameterizations have also done (e.g., Gillette and Passi, 1988;Shao et al., 1996;Marticorena and Bergametti, 1995).
In addition to u * t , we define the standardized threshold friction velocity (u * st ) as the value of u * t at standard atmospheric density at sea level (ρ a0 = 1.225 kg m −3 ).Consequently, u * st is not only independent of the presence of roughness elements, but is also invariant to variations in ρ a , and is thus equal for similar soils at different elevations.Therefore, u * st is a measure of the soil's susceptibility to wind erosion that depends on the state of the bare soil only.Since u * t ∝ √ ρ a (e.g., Bagnold, 1941), We hypothesize that u * st is a proxy for many of the soil properties known to affect dust emission, including soil cohesion, size distribution, and mineralogy (Fecan et al., 1999;Alfaro and Gomes, 2001;Shao, 2001).That is, although we do not understand in detail the effect of each of these soil properties on the dust flux (Shao, 2008), changes in soil properties that decrease the dust flux tend to also increase u * st .Consequently, it is possible that u * st can be used to partially account for the poorly understood effect of these soil properties on the dust flux.

The mean saltator impact speed (v imp )
After saltation has been initiated by the aerodynamic lifting of surface particles, new particles are brought into saltation primarily through the ejection, or splashing, of surface particles by impacting saltators (Ungar and Haff, 1987;Duran et al., 2011;Kok et al., 2012).(Note that this is only correct for soils with a sufficient supply of loose sand particles.The present theory is not valid for soils that instead are supplylimited, which we discuss in further detail in Sect.3.6.)Saltation is thus in steady state when exactly one particle is ejected from the soil bed for each particle impacting it.Since the number of splashed particles increases with the impacting saltator's speed (Kok et al., 2012), this condition for steady state is met at a particular value of v imp .Consequently, theory and measurements indicate that, while the shape of the probability distribution of v imp changes with u * (Fig. 1), v imp is independent of u * for steady-state saltation (Ungar and Haff, 1987;Duran et al., 2011;Kok, 2011a;Kok et al., 2012) (Supplement Fig. S1).Although v imp is independent of u * , it does depend on soil properties.In particular, the soil's saltation threshold sets the wind speed in the near-surface layer (Bagnold, 1941), which in turn determines the particle speed (Duran et al., 2011;Kok et al., 2012).Then, to first order, where C v ≈ 5 since v imp ≈ 1 m s −1 for loose sand with u * st ≈ 0.20 m s −1 (Supplement Fig. S1).

The fragmentation fraction (f frag )
An impacting saltator can fragment a dust aggregate in the soil if its impact energy exceeds a certain threshold (Kun and Herrmann, 1999;Kok, 2011b).The threshold impact energy per unit area ψ (J m −2 ) required to fragment a soil dust aggregate scales with the sum of the energetic cohesive bonds E coh between the constituent particles that make up the aggregate (Kun and Herrmann, 1999).That is, where D s is the saltator size, and the sum is over all interparticle bonds in the aggregate.Measurements and theory suggest that (Shao, 2001) where D c is the typical size of a constituent particle of the dust aggregate.The parameter β (J m −2 ) scales the interparticle force, which is the sum of a complex collection of individual forces, including van der Waals, water adsorption, and electrostatic forces (Shao and Lu, 2000).Consequently, β depends on the state of the soil, including soil moisture content, mineralogy, and size distribution.Since the number of bonds in the aggregate scales with D 3 ag D 3 c , where D ag is the aggregate size, Eq. ( 8) becomes For highly erodible, dry soils, β = β 0 ≈ 1.5 × 10 −4 J m −2 (Shao and Lu, 2000;Kok and Renno, 2006).Experiments suggest that most typical saltator impacts (i.e., D s = 100 µm and v imp = 1 m s −1 ) eject dust for such highly erodible, dry soils (Rice et al., 1996), yielding ψ 0 ≈ 0.1 J m −2 .Thus, where ψ = ψ / ψ 0 and β = β / β 0 .The dimensionless parameter c ψ is of order unity and depends on the soil size distribution since it scales with D 3 ag D 2 s D c .In particular, because saltators are often aggregates (Shao, 2008), with both D ag and D s having typical sizes of the order of 100 µm (Shao, 2001), the leading order scaling is likely c ψ ∼ D ag / D c .Here we take c ψ as a constant, both because there are insufficient vertical dust flux data sets available that report a detailed soil size distribution, and because global soil data sets are not nearly detailed enough to represent spatial and temporal variability in the soil size distribution.
Since the soil's standardized threshold friction velocity (u * st ) depends on the strength of interparticle forces (Shao and Lu, 2000), ψ must increase monotonically with u * st (Shao et al., 1996).This is intuitive: soils that are more erosion resistant, for example with strongly bound soil aggregates due to surface crusts or high moisture content, require a larger impact energy to fragment (Rice et al., 1996(Rice et al., , 1999)).For such soils, wind tunnel experiments show that only a small fraction of saltator impacts produce dust emission (Rice et al., 1996).
We calculate the fragmentation fraction f frag from the overlap between the probability distributions of ψ and the saltator impact energy per unit area E imp .Since ψ is the sum of a large number of individual cohesive bonds, its probability distribution P ψ (ψ) is normally distributed per the central limit theorem (Kallenberg, 1997), with a mean ψ and standard deviation σ ψ .The total fraction of saltator impacts that produces dust emission through fragmentation then equals where erf is the error function, which results from the integration of the normally distributed ψ.

Determining P E imp with the numerical saltation model COMSALT
In order to calculate f frag with Eq. ( 12), we require the probability distribution of saltator impact energies (P E imp ) for given values of u * , β, and D s , which we obtain through simulations with the numerical saltation model COMSALT (Kok and Renno, 2009).This model explicitly simulates the trajectories of saltators due to gravitational and fluid forces, and accounts for the stochasticity of individual particle trajectories due to turbulence and collisions with the irregular soil surface.Moreover, COMSALT simulates the retardation of the wind profile by the drag of saltating particles, which is the process that ultimately limits the number of particles that can be saltating at any given time.Finally, in contrast to many previous models, COMSALT includes a physically based parameterization of the ejection ("splashing") of surface particles, based on conservation of energy and momentum (Kok and Renno, 2009).Because of this explicit inclusion of splash, as well as other improvements over previous studies, COMSALT is the first numerical model capable of reproducing a wide range of measurements of naturally occurring saltation.
Since COMSALT was developed for saltation of soils made up of loose sand, it must be adapted in order to simulate saltation over dust-emitting soils.For soils made up of loose sand, the splashing of new saltating particles is constrained predominantly by the momentum transferred by impacting saltators (Kok and Renno, 2009).That is, the total momentum of splashed particles scales with the impacting saltator momentum (Beladjine et al., 2007;Oger et al., 2008).For dust emitting soils, this situation is likely different, because saltating particles are more strongly bound in the soil by co-hesive forces (Shao and Lu, 2000;Kok and Renno, 2009).We therefore assume that, for dust emitting soils, the number of particles splashed by an impacting saltator scales with its impacting energy (Shao and Li, 1999).Furthermore, in order for a saltating particle to eject another saltator from the soil, the impact must be sufficiently energetic to overcome the cohesive the bonds with other soil particles.Therefore, the larger the soil cohesive forces, the stronger the cohesive binding energy E coh,s with which sand-sized particles are bonded to other soil particles, resulting in a smaller number of splashed saltating particles N.That is, Since E coh,s scales with βD 2 s (see Eq. 9 and Shao, 2001), Eq. ( 13) becomes where ρ p ≈ 2650 kg m −3 is the density of the saltating particle (Kok et al., 2012), and the dimensionless parameter a E scales the number of splashed particles.We obtain a E = 6.1 × 10 −5 by forcing the minimum u * for which saltation can occur in COMSALT with β = β 0 to equal the minimal value of u * st for an optimally erodible soil.We define this minimal value as u * st0 , and measurements show that u * st0 ≈ 0.16 m s −1 for a bed of 100 µm loose sand particles (Bagnold, 1941;Iversen and White, 1982;Kok et al., 2012).
Other parameters of the splash process, such as the speed of splashed particles, the coefficient of restitution, and the probability that an impacting saltator does not rebound, are treated as described in Kok and Renno (2009).We thus neglect any change in these parameters with changes in soil cohesion since there is very little experimental data available to account for any such dependences (O'Brien and McKenna Neuman, 2012).COMSALT also computes the soil's standardized threshold friction velocity u * st as the minimum value of u * at which saltation can be sustained for a given value of β, following the procedure outlined in Kok and Renno (2009).
COMSALT simulations of P E imp show that, although the mean saltator impact speed (v imp ) remains approximately constant with u * (see above), the distribution of E imp does not (Fig. 1).Because the total drag exerted by saltators on the flow increases with u * , the wind profile lower in the saltation layer is relatively insensitive to u * (Owen, 1964;Ungar and Haff, 1987;Duran et al., 2011;Kok et al., 2012).Conversely, the wind speed higher up in the saltation layer does increase with u * (Bagnold, 1941), which causes the speed and abundance of energetic particles moving higher in the saltation layer to also increase.This causes a nonlinear increase in the high-energy tail of P E imp with u * (Fig. 1; also see Duran et al., 2011 andKok et al., 2012).).The value of f frag increases with u * for erosion-resistant soils, but not for highly erodible soils, as shown explicitly in (c).All plotted energy values are normalized by ψ 0 , the energy per unit area of a 100 µm saltator impacting at 1 m s −1 , and P ψ (ψ) was calculated using c ψ = 2 and σ ψ = 0.2ψ.

Dependence of f frag on u * and u * st
Since we can obtain P E imp for given values of u * , D s , and β (and thus u * st ) from COMSALT simulations, we can use Eq. ( 12) to determine f frag for given values of c ψ and σ ψ .Considering that the exact values of c ψ and σ ψ for any particular soil are unknown, our objective in using Eq. ( 12) is to understand the functional form of the dependence of f frag , and thus F d , on u * and u * st .To understand these dependencies, we consider the distributions of E imp and ψ for two limiting cases: a highly erodible and an erosion-resistant soil (Fig. 1).For a highly erodible soil, a large fraction of saltator impacts can be expected to produce fragmentation (Rice et al., 1996 and Fig. 1a), such that E imp ∼ ψ.In this case, the value of f frag is thus approximately constant with u * (Fig. 1c).Conversely, when the soil is erosion-resistant, E imp ψ, and only the high-energy tail of the impact energy distribution results in dust emission through fragmentation (Fig. 1b).Since this high-energy tail increases sharply with u * , f frag also increases sharply with u * (Fig. 1c).Consequently, F d scales more strongly with u * for erosion-resistant than for highly erodible soils.Our results thus show that f frag depends on both u * and u * st (Fig. 1c).Since f frag is dimen-sionless, its dependency on u * and u * st should take the form of the nondimensional ratios that capture the physical processes determining f frag (Buckingham, 1914).That is, f frag should depend only on (i) the dimensionless friction velocity u * / u * t , which sets the increase of the high-energy tail (Fig. 1), and (ii) the dimensionless standardized threshold velocity u * st / u * st0 , which sets the soil's susceptibility to wind erosion.From Fig. 1c, we infer Since this power law accounts for the dependence of f frag on u * / u * t , the dimensionless fragmentation constant C fr and exponent α must depend only on the other dimensionless number, u * st / u * st0 (Buckingham, 1914).Since highly erodible soils with u * st = u * st0 have α ≈ 0 (Fig. 1), we hypothesize that where C α is a dimensionless constant.Equation ( 16) is supported by numerical simulations of f frag for a range of plausi-  16) and ( 17) to the corresponding simulation results, and the solid black lines represents the best fit to the experimental data in Fig. 4.
ble values of the saltator diameter D s and the threshold fragmentation energy's normal distribution parameters (Fig. 2a).
The proportionality constant C fr in Eq. ( 15) must decrease sharply with u * st (Fig. 1c), because increases in u * st are primarily driven by increases in soil (aggregate) cohesion (Shao and Lu, 2000;Shao, 2008;Kok et al., 2012), for instance due to increases in soil moisture.Such increases in aggregate cohesion reduce the fragmentation fraction f frag , and numerical simulations indicate that (Fig. 2b) where C fr0 ≈ 0.5 is the fragmentation fraction for highly erodible soils (Fig. 1c), and C e is a dimensionless constant.

Full theoretical expression for the vertical dust flux
We complete our theoretical expression by substituting Eqs.
(2), (5), and ( 15)-( 17) into Eq.( 3), yielding where with C d0 = γ εC ns C fr0 C v .Equation ( 18) thus predicts that the dust flux (F d ) scales with the soil friction velocity (u * ) to the power a ≡ α + 2. We determine the dimensionless coefficients C α , C e , and C d0 through comparison against a quality-controlled compilation of vertical dust flux data sets in Sect.3. The dimensionless dust emission coefficient C d is independent of the soil friction velocity u * , and is thus a measure of a soil's ability to produce dust under a given wind stress.This susceptibility to dust emission is termed the soil erodibility in the dust modeling literature (e.g., Zender et al., 2003b), which is not to be confused with the identical term in the soil erosion literature referring more generally to the susceptibility of soil particles to detachment by erosive agents (e.g., Webb and Strong, 2011).
The increase in the dust emission coefficient C d with decreasing u * st accounts for a soil's increased ability to produce dust under saltation bombardment as the soil becomes more erodible (i.e., its threshold friction velocity decreases).This is an important result, as this process is not included in the previous dust flux parameterizations of Gillette and Passi (1988) and Marticorena and Bergametti (1995) that dominate dust modules in current climate models (e.g., Ginoux et al., 2001;Zender et al., 2003a;Huneeus et al., 2011).In particular, this result implies that the dust flux is more sensitive to the soil's threshold friction velocity than climate models currently account for.We further discuss this result and its implications in Sect. 4 and in the companion paper (Kok et al., 2014).
Note that the dust flux parameterization of Eq. ( 18) is considerably simpler than previous physically based dust emission models (Shao et al., 1996;Shao, 2001).This was achieved in large part by using u * st as a measure of soil erodibility, which allowed us to substantially simplify the energetics of dust emission.Furthermore, since our parameterization's main variables (u * , u * t , and f clay ) are available in weather and climate models, its implementation is relatively straightforward, in contrast to these more complex models (Darmenova et al., 2009).3 Assessment of parameterization performance using a quality-controlled compilation of dust flux measurements We test our proposed dust emission parameterization using a compilation of quality-controlled literature data sets.We do so by first separately testing the two main improvements of Eq. ( 18) over previous theories: the linear increase of the dust emission coefficient a with u * st , and the exponential decrease of the dust emission coefficient C d with u * st .This procedure also yields estimates of the dimensionless parameters C α , C d0 , and C e , subsequently allowing us to directly compare the measured dust flux against the predictions of Eq. ( 18).
The following section discusses the quality-control criteria that data sets need to meet in order to allow for an accurate comparison against our theoretical expression.Section 3.2 then describes the various corrections applied to bring all data sets on an equal footing, after which Sect. 3.3 describes the procedure for determining the dust emission coefficient (C d ) and fragmentation exponent (α) from literature data sets of dust flux measurements.We then test the functional form of the parameterization against the estimates of C d and α extracted from the literature data sets in Sect.3.4, and test the parameterization's predictions of the vertical dust flux against our dust flux compilation in Sect.3.5.Finally, we discuss the limitations of our parameterization in Sect.3.6.

Data set quality-control criteria
We strive to obtain a compilation of high-quality vertical dust flux measurements that we can use to test our new parameterization.We thus apply several quality-control criteria that data sets need to meet in order to be included in our compilation; these criteria are designed to ensure that the measured dust flux is governed by a soil in an approximately constant state.This is critical, because any changes in the soil state affects u * t , which is one of the main parameters in our parameterization.Since changes in the threshold friction velocity can occur on timescales as short as an hour (Wiggs et al., 2004;Barchyn and Hugenholtz, 2011), we only use data sets for which all data were taken within a limited time period of up to ∼12 h.This requirement excludes many of the data sets on which previous dust flux schemes were based, in particular data sets by Gillette (1979), Nickling and colleagues (Nickling, 1978(Nickling, , 1983;;Nickling and Gillies, 1993;Nickling et al., 1999), and Gomes et al. (2003).In addition, we require that a data set contains sufficient measurements to reliably determine the threshold friction velocity for the measurements.Furthermore, we only use data sets of natural dust emission taken in the field, because the characteristics of saltation and dust emission simulated in (portable) wind tunnels have been shown to, in some cases, be substantially different from the characteristics of natural saltation (Sherman and Farrell, 2008;Kok, 2011a).Finally, the measurements should be made for relatively homogeneous terrain, such that the soil state is spatially approximately constant.This last constraint is only required for predicting the dust emission coefficient C d .Therefore, data sets that meet all criteria except that of homogeneous terrain (i.e., the data sets of Fratini et al., 2007 andPark et al., 2011) are not used for comparison against the theoretical equations for C d and F d , but are still used for assessing the fragmentation exponent α.
Our literature search for vertical dust flux measurements that met the above quality-control criteria resulted in the identification of six studies: Gillies and Berkofsky, 2004 (hereinafter referred to as GB04), Zobeck and Van Pelt, 2006 (ZP06), Fratini et al., 2007 (FC07), Sow et al., 2009 (SA09), Shao et al., 2011 (SI11), and Park et al., 2011 (PP11).Images of the experimental sites of these six studies are shown in Fig. 3, and the main properties of each data set are summarized in Table 1.We used the original data for each of these six studies, and extracted 11 individual data sets from them.We describe the general procedures for correcting for differences between data sets and for extracting estimates of u * t , α, and C d in the next two sections.A detailed description of the analysis of each individual data set is provided in the Supplement.

Correcting for differences in averaging period and measured size range
A critical property of dust flux data sets is the time period over which measurements are averaged.In particular, since the vertical dust flux is nonlinear in the friction velocity, the averaging period needs to be consistent among the data sets (Sow et al., 2009;Martin et al., 2013).In setting the averaging period, an important consideration is that the friction velocity, being a turbulence parameter, is only meaningful when obtained over averaging periods long enough to sample a sufficient range of the turbulent eddies contributing to the downward flux of horizontal fluid momentum (Kaimal and Finnigan, 1994;Namikas et al., 2003;van Boxel et al., 2004).Moreover, the averaging period needs to be short enough such that the meteorological forcing of the boundary layer, which partially sets the downward momentum transfer, remains approximately constant.A compromise between these constraints is an averaging period of 30 min et al., 1996; Aubinet et al., 2001;van Boxel et al., 2004;Fratini et al., 2007), which conveniently is also of the order of the typical time step in global models.We thus reanalyzed each data set using a 30 min averaging period.In order to get maximum use out of each data set, the data were averaged over 30 min with a running average (e.g., a 60 min continuous data set with 1 min resolution yielded 31 data points).In addition to using the same averaging period for each data set, we also need to correct for differences in the measured dust size range between the data sets.We therefore corrected each data set to represent the mass flux of dust aerosols with a geometric diameter D d between 0 and 10 µm, which is a size range commonly represented in atmospheric circulation models (Mahowald et al., 2006).Several of the dust flux data sets (e.g., GB04, ZP06) reported size ranges not in terms of the geometric diameter D d , which is defined as the diameter of a sphere having the same volume as the irregularly shaped dust aerosol, but in terms of the aerodynamic diameter, D ae , which is defined as the diameter of a spherical particle with density ρ 0 = 1000 kg m −3 with the same aerodynamic resistance as the dust aerosol (Hinds, 1999).Therefore, depending on the data set, two separate corrections need to be made: one to correct from aerodynamic diameter to geometric diameter, and one to correct the measured geometric size range to 0-10 µm.
The geometric and aerodynamic diameters are related by Hinds (1999) and Reid et al. (2003) as where ρ p ≈ 2.5 ± 0.2 × 10 3 kg m −3 is the typical density of a dust aerosol particle (Kaaden et al., 2009), and χ is the dynamic shape factor, which is defined as the ratio of the drag force experienced by the irregular particle to the drag force experienced by a spherical particle with diameter D d (Hinds, 1999).Measurements of the dynamic shape factor for mineral dust particles with a geometric diameter of ∼10 µm find χ ≈ 1.4 ± 0.1 (Cartwright, 1962;Davies, 1979;Kaaden et al., 2009).Inserting this into Eq.( 19) then yields that D d ≈ (0.75 ± 0.04) D ae , where the standard error was obtained using error propagation (Bevington and Robinson, 2003).After converting each data set's measured aerodynamic particle size range to a geometric size range as necessary, we corrected the measured dust flux by assuming that the size distribution at emission is well-described by the theoretical dust size distribution expression of Kok (2011b) For instance, Eq. ( 6) in Kok (2011b) predicts that 71 ± 5 % of emitted dust in the geometric 0-10 µm size range lies in the aerodynamic 0-10 µm size range (which is equivalent to the geometric 0-7.5 ± 0.4 µm size range).We thus apply a correction factor of (0.71 ± 0.05) −1 = 1.42 ± 0.10 in order to correct a measured aerodynamic PM 10 flux (e.g., GB04, ZP06) to a geometric ≤ 10 µm flux.Note that the uncertainty in the correction factor is propagated into the uncertainty on the value of C d extracted from each data set (see the Supplement).
In addition to correcting for differences between data sets in the averaging time and the measured size range, we also corrected for differences in the fetch length when possible (see the Supplement).

Procedure for obtaining u * t , α, and C d
After putting all data on an equal footing using the above procedures, we extracted the parameters u * t , α, and C d from the dust flux data sets.Because u * t is required to determine the other parameters, we first determined the soil's threshold friction velocity for each data set.
Since many field experiments did not report the threshold friction velocity, and because of differences in the definition of threshold between data sets that did report a threshold friction velocity, we estimated u * t in a similar manner for each data set as described in detail in Sect.B in the Supplement.In brief, we estimated u * t using least-squares fitting of a second-order Taylor series of Eq. ( 22) below to saltation flux measurements within a limited range around the threshold (Barchyn and Hugenholtz, 2011).If the data set did not contain sand flux measurements, we instead used a least-squares fit of a second-order Taylor series of Eq. ( 18) to measurements of the dust flux.
After determining u * t in this manner, we used the following procedure to extract C d and α from each data set's dust flux measurements.Following Eq. ( 5), we start by calculating the dimensionless dust flux for each measurement of F d at given values of u * and u * t (obtained as described below) as Through substitution of Eq. ( 18) we now obtain an analytical expression for F d as a function of C d and α: We then use least-squares fitting of Eq. ( 21) to the values of F d calculated from dust flux measurements to determine the dust emission coefficient C d and the fragmentation exponent α, as well as their uncertainties, for each data set.The leastsquares fitting procedure and the calculation of uncertainties is described in more detail in the Supplement.
In addition, we obtain an independent estimate of the fragmentation exponent α, and thus the dust emission exponent a = α + 2, by using measurements of the sandblasting efficiency, which is defined as the ratio of the vertical dust flux to the horizontal saltation flux (Gillette, 1979).The sandblasting efficiency is thus defined for the data sets that reported measurements of both the dust flux and the (impact) flux of saltators at a certain height (i.e., ZP06, SA09, and SI11).This latter variable was usually measured with the Sensit piezoelectric instrument (Stockton and Gillette, 1990), which has been shown to provide a good measure of the horizontal saltation flux (Gillette et al., 1997;van Donk et al., 2003).
We extract α from measurements of the sandblasting efficiency as follows.We start with the saltation mass flux, which is given by (Bagnold, 1941;Kok et al., 2012) where L is the typical saltation hop length, and v is the average difference between the saltators' impact and lift-off speeds.The ratio L / v is thought to scale with the friction velocity, where the exponent r ranges from 0 (Ungar and Haff, 1987;Duran et al., 2011;Ho et al., 2011;Kok et al., 2012) to 1 (Owen, 1964;Shao et al., 1993), such that we take r = 0.5 ± 0.5.We now obtain an analytical expression for the sandblasting efficiency by combining equations (Eqs. 18,22,23): where the dimensional constant C s contains all parameters that do not depend on u * .We then obtain α and its uncertainty by fitting measurements of the sandblasting efficiency to the power law in u * of Eq. ( 24); this procedure is described in more detail in the Supplement.Note that an important advantage of the calculation of α from the sandblasting efficiency is that, unlike the calculation of α from the dimensionless dust flux described above, the result does not depend on the determination of the threshold friction velocity u * t .Therefore, errors that arise due to the procedure for assessing u * t do not affect the estimate of α derived from the sandblasting efficiency.

Test of the parameterization's functional form with dust flux measurements
All 11 data sets from the six studies that met the qualitycontrol criteria discussed in Sect.3.1 were used to determine the fragmentation exponent α through nonlinear leastsquares fitting of Eq. ( 21) to the vertical dust flux (see Supplement Fig. S5).Furthermore, five data sets featured simultaneous dust flux and saltation flux measurements, which we used to determine α by fitting Eq. ( 24) to the ratio of the vertical dust and horizontal saltation (impact) fluxes (see Supplement Fig. S6), and seven data sets were taken over spatially homogeneous terrain and thus were used to determine the dust emission coefficient C d (see Supplement Fig. S5).
The resulting analysis of the compilation of qualitycontrolled dust flux data sets shows an approximately linear increase in the dust emission exponent α with u * st (Fig. 4a), as predicted by Eq. ( 16).We obtain the dimensionless constant C α using least-squares fitting of Eq. ( 16), yielding C α = 2.7 ± 1.0.Moreover, the literature-extracted data sets show an approximately exponential decrease of the dust emission coefficient C d with u * st , as also predicted from our theory (Eq.18) and numerical simulations (Fig. 4b).We obtain C e = 2.0 ± 0.3 and C d0 = (4.4± 0.5) × 10 −5 from least squares fitting of Eq. (18b).

Test of the parameterization's predictions with dust flux measurements
After testing the parameterization's functional form and determining the values of its dimensionless coefficients, we can compare the predictions of Eq. ( 18) against our qualitycontrolled compilation of dust flux measurements.To avoid testing the model with the same data used to obtain its dimensionless coefficients (see previous section), we use the crosscorrelation method (e.g., Wilks, 2011;p. 252-253).That is, we use the following method for each data set: first, we obtain the dimensionless coefficients using the procedure in the previous section, but without using that particular data set or any other data sets from the same study.We then use the obtained dimensionless coefficients, which are thus specific for each of the six studies in our compilation, to predict the dust flux for each of the 11 data sets in our compilation.The resulting comparison between model and measurements is reported in Fig. 5c and Table 2.
For reference, we also compare against the predictions of the previous dust flux parameterizations GP88 (Gillette and Passi, 1988) and MB95 (Marticorena and Bergametti, 1995).Note that we unfortunately cannot compare our measurements compilation against the physically explicit dust flux parameterizations of Shao and colleagues (Shao et al., 1993(Shao et al., , 1996;;Shao, 2001), because these parameterizations use detailed soil properties that are unavailable for most data sets.
The MB95 dust flux parameterization is given by where the dimensionless parameter C MB is a proportionality constant, and the sandblasting efficiency η (units of m −1 ) depends on the clay fraction following η = 10 13.4f clay −6 .Note that Eq. ( 25) simplifies Eq. ( 34) in Marticorena and Bergametti (1995) by using a single value of u * t for the soil rather than different thresholds for different soil particle size bins.This is a common simplification necessary for the implementation of MB95 into most large-scale models (e.g., Zender et al., 2003a).Moreover, measurements, numerical models, and theory indicate that this simplification is actually more realistic (Bagnold, 1938;Rice et al., 1995;Namikas, 2006;Kok et al., 2012).Also, note that u * t in MB95 is calculated through a drag partition parameterization (Eq.20 in MB95), which we use for consistency for the comparison of MB95 against the measurement compilation (see the Supplement).The GP88 parameterization is given by where C GP (kg m −6 s 3 ) is a proportionality constant.Note that GP88 is thus formulated in terms of the soil friction velocity u * since it converts wind speed measurements taken  Gillette and Passi, 1988 (GP88), Marticorena and Bergametti, 1995 (MB95), and Eq. ( 18).RMSE values were calculated for two separate cases and the lowest RMSE for the three different parameterizations is underlined for each case.For each parameterization's first case, the proportionality constant was tuned to a single value that minimized the mean RMSE for all data sets.The resulting RMSE for this case is thus a measure of the parameterization's ability to reproduce variations in the dust flux due to variations in both u * and soil conditions (u * st and f clay ).For the second case, the proportionality constant in each parameterization was tuned separately for each data set.The resulting RMSE is thus a measure of a parameterization's ability to reproduce the dust flux's dependence on u * for each individual data set.Data set names are defined in Sect.3.1.
Our new parameterization reproduces the compilation of dust flux measurements with substantially less error than the parameterizations of GP88 and MB95 (Figs. 5a-c, S3 in the Supplement, Table 2).Equation ( 18) also produces better agreement when each parameterization's proportionality constant is tuned to each individual data set (Table 2, Fig. S2).

Limitations of the dust emission theory and parameterization
We derived the dust emission parameterization of Eq. ( 18) for dust emission occurring primarily through the fragmentation of either soil dust aggregates or saltating aggregates by the energetic impacts of saltators.Nonetheless, the main assumption used in deriving Eq. ( 18) is the existence of a normally distributed threshold controlling dust emission.Consequently, Eq. ( 18) theoretically applies to any dust emission processes controlled by an approximately normally distributed threshold.This point is underscored by the insensitivity of the functional form of Eqs. ( 16) and ( 17) to the threshold's normal distribution parameters and the saltator size (Fig. 2).Examples of dust emission processes other than aggregate fragmentation that are controlled by a normally distributed threshold could include dust emission from crusted soils (Rice et al., 1996) and from sand particles with clay coatings (Bullard et al., 2004).Since we do not know what the relative contribution of different dust emission pro-cesses is to each of the dust flux data sets used to calibrate the dimensionless coefficients in Eq. ( 18), it is likely that the obtained values of these coefficients represents some weighted average of the relative contribution of each dust emission process.As discussed in Sect.2.2, we consider it most likely that the fragmentation process contributes the largest fraction of the dust flux for each data set.Thus, although our parameterization theoretically applies to dust emission from soils dominated by processes other than fragmentation, the dimensionless coefficients in Eq. ( 18) could be quite different for such soils.We are not aware of any experimental data sets that meet our quality-control criteria that could be used to estimate the dimensionless coefficients for soils for which dust emission is dominated by any specific process other than fragmentation.Furthermore, as mentioned in Sect.2.2.1, our theory applies only to soils for which the saltation flux is limited by the availability of wind momentum, and are thus transport limited (e.g., Nickling and McKenna Neuman, 2009).The present theory is thus not valid for soils for which the horizontal saltation flux at a given point in time is limited by the availability of sand-sized sediment.Such supply-limited soils are inherently inefficient sources of dust aerosols (Rice et al., 1996), and are thus probably less important in the global dust budget.Note that dust emission from some prominent sources can be limited by the sediments supplied to these sources, for instance through the deposition of fluvially eroded sediment (Bullard et al., 2011;Ginoux et al., 2012).However, when substantial emission occurs from such regions, the soil is generally not supply limited at that point in time (Bullard et al., 2011), such that Eq. ( 18) could be used to parameterize the dust flux.
Our parameterization attempts to include only the most important processes affecting the dust flux.Thus, Eq. ( 18) does not explicitly account for many other processes that might affect dust emission, including changes in the parameters γ and ε with u * and u * t , and the dependence of c ψ and σ ψ on the soil size distribution, mineralogy, and other soil properties.Future studies should consider these effects, especially if more extensive global (or regional) soil data sets become available, or if more dust flux data sets that sufficiently characterize these soil properties become available.However, as mentioned above, many of these processes partially affect the dust emission flux F d by increasing or decreasing u * st , such that some of their effect might be captured in the calibration of the dimensionless coefficients of Eq. ( 18) to our compilation of vertical dust flux data sets.
Another limitation of our theory is that it does not account for dust emission due to saltator impacts that do not produce fragmentation but that nonetheless produce dust by "damaging" the dust aggregate (Kun and Herrmann, 1999).It also does not account for the lowering of an aggregate's fragmentation threshold through the rupturing of cohesive bonds by impacting saltators.These effects might dominate for very erosion-resistant soils, such as crusted soils.A further limitation of our theory is that it simplifies the energetics of dust emission by considering u * st the prime determinant of soil erodibility (Shao and Lu, 2000).Although the threshold for saltation (u * st ) and the threshold energy required to fragment dust aggregates (ψ) are likely strongly coupled for many soils (Shao et al., 1993;Rice et al., 1996Rice et al., , 1999)), increases in ψ might not produce corresponding increases in u * st for some soils.An example of such a soil is a sandy soil for which dust emissions occur primarily from the removal of dust coatings on sand grains (Bullard et al., 2004); thus, emission from such soils might be poorly captured by the present theory.
We have used a combination of theory and numerical simulations to derive a physically based parameterization of the vertical dust flux emitted by an eroding soil.Our new dust flux parameterization includes two main improvements over previous schemes used in large-scale models.First, it accounts for the predicted (Figs. 1, 2a) and observed (Fig. 4a) increasing scaling of F d with u * that occurs with increasing threshold friction velocity t ; this advance helps explain the numerous observed scalings of F d with u * (Shao, 2008;Kok et al., 2012).Second, our parameterization accounts for a soil's increased ability to produce dust under saltation bombardment as the soil becomes more erodible (Figs. 1,2b,4b).This second improvement is especially important, as it implies that previous parameterizations have underestimated the sensitivity of the dust flux to the soil's dust emission threshold (u * st ) (also see Fig. 1 in Kok et al., 2014).This underestimation is not sensitive to the details of our parameterization because it follows directly from the energetics of dust emission: increases in soil cohesion both raise the dust emission threshold and cause dust emission to require more energy, thereby reducing the dust flux for a given saltator kinetic impact energy.Previous work by Shao and colleagues (Shao et al., 1993(Shao et al., , 1996;;Shao, 2001) has noted that soils with stronger interparticle forces should produce less dust per saltator impact, but this insight had not been included in dust emission parameterizations commonly implemented in large-scale models (e.g., Ginoux et al., 2001;Zender et al., 2003a;Cakmur et al., 2006;Menut et al., 2013;Zhao et al., 2013).
Partially as a result of the inclusion of these two additional physical processes, our parameterization is in better agreement with a quality-controlled compilation of dust flux measurements than the previous dust flux parameterizations of Gillette and Passi (1988) and Marticorena and Bergametti (1995) (see Fig. 5).Although our parameterization thus appears to account for more of the processes driving the dust flux than these previous parameterizations, it is straightforward to implement as it uses only variables that are readily available in weather and climate models (note that the code to implement the parameterization in the Community Earth System Model is freely available from the main author).This is made possible because of several advances and simplifications over previous theories.Arguably the main advance is that we use the soil's standardized threshold friction velocity (u * st ) as a measure of soil erodibility (i.e., the soil's ability to emit dust), allowing us to substantially simplify the energetics of dust emission relative to previous physically explicit schemes (Shao et al., 1996;Shao, 2001).Furthermore, many previous parameterizations used a different threshold friction velocity for each soil particle size bin.However, experiments, numerical modeling, and theory all indicate that, once the saltation threshold is exceeded, particles of a wide range of sizes are set into motion (e.g., Bagnold, 1938;Rice et al., 1995;Namikas, 2006;Kok and Renno, 2009;Kok et al., 2012).We therefore characterized the threshold friction velocity with a single value, which can for instance be calculated using the models of Iversen and White (1982), Fecan et al. (1999), or Shao and Lu (2000).
Our result that the dust flux is more sensitive to the soil's threshold friction velocity than most current parameterizations account for emphasizes the importance for models to accurately represent spatial and temporal variations in soil erodibility.Our parameterization provides a convenient way of doing so through the exponential dependence of the dust emission coefficient C d on the standardized dust emission threshold u * st .However, the parameterization of u * st in most models is relatively primitive (e.g., Zender et al., 2003a).For instance, one of the main determinants of u * st is the moisture content of the top layer of soil particles.Yet, the most commonly used parameterization of the effect of soil moisture on u * st (Fecan et al., 1999) is found to produce unrealistic results in some models, requiring the use of a tuning constant (Zender et al., 2003a;Mokhtari et al., 2012).Furthermore, effects of soil aggregation and crust formation on u * t are not included in the most widely used global dust modules (Ginoux et al., 2001;Zender et al., 2003a;Huneeus et al., 2011).Considering the paramount importance of u * st in determining dust fluxes (see Eq. 18), an effective way to improve the fidelity of dust cycle simulations would be to develop improved parameterizations of u * st as a function of soil properties, precipitation events, atmospheric relative humidity, and other relevant parameters.Alternatively, for simulations of the current dust cycle, u * st could be remotely sensed (Chomette et al., 1999;Chappell et al., 2005;Draxler et al., 2010).Doing so requires the simultaneous determination of the threshold wind speed and the surface roughness (Marticorena et al., 2004), such that the remotely sensed threshold wind stress can be partitioned between the portion causing dust emission (τ s ) and that absorbed by nonerodible elements (τ R ) (Raupach et al., 1993;Marticorena and Bergametti, 1995).
Current large-scale models commonly use semiempirical dust source functions (e.g., Ginoux et al., 2001;Tegen et al., 2002;Zender et al., 2003b) to help parameterize spatial variability in soil erodibility and the consequent dust emissions.The use of these source functions usually shift emissions towards the most erodible regions.Because our parameterization accounts for a soil's increased ability to produce dust under saltation bombardment as the soil becomes more erodible, its implementation in models would also result in a shift of emissions to the most erodible regions.We therefore hypothesize that our parameterization reduces the need for empirical source functions in dust modules.We test this hypothesis in our companion paper (Kok et al., 2014).probability distribution of threshold fragmentation energy ψ, J −1 m 2 ρ 0 standard density of aerosol particle, kg m −3 Q saltation mass flux, kg m −1 s −1 r exponent of u * scaling the ratio of saltating particle hop length and impact and rebound speed differential ρ a air density, kg m −3 ρ a0 air density at standard atmosphere, kg m −3 ρ d density of a dust aerosol particle, kg m −3 ρ p density of a saltating particle, kg m −3 σ ψ standard deviation of normal distribution of the threshold impact energy ψ, J m −2 τ total wind stress exerted on surface, N m −2 τ R wind stress exerted on nonerodible roughness elements only, N m −2 τ s wind stress exerted on bare soil only, N m −2 τ st threshold wind stress exerted on bare soil above which saltation occurs, N m −2 u * soil friction velocity, derived from shear stress on bare soil τ s , m s −1 u * friction velocity, derived from total shear stress on surface τ , m s −1 u * st the soil threshold friction velocity standardized to standard atmospheric density, m s −1 u * st0 the standardized threshold friction velocity of an optimally erodible soil, m s −1 u * t soil threshold soil friction velocity above which saltation occurs, m s −1 u * t threshold friction velocity above which saltation occurs, m s −1 v imp mean saltator impact speed, m s −1 χ dynamic shape factor for irregularly shaped aerosol particles ψ threshold impact energy per unit area required to fragment a soil dust aggregate, J m −2 ψ mean value of normal distribution of the threshold impact energy ψ, J m −2

Figure 1 .
Figure 1.Probability distributions of the threshold impact energy per unit area (P ψ ) required for aggregate fragmentation (solid black line), and of the saltator impact energy per unit area (P E imp ) for saltation of 100 µm particles at different values of u * (colored lines).Shown are results for (a) a highly erodible soil (u * st = 0.16 m s −1 ) and (b) an erosion-resistant soil (u * st = 0.40 m s −1).The value of f frag increases with u * for erosion-resistant soils, but not for highly erodible soils, as shown explicitly in (c).All plotted energy values are normalized by ψ 0 , the energy per unit area of a 100 µm saltator impacting at 1 m s −1 , and P ψ (ψ) was calculated using c ψ = 2 and σ ψ = 0.2ψ.

Figure 2 .
Figure 2. Simulations of the fragmentation exponent α (a) and fragmentation constant C fr (b) with the numerical saltation model COMSALT (Kok and Renno, 2009) for different values of the saltating particle size (D s ) and the threshold fragmentation energy's normal distribution parameters (c ψ and σ ψ ).The colored dashed lines represent the best fits of the functional forms of Eqs.(16) and (17) to the corresponding simulation results, and the solid black lines represents the best fit to the experimental data in Fig.4.

Figure 4 .
Figure 4. Values of (a) the dust emission exponent a (= α + 2) and (b) the dust emission coefficient C d as a function of the standardized threshold friction velocity u * st , determined from the analysis of available quality-controlled data sets.Open symbols refer to estimates of C d and a from the least-squares fit of the measured dust flux to Eq. (18), whereas filled symbols refer to estimates of a from a least-squares fit to ratios of the measured vertical dust flux and the horizontal saltation flux (see text for details).The dashed line indicates the best-fit forms of Eqs.(16) and (18b), and the grey shaded area denotes one standard error from the fitted relation.Data set names are defined in Sect.3.1.

Figure 5 .
Figure 5.Comparison of measured dust fluxes with the predictions of the parameterizations of (a) Gillette and Passi (1988), (b) Marticorena and Bergametti (1995), and (c) this study.The proportionality constant in each parameterization was adjusted to maximize agreement with the compilation of measurements.To prevent cluttering of the graph, only 15 representative measurements are shown for each data set.Error bars denote uncertainty arising from the measurement of u * t , u * , and F d (see the Supplement).Data set names are defined in Sect.3.1.

Table 1 .
Summary of the main characteristics of the quality-controlled data sets used in this study.Data set names are defined in Sect.3.1.

Table 2 .
Root mean square error (RMSE) of the vertical dust flux predicted by the parameterizations of

Table A1 .
Notation.energeticbond between constituent particles in dust aggregate (E coh ), J m −2 β 0 approximate value of β for an optimally erodible soil, J m −2 c ψ dimensionless constant linking ψ to β C d dimensionless dust emission coefficient, scaling the vertical dust flux C d0 dimensionless constant scaling the dust emission coefficient C e dimensionless constant scaling the exponential decrease of C d with u * st C fr dimensionless constant scaling the fragmentation fraction (f frag ) C fr0 value of C fr for an optimally erodible soil C GP constant scaling the dust flux in GP88, kg m −6 s 3 C MB dimensionless constant scaling the dust flux in MB95 C ns dimensionless constant scaling the number of saltator impacts (n s ) C s dimensional constant scaling the ratio of vertical dust flux to horizontal saltation flux C v dimensionless constant scaling the mean saltator impact speed (v imp ) C α dimensionless constant scaling the fragmentation exponent (α) γ dimensionless constant scaling the emitted dust per saltator fragmenting impact (m frag ) D ae dust aerosol aerodynamic diameter, m D ag size of soil dust aggregate, m D c typical size of constituent particles of a soil dust aggregate, m D d dust aerosol geometric diameter, m D s size of saltating particle, m vaverage difference between the saltators' impact and lift-off speeds E coh energy of the energetic bond between constituent particles of soil dust aggregate, J E coh,s energy of the energetic bond between sand particles and other soil particles, J ε mass fraction of emitted dust that does not reattach to the surface f bare fraction of the surface consisting of bare soil impacts on soil surface per unit area and time, m −2 s −1 N number of particles splashed by impacting saltator P E imp probability distribution of saltator impact energy E imp , J −1 P ψ