Articles | Volume 23, issue 11
Research article
14 Jun 2023
Research article |  | 14 Jun 2023

A new process-based and scale-aware desert dust emission scheme for global climate models – Part I: Description and evaluation against inverse modeling emissions

Danny M. Leung, Jasper F. Kok, Longlei Li, Gregory S. Okin, Catherine Prigent, Martina Klose, Carlos Pérez García-Pando, Laurent Menut, Natalie M. Mahowald, David M. Lawrence, and Marcelo Chamecki

Desert dust accounts for most of the atmosphere's aerosol burden by mass and produces numerous important impacts on the Earth system. However, current global climate models (GCMs) and land-surface models (LSMs) struggle to accurately represent key dust emission processes, in part because of inadequate representations of soil particle sizes that affect the dust emission threshold, surface roughness elements that absorb wind momentum, and boundary-layer characteristics that control wind fluctuations. Furthermore, because dust emission is driven by small-scale ( 1 km or smaller) processes, simulating the global cycle of desert dust in GCMs with coarse horizontal resolutions ( 100 km) presents a fundamental challenge. This representation problem is exacerbated by dust emission fluxes scaling nonlinearly with wind speed above a threshold wind speed that is sensitive to land-surface characteristics. Here, we address these fundamental problems underlying the simulation of dust emissions in GCMs and LSMs by developing improved descriptions of (1) the effect of soil texture on the dust emission threshold, (2) the effects of nonerodible roughness elements (both rocks and green vegetation) on the surface wind stress, and (3) the effects of boundary-layer turbulence on driving intermittent dust emissions. We then use the resulting revised dust emission parameterization to simulate global dust emissions in a standalone model forced by reanalysis meteorology and land-surface fields. We further propose (4) a simple methodology to rescale lower-resolution dust emission simulations to match the spatial variability of higher-resolution emission simulations in GCMs. The resulting dust emission simulation shows substantially improved agreement against regional dust emissions observationally constrained by inverse modeling. We thus find that our revised dust emission parameterization can substantially improve dust emission simulations in GCMs and LSMs.

1 Introduction

Desert dust accounts for more than half of the atmospheric mass loading of particulate matter (PM) (Kinne et al., 2006; Kok et al., 2017) and produces a wide range of important impacts on multiple components of the Earth system (Shao et al., 2011a; Kok et al., 2023). Like other aerosols, dust changes Earth's radiative budget and atmospheric dynamics directly by scattering and absorbing radiation (Sokolik and Toon, 1996; Miller and Tegen, 1998) and indirectly by mediating cloud formation (Rosenfeld et al., 2001; Shi and Liu, 2019; McGraw et al., 2020; Froyd et al., 2022). These dust–radiation interactions and dust–cloud interactions also drive day-to-day variability in large-scale circulation patterns and local weather events such as monsoons and rainfall (Jin et al., 2021; Parajuli et al., 2022). Dust further impacts biogeochemistry by delivering nutrients such as iron and phosphorus to ocean and land ecosystems (Mahowald et al., 2010; Hamilton et al., 2020).

In recent decades, modelers have made substantial progress in developing various parameterizations for the main dust cycle processes including emission (e.g., Shao et al., 1993; Marticorena and Bergametti, 1995; Marticorena et al., 1997; Tegen and Fung, 1995; Klose and Shao, 2013), advection (e.g., Prospero, 1999; Lin, 2004; Van der Does et al., 2018), deposition (e.g., Barth et al., 2000; Liu et al., 2001; Zhang et al., 2001; Petroff and Zhang, 2010) and its biogeochemical effects (e.g., Jickells et al., 2005; Mahowald et al., 2005, 2010; Hamilton et al., 2020), optics (e.g., Sokolik and Golitsyn, 1993; Linke et al., 2006; Adebiyi and Kok, 2020), and radiative effects (e.g., Di Biagio et al., 2020; Li et al., 2021). Despite substantial progress in dust modeling, current global climate models (GCMs) and Earth system models (ESMs) still struggle to adequately simulate the dust cycle, which impedes an accurate assessment of dust impacts (Huneeus et al., 2011; Wu et al., 2020; Zhao et al., 2022). For instance, model simulations still show large discrepancies when compared against observations of the spatial and temporal characteristics of the dust cycle, including dust emission (Kok et al., 2014a; Pierre et al., 2014a), dust PM concentration (Wu et al., 2019; Pu et al., 2020; Li et al., 2022), dust aerosol optical depth (DAOD or DOD) (Ridley et al., 2012; Kok et al., 2014b; Pu and Ginoux, 2018; Parajuli et al., 2019), dust deposition (Ginoux et al., 2001; Albani et al., 2014; Kok et al., 2014b; Li et al., 2022), and dust size distributions (Parajuli et al., 2019; Adebiyi and Kok, 2020; Li et al., 2022). Also, models struggle to capture the observed interannual and decadal variability of dust (Ridley et al., 2014; Smith et al., 2017; Evan, 2018; Kok et al., 2018) and the sensitivity of dust to climate changes (Evan, 2018; Kok et al., 2018). An improved quantification of dust impacts on the Earth system thus requires improvements on how dust is simulated in models.

One key piece of physics that models struggle to parameterize is the dust emission threshold. The dust emission threshold u∗t is defined as the threshold wind stress/speed above which winds initiate, or below which winds cease, the lifting of sand particles whose impacts on the soil surface emit dust aerosols (Kok et al., 2012; Comola et al., 2019b). The dust emission threshold is a function of soil properties and atmospheric conditions like particle size distribution, soil moisture, and air density. There are various reasons for the inadequate parameterization of the dust emission threshold. First, many models assume a globally constant soil particle size in calculating a spatially varying dust emission threshold (Zender et al., 2003a; Darmenova et al., 2009; Kok et al., 2014b), whereas the actual soil particle size is likely a function of space and time and could depend on soil properties, such as texture, pH, and organic matter content, since these variables modulate the cohesion between soil particles (Webb et al., 2016). Some models proposed that soil particle sizes are related to the soil texture and therefore represent the soil particle size as a global map (Tegen et al., 2002; Darmenova et al., 2009; Menut et al., 2013; Klose et al., 2021), but these maps either have not yet been thoroughly validated against observations or are based upon extrapolation of a limited number of observations. Second, most current models use the fluid threshold (also named static threshold or initiation threshold) above which saltation is initiated as the dust emission threshold, but it is well known that dust emission is governed by both the larger fluid threshold and the smaller impact threshold (also named dynamic threshold or cessation threshold) below which saltation is terminated (Bagnold, 1941; Shao, 2008; Martin and Kok, 2018; Comola et al., 2019a, b; Pähtz et al., 2020). Moreover, dust emission is a nonlinear process (i.e., it varies with the wind speed to the second to fifth power, per Kok et al., 2014a), and the emission flux is particularly sensitive to the magnitude of the emission threshold (Kawai et al., 2021). Thus, land-surface models (LSMs) within GCMs and ESMs need to parameterize the emission threshold correctly to get an adequate spatiotemporal variability of the modeled atmospheric dust.

The second key dust emission physics that LSMs struggle to represent is the partitioning of the wind stress. Wind drag is partitioned into the part absorbed by surface roughness elements (mainly rocks and plants) and the part exerted on the bare soil that drives dust emissions. This drag partitioning effect is modeled by several dynamical schemes (Raupach et al., 1993; Marticorena and Bergametti, 1995; Okin, 2008), and it is accounted for in some models (LeGrand et al., 2019; Klose et al., 2021; Tai et al., 2021) but not others (e.g., Kok et al., 2014b; Evans et al., 2016). One major challenge in modeling drag partition is to quantify the abundance of rocks (which includes rocks, pebbles, and gravel in this study) and their corresponding partition effect, because there are few measurements of rock roughness. To cope with this issue, studies have used in situ and/or remote sensing scatterometer measurements to quantify the small-scale land-surface roughness (e.g., Greeley et al., 1997; Roujean et al., 1997; Marticorena et al., 2004; Prigent et al., 2005, 2012), especially over arid desert regions over which rocks, pebbles, and gravel dominate the roughness. However, with a few recent notable exceptions that attempted to represent the roughness effects of both rocks and vegetation (e.g., Darmenova et al., 2009; Foroutan et al., 2017; Klose et al., 2021), studies often omitted the drag partition effect either due to vegetation (e.g., Menut et al., 2013) or due to rocks (e.g., Wu et al., 2016; LeGrand et al., 2019; Tai et al., 2021). To resolve these issues, we propose a new approach that combines the drag partition effects of both elements, leveraging satellite scatterometer measurements to quantify the surface rock roughness and using observable vegetation and land-surface variables to quantify the surface vegetation roughness.

The third key piece of fundamental dust emission physics not accounted for by many models is the effect of turbulence-driven high-frequency wind variability on dust emissions. Most current GCMs assume a constant wind speed (and thus a constant emission flux) within the relatively large model time step, e.g., 30 min (e.g., Rahimi et al., 2019; Dunne et al., 2020). However, in reality, high-frequency turbulent fluctuations cause the wind speed to fluctuate within a time step (from seconds to minutes). Because dust emissions scale nonlinearly with wind speed, this causes highly uneven and fluctuating dust emission fluxes (Durán et al., 2011). Even more importantly, turbulent wind fluctuations can sweep across the dust emission threshold multiple times and shut off dust emissions intermittently within one model time step, resulting in strong dust emission intermittency (Comola et al., 2019b). Even regional climate models (RCMs), which typically use a smaller time step (e.g., <1 min), do not resolve turbulence unless they are run in the computationally expensive large-eddy-simulation (LES) mode (e.g., WRF–LES). Omitting turbulence by GCMs and RCMs thus causes either an overestimate or an underestimate of dust emissions, especially over marginal source regions where winds fluctuate around the high emission threshold, as models do not account for the cessations or initiations of dust emissions due to turbulent fluctuations. To account for the instantaneous wind fluctuations, a dynamical approach is to derive a probability density function (PDF) for the instantaneous momentum flux using LES, which is then used for quantifying instantaneous dust emission fluctuations (Klose and Shao, 2012; Klose et al., 2014). A parameterization approach is to use the Monin–Obukhov similarity theory (MOST) to relate the standard deviation of the instantaneous wind to the boundary-layer dynamical variables (Comola et al., 2019b). In this study, we will account for turbulent dust emissions by following Comola et al. (2019b), which showed significant improvements in representing the small-magnitude saltation and dust fluxes that are particularly important over marginal source regions.

In addition to these issues of models missing some of the fundamental physics of dust emission, a central issue in modeling the global dust cycle is that dust emissions are grid-resolution-dependent because of the nonlinear dependence of dust emissions to meteorological fields and land-surface variables. Since dust emission varies nonlinearly with wind speed and has an even more complex relation to the soil moisture (Gillette and Passi, 1988; Fécan et al., 1999; Shao, 2001; Kok et al., 2014a), the total regional and global emissions can vary significantly with grid resolution (Ridley et al., 2013; Meng et al., 2021). For instance, modeled emissions were found to increase by  29 % from a 1× 1 to 0.25× 0.25 resolution (Ridley et al., 2013; Feng et al., 2022). As a consequence, GCMs and ESMs often need to tune emissions separately for different grid resolutions to match observational dust budgets (Ginoux et al., 2001; Zender et al., 2003a; Albani et al., 2014; Kok et al., 2014b; Chappell et al., 2021). This issue occurs because current GCM grid sizes of  1 or 100 km cannot resolve the spatial scales of  1 m to  1 km over which soil properties and wind speeds change (Ridley et al., 2013). When adopting coarse grid resolutions, coarser modeled meteorological fields (for GCMs) or spatially averaged input meteorological fields (for chemical transport models or CTMs) will smooth out the local wind extrema, possibly causing wind speeds to fall below the dust emission threshold. As a result, the coarse modeled winds usually result in strong GCM emissions underestimations (Ridley et al., 2013). The same smoothing problem also occurs for soil moisture for instance, with its maxima smoothed out leading to an overestimation of dust emissions. Thus, although some GCMs and ESMs recently implemented more physical schemes (Zhao et al., 2022), their inability to resolve the small scales still causes challenges for capturing the accurate spatial distributions of dust emissions (Meng et al., 2021). For the same reason, GCMs tend to neglect small-scale emissions over marginal source regions. In this study, we will analyze the scale dependence of our dust emission scheme given specific input datasets and propose a method of upscaling the coarse dust emissions to alleviate the scale-dependence problem.

To tackle the above problems and improve simulations of the global dust cycle, we propose a new emission scheme for global models that includes key dust emission physics missing from current models. Specifically, we (1) account for the effects of the soil particle size distribution (PSD) on the dust emission threshold, (2) draw on satellite data and physically explicit models to account for wind momentum absorption by both rocks and vegetation, and (3) account for turbulence-driven intermittency in dust emission fluxes. After we review current dust emission schemes in Sect. 2, we present our new scheme in Sect. 3. In Sect. 4, we code the new dust emission scheme as a standalone sandbox model (see Sect. 2.4) and examine the resulting spatiotemporal variability of the new dust emissions. We then examine the grid-scale dependence of dust emissions and derive a correction map for coarser dust emission simulations to correct their spatial variability to high-resolution simulations. Sections 5 and 6 discuss and summarize the main findings of this paper. In our companion paper (Leung et al., 2023b), we will implement this new scheme in the state-of-the-art Community Earth System Model version 2 (CESM2) and evaluates its performance against observations.

2 Current dust emission schemes and their input variables

This section provides a basic review of current GCM's dust emission schemes and their required input meteorological variables. Broadly speaking, dust emission schemes consist of a parameterization of the dust emission threshold (Sect. 2.1); a parameterization of the reduction of the wind drag on the bare soil surface due to momentum absorption by surface roughness elements (Sect. 2.2); and a parameterization of dust emission flux given wind speed, threshold wind speed, and drag absorption (Sect. 2.3). We will develop an improved emission scheme for GCMs in Sect. 3 by improving upon each of these three core ingredients of dust emission schemes. We also provide a brief description of meteorological data used for computing the dust emission schemes in Sect. 2.4.

2.1 Parameterization of the dust emission thresholds

There have been extensive studies on the dust emission threshold u∗t, defined as the wind speed or drag that corresponds to the initiation or cessation of dust emission. Dust emission is caused by saltation, a process by which sand particles (geometric diameter >63µm) on the surface are lifted by wind drag into the airstream (known as aerodynamic entrainment) and undergo ballistic trajectories (Anderson, 1989; Kok et al., 2012). The minimum wind friction velocity u required for initiating saltation through aerodynamic entrainment is called the fluid threshold u∗ft (McKenna Neuman and Sanderson, 2008). Once saltation is initiated, a smaller u is needed to maintain saltation because saltation bombardment (saltating particles impacting on the granular bed) can create further saltation, which is more efficient than creating saltation solely through the wind drag. Saltation will be maintained at a slightly smaller u called the impact threshold u∗it (Bagnold, 1941; Martin and Kok, 2018). The ratio uit/uft is about 0.8–0.85 for loose dry sand and less for soils with other sources of cohesion (e.g., moisture, organic matter), because cohesion rapidly increases u∗ft, but low-to-moderate levels of cohesion do not increase u∗it as indicated by numerical simulations (Comola et al., 2019a; Ralaiarisoa et al., 2022). It follows that dust emission can occur below u∗ft, especially in marginal dust source regions with high soil moisture for which u∗it can be much smaller than u∗ft. u∗t is thus a general concept comprised of both u∗ft and u∗it. However, u∗it is not currently accounted for in most current GCMs, which simply use u∗ft as u∗t.

One challenge in parameterizing u∗ft and u∗it lies in the representation of the effect of the soil particle size Dp on both thresholds. These two thresholds are mainly governed by soil particle diameter Dp, air density ρa, and soil moisture w (Greeley et al., 1997; Shao and Lu, 2000). Although there are multiple data sources of globally gridded products for ρa and w, there are relatively few efforts on obtaining globally gridded Dp, since there are no methods for satellites to observe and derive surface Dp observations. With few comprehensive field studies of saltation dynamics over polydisperse soils, past saltation studies either assumed that particles of different sizes saltate independently of each other (Marticorena and Bergametti, 1995; Shao et al., 1996; Alfaro and Gomes, 2001; Zender et al., 2003a) or assumed that a single grain size (e.g., the median) could be used to represent the whole PSD of the soil bed (Elbelrhiti et al., 2005; Andreotti et al., 2010). The dispute of whether the assumption of “independent” saltation (Shao, 2008) or “representative” saltation (Claudin and Andreotti, 2006) is more appropriate was informed by Martin and Kok (2019), who showed that modeling the threshold using a single particle size (representative saltation) is more realistic than assuming no interactions between saltation of different particle sizes. They argued that the median particle diameter Dp of the PSD should be used for calculating the threshold of a mixed soil, and the emission should then be calculated for the whole soil bed using the median instead of using a spectral/independent approach of calculating emissions of different particle sizes separately.

To model u∗ft, current models assume that u∗ft is mainly dependent on the soil PSD and soil moisture (Iversen and White, 1982; Marticorena and Bergametti, 1995; Zender et al., 2003a):

(1) u ft = u ft 0 ( D p , ρ a ) f m ( w ) ,

where u∗ft0 is the “dry” fluid threshold friction velocity (in m s−1) on a smooth and bare surface as a function of air density ρa(long,lat,t) (longitude, latitude, time) and Dp(long,lat), which in this study will be the median diameter Dp of a polydispersed, mixed soil PSD. fm is the correction factor for the presence of soil moisture w(long,lat,t); fm≥1, such that soil moisture protects soil particles from being lifted. u∗ft is the “wet” fluid threshold accounting for the moisture effect. Other factors can also affect u∗ft, such as salt concentration, organic matter, electrostatic effects, and surface crusts, but they are not included in most studies because they are not well understood and modeled (Shao et al., 2011; Foroutan et al., 2017).

u∗ft0 is parameterized by considering the balance between aerodynamic drag and lift against gravity and interparticle cohesion on a soil particle. The Shao and Lu (2000) (hereafter S&L00) scheme derived a simple solution to the force balance (also see Kok et al., 2012), assuming that the cohesive force is proportional to particle size. Using wind tunnel measurements (e.g., Greeley and Iversen, 1985), they obtained the following equation with fitting parameters:

(2) u ft 0 = A ( ρ p g D p + γ / D p ) ρ a - 0.5 ,

where g=9.81 m s−2 is the gravitational acceleration, ρp= 2650 kg m−3 is the typical soil particle density, ρa is in kg m−3, and A=0.0123 and γ=1.65×10-4 kg s−2 are empirical constants accounting for the aerodynamic forces and interparticle forces, respectively. Assuming an air density ρa=1.225 kg m−3, Eq. (2) will yield the smallest u∗ft0 of 0.21 m s−1 at Dp=80µm. For larger sizes, the particles are heavier to lift; for smaller sizes, the particles are more strongly bound by interparticle forces.

An alternative parameterization for u∗ft0 is the Iversen and White (1982) scheme (hereafter I&W82). They derived a similar solution as S&L00 but further considered the effects of soil particle size to the airflows, characterized by the particle Reynolds number Rep. I&W82 calculates uft0=uft0ρp,ρa,g,Dp,Rep in a similar form to S&L00 (see the detailed solution in I&W82 or Kok et al., 2012), with

(3) Re p = u ft 0 D p / ν ,

where ν is the kinematic viscosity of air. Since u∗ft0 is a function of Rep, which itself is a function of u∗ft0, u∗ft0 is an implicit function of itself, and the calculation needs to be iterated a few times given ρa and Dp in calculation.

Regardless of whether S&L00 or I&W82 is used, different models make different assumptions for soil particle sizes Dp, including a globally constant value (e.g., Zender et al., 2003a), a function of soil texture (e.g., Menut et al., 2013), or other forms. For instance, the Community Land Model (CLM), the land component of CESM, uses a global optimal soil diameter of Dp=75µm for the threshold calculation following Zender et al. (2003a) (Oleson et al., 2013; also see the latest version of CLM5.0 technical note on, last access: 3 October 2022). Thus u∗ft0 becomes solely a function of ρa in CESM (e.g., uft0=0.21ms-1 for ρa=1.225 kg m−3 at Dp=75µm).

Most models parameterize the effect of soil moisture fm on u∗ft following Fécan et al. (1999):



where w(long,lat,t) is the gravimetric soil moisture (kg kg−1) in the shallowest soil layer (see Sect. S1 and Fig. S1 for the relation between volumetric and gravimetric moisture); wt(long,lat) is the threshold gravimetric water content above which u∗ft increases; fclay(long,lat) is the fraction of clay content in the topmost layer of soil between zero and one; %clay=100fclay is the corresponding clay percentage; and a, a tunable constant usually of order 1, was introduced by Zender et al. (2003a) to account for the mismatch in the small scales for which Fécan et al. (1999) obtained their parameterization and the large scales on which it is used in climate models (e.g., Zender et al., 2003a; Mokhtari et al., 2012; Kok et al., 2014b). wt increases with soil clay content as water adsorbs onto clay such that more moisture is needed to enhance u∗ft.

Another essential dust emission threshold for this study is the dynamic or impact threshold u∗it, which is the lowest wind speed or stress to maintain saltation (Kok et al., 2012; Comola et al., 2019b):

(5) u it = B it u ft 0 ,

where Bit=0.82 is approximately constant with soil properties and particle size (Bagnold, 1937; Kok et al., 2012). Eqs. (2)–(5) imply that uftuft0>uit and that u∗ft and u∗it have different spatiotemporal variability. Also, the difference between u∗ft and u∗it could be much larger in nonarid regions because fm is much larger than one. In this study, we propose that dust emission models should use u∗it instead of u∗ft for dust emission equations (e.g., Eqs. 10 and 13), which will cause substantial changes in the simulated spatiotemporal variability of dust emission (see Sect. 4.1). This is needed to allow dust emission when the u is intermediate between u∗it and u∗ft, which is especially common in marginal dust source regions. Additionally, this is more physically correct as the dust emission threshold is the minimum friction velocity at which the saltation and dust emission fluxes are nonzero, which is true at u∗it but not true at u∗ft (Martin and Kok, 2018; Comola et al., 2019b; Pähtz et al., 2020).

2.2 Parameterization of drag partition effects

Apart from the dust emission threshold, another essential parameter for determining the dust emission flux is the wind drag partition effect, Feff, due to the existence of land-surface roughness elements covering the desert surfaces (Raupach, 1992; Marticorena and Bergametti, 1995). It is crucial to account for this effect for accurately simulating the magnitude and spatial pattern of dust emissions. Many past modeling studies treated this effect as increasing the dust emission threshold u∗ft (e.g., Raupach, 1992; Marticorena and Bergametti, 1995; Darmenova et al., 2009; Menut et al., 2013), such that the relation is expressed as (Raupach et al., 1993; Marticorena and Bergametti, 1995; Marticorena et al., 1997, 2006; Foroutan et al., 2017; Webb et al., 2020)

(6a) u ft = u ft 0 ( D p , ρ a ) f m ( w ) / F eff ,

where Feff<1 when roughness elements are present, such that roughness elements increase u∗ft and decrease the dust emission. However, this approach is physically incorrect because roughness elements reduce the wind stress exerted on the bare soil and do not increase the forces resisting particle lifting that determine u∗ft (Kok et al., 2014a; Webb et al., 2020). As a consequence, Webb et al. (2020) showed that dust models using Eq. (6a) will overestimate the dust emission flux compared to those using Eq. (6b). A correct emission modeling approach should instead combine Eq. (1) for u∗ft with the effect of drag partition Feff to u:

(6b) u s = u F eff ,

where u∗s is called the soil surface friction velocity (Webb et al., 2020). Dust emissions should thus be a function of u∗s instead of u.

There are different schools of drag partition schemes. A major school of drag partition parameterization originated from Arya (1975) and later Marticorena and Bergametti (1995) (hereafter M&B95), who primarily used the roughness length z0 to quantify roughness. Because of the large differences in the length scales between mountains/orography, rocks, and plants, as well as down to soil particles, Menut et al. (2013) distinguished three distinct roughness lengths describing different sizes of roughness. First, the aerodynamic momentum roughness length z0m mainly represents roughness due to large-scale orography, forests, and/or urbanization (with sizes of 10–103 m) with values ranging from  1 cm to 1 m (Menut et al., 2013). Second, the aeolian roughness length z0a quantifies the roughness due to smaller elements such as rocks and vegetation, with a typical order of magnitude of 10−3 to 10 cm (Prigent et al., 2005; Prigent et al., 2012). z0a is the relevant roughness length that informs the partition of the wind stress when considering the near-surface ( 1 m) flows in which saltation occurs. Third, the smooth roughness length z0s quantifies the roughness of a bed of fine soil particles in the absence of roughness elements. z0s characterizes the roughness of mobile, erodible soil particles over an exposed surface. z0s is directly related to the particle diameter Dp by (Nikuradse, 1950; Sherman, 1992; Pierre et al., 2014b)

(7) z 0 s = 2 D p / 30 .

M&B95 proposed their drag partition scheme by arguing that behind a roughness element (obstacle), an internal boundary layer (IBL) grows, and the wind within the IBL follows the log law of the wall as a function of u∗s and the local roughness length z0s. They then pointed out that without the obstacle, the planetary-boundary-layer (PBL) wind profile would follow the log law as a function of u and z0a. By arguing that the two wind speeds must be equal at the IBL height δ, they derived Feff as a function of z0a and z0s:

(8) F eff = 1 - ln z 0 a z 0 s ln δ z 0 s .

Later studies improved this equation based on more observations for calibrating several parameters (MacKinnon et al., 2004; King et al., 2005; Darmenova et al., 2009; see Eq. 15 in Sect. 3.2). Historically this scheme has been employed by Marticorena and others to represent the roughness due to rocks (e.g., Marticorena et al., 1997; Darmenova et al., 2009; Menut et al., 2013).

Another major school of drag partition parameterization originated from Raupach (1992) and Raupach et al. (1993) (hereafter R93), which primarily used the roughness density λ to quantify roughness. λ is defined as the total frontal area of roughness elements divided by the area of land A:λ=nhb/A, where h and b are the obstacle height and width, and n is the number of obstacles within the area. Knowing the geometric and aerodynamic properties of the roughness elements, R93 showed that the drag force of the exposed area τS is related to the total drag force τ, given λ, the roughness-element basal area-to-frontal area ratio σ, and the ratio of the roughness element-to-surface drag coefficient β:

(9a) τ S τ = 1 ( 1 - m σ λ ) ( 1 + m β λ ) ,

where m is a geometric parameter to account for the spatial variability of τS on the erodible surface. Raupach then applied this ratio to the dust emission threshold (per Eq. 6a).

Many later studies used the R93 parameterization for plants (specifically shrubs) with prescribed σ, m, and β (Darmenova et al., 2009; Xi and Sokolik, 2015). λ, however, is related to the abundance of obstacles and is thus spatially variable, and thus far there are no globally gridded datasets of λ available. Most studies thus related grid-scale λ to other grid-scale properties; for instance, Shao et al. (1996) linked λ to the vegetation cover fraction fv using in situ observations:

(9b) λ = c λ ln 1 - f v ,

where cλ is a proportionality constant. Gridded λ could thus be obtained from gridded satellite retrievals of vegetation cover (Gutman and Ignatov, 1998; Wu et al., 2016; Foroutan et al., 2017) or parameterized as a function of other gridded land-surface variables such as the leaf area index (LAI) (e.g., Klose et al., 2021). Later studies have attempted to improve Raupach's parameterization, and newer schemes relating Feff and λ have emerged (e.g., Okin, 2008).

Many previous modeling studies have not accounted for the drag partition effects of both rocks and vegetation on dust emissions (e.g., Ginoux et al., 2001; Tegen et al., 2002; Zender et al., 2003a; Kok et al., 2014b). Many past studies either accounted for only the drag partitioning by rocks (e.g., Marticorena et al., 2006; Menut et al., 2013) or by vegetation (e.g., Shao et al., 2011b; Wu et al., 2016), mainly because it is very challenging to use proxies of both rocks and vegetation in either the M&B95 or R93 scheme. For instance, R93 was historically mostly used for modeling vegetation but not rock drag partitioning because there was no dataset of the λ due to rocks. Similarly, vegetation roughness is historically mostly represented by λ rather than z0a, so there is no globally gridded z0a observations for vegetation that can be fed into the M&B95 scheme. Some modeling studies (e.g., Klose et al., 2021) generated globally gridded vegetation z0a by relating plant λ with z0a (e.g., Minvielle et al., 2003; Shao and Yang, 2005; Marticorena et al., 2006; Foroutan et al., 2017; Klose et al., 2021), but these studies all found slightly different relations between λ with z0a, and often the in situ obstacle height h is required by the relation. It is thus very challenging to model vegetation drag partitioning using M&B95 by converting λ to z0a when globally gridded h (short vegetation height but not canopy height) is mostly unknown in GCMs or possesses strong subgrid variability. A more recent approach quantifies surface roughness by detecting the shadow (sheltered area) behind a roughness element using satellite-derived albedo (Chappell and Webb, 2016). This approach could potentially capture both rock and vegetation roughness and was also employed by later dust modeling studies (e.g., LeGrand et al., 2023). To our knowledge, there are a few studies that attempted to represent both rock and vegetation roughness in one drag partition scheme (e.g., Darmenova et al., 2009; Foroutan et al., 2017; Klose et al., 2021), but all were affected by important limitations (see Sect. S2 for a discussion on their approaches). In Sect. 3.2, we will propose a novel approach that incorporates roughnesses of both rocks and plants and equally respects the z0a and λ from both schools of drag partition parameterizations, quantifying the drag partitions of rocks and plants into one hybrid drag partition factor Feff.

2.3 Parameterization of dust emission flux

There are multiple available dust emission equations (e.g., Gillette and Passi, 1988; Shao et al., 1996; Alfaro and Gomes, 2001; Ginoux et al., 2001; Tegen et al., 2002; Zender et al., 2003a; Shao, 2004; Kok et al., 2014b; Evans et al., 2016; and more) implemented in GCMs and ESMs to calculate dust emission fluxes. For example, the Zender et al. (2003a) scheme (hereafter Z03) is based on the Marticorena and Bergametti (1995) dust emission equation and is a popular dust emission scheme adopted by many GCMs (e.g., Miller et al., 2004; Oleson et al., 2013; Meng et al., 2021). Z03 calculates dust emission as follows:

(10) F d = STC MB φ f bare ρ a g u s 3 1 - u t 2 u s 2 1 + u t u s for u s > u t ,

where Fd is the emission flux (in kg m2 s−1); CMB is a proportionality constant for bridging the gap between local-scale and large-scale dust fluxes; φ=1013.4fclay-6 is the sandblasting efficiency, the vertical dust emission flux produced per unit of horizontal saltation flux as a function of soil clay fraction fclay; u∗t is the dust emission threshold (in m s−1; Z03 used u∗ft as u∗t); and fbare characterizes the fraction of land not covered by vegetation. S(long,lat) is an empirical “source function” (Ginoux et al., 2001; Zender et al., 2003a, b; Koven and Fung, 2008) to characterize soil erodibility and thus preferential source regions where fluvial sediment accumulates and scale down emission flux out of desert regions. For fbare, Mahowald et al. (2006) used a simple parameterization in which fbare is a pure function of LAI (neglecting the effects of other objects such as snow, rocks, and buildings):


While Mahowald et al. (2006) took LAIthr=0.3, we take LAIthr= 1 in this study instead because (1) observations show that there could be dust emitted from semiarid regions with LAI > 0.3 (Okin, 2008); and (2) Mahowald et al. (2006) did not account for wind drag partitioning due to plants, and thus by setting a small LAIthr, emission (Fdfbare) drops more rapidly with LAI such that the drag partition effect is also incorporated in the fbare term. However, since we are considering Feff in this study, we can set a more realistic LAIthr value such that fbare becomes less sensitive to LAI.

In this study, we use the Kok et al. (2014b) dust emission equation (hereafter K14) which is increasingly adopted by more GCMs (e.g., Evan et al., 2015; Ito and Kok, 2017; Mailler et al., 2017; Tai et al., 2021; Li et al., 2021; Klose et al., 2021; Li et al., 2022). One key advance of K14 over Z03 is that K14 eliminated the need to use an empirical, time-invariant source function S to tune the spatial variability of dust emissions. K14 proposed that a dynamical and time-varying soil erodibility (named Cd in K14) can be physically parameterized using the standardized fluid threshold ust=uftρa/ρa0, which is u∗ft scaled to the standard air density of ρa0=1.22 kg m−3:

(12) C d = C d 0 exp - C e u st - u st 0 u st 0 ,

where Cd(long,lat,t) is the time-varying dust emission coefficient or soil erodibility coefficient, Cd0=(4.4±0.5)×10-5, Ce=2.0±0.3, and ust0=0.16 m s−1. Furthermore, K14 derived a new dust emission equation for Fd (kg m−2 s−1):

(13a) F d = C tune C d f bare f clay ρ a u s 2 - u t 2 u st u s u t κ for u s > u t

(13b) κ = C κ u st - u st 0 u st 0 ,

where Cκ=2.7±1.0, Ctune=0.05 is the proportionality constant, fbare is modeled by Eq. (11), us=u in K14, and u∗t is again the emission threshold (K14 assumed for simplicity that ut=uit=uft). κ is the fragmentation exponent which quantifies the sensitivity of Fd to u∗s. Here we limit the value of κ to 3 in order to prevent excessive sensitivity of the model to wind speeds, which can be problematic around topography. From Eq. (12), Cd increases exponentially with u∗st, and thus K14 dust emission is very sensitive to u∗ft. K14 showed improvements compared with Z03 when evaluated against ground-based DAOD measurements (Kok et al., 2014a, b; Li et al., 2022).

2.4 Input required by dust emission schemes

Calculating the above dust emissions and aeolian processes requires meteorological and land-surface variables as inputs. We employ the required input data from the Modern-Era Retrospective Analysis for Research and Applications version 2 (MERRA-2) (Gelaro et al., 2017). MERRA-2 is a reanalysis dataset provided by NASA's Global Modeling and Assimilation Office (GMAO). MERRA-2 has a native resolution of 0.5× 0.625 and hourly data assimilation. All MERRA-2 modeled fields and other input variables in this study are listed in Table 1. In this study, we code the dust emission scheme with all new aeolian processes in the statistical programming language R (v4.2.1) as an offline (outputs do not feedback onto input forcings), standalone sandbox model. In this study, we use the standalone model to read in all input atmospheric and land surface forcings for the year 2006 and employ equations in Sects. 2 and 3 to compute 2006 dust emissions as outputs and results for Sects. 3 and 4.

Table 1Input meteorological and land-surface variables employed for running the standalone dust emission model in this study.

Download Print Version | Download XLSX

3 Physics-based parameterization of dust emission threshold

In this section, we propose additions and improvements to several parameterizations of dust emission physics, which include (1) deriving a more realistic soil median diameter map and including it in the u∗ft calculation, (2) proposing a new hybrid approach to incorporate the drag partition parameterizations of both rocks and vegetation, and (3) implementing a parameterization of the effects of turbulence on the intermittency of dust emissions. We will use the improved model from this section to compute hourly dust emissions in Sect. 4, driven by meteorological and land-surface fields.

3.1 Improving the description of soil particle size parameter

The PSD of the soil bed is a critical factor to determine the dust emission threshold. In this section, we focus on deriving a new global soil median diameter (a good proxy for the soil PSD) (Martin and Kok, 2019) as a parameter for computing the dust emission thresholds. Section 5 discusses the caveats and limitations of this approach.

3.1.1 Motivation and literature compilation of soil particle size distribution

As discussed in Sect. 2.1, Martin and Kok (2019) argued that u∗ft of a mixed soil should be determined by the median diameter Dp of the soil PSD. Thus, we ideally need a global gridded map of Dp to calculate u∗it and u∗ft over the globe. However, there are only very limited in situ measurements of soil PSDs (e.g., see Table S1) that are insufficient to compile a global Dp map. Meanwhile, extensive studies have compiled global maps of many other soil properties, such as soil texture, soil bulk density, pH value, soil organic carbon (SOC), and cation exchange capacity (e.g., FAO/IIASA/ISRIC/ISS-CAS/JRC, 2012; Shangguan et al., 2014; Hengl et al., 2017; Dai et al., 2019). Therefore, to determine and predict Dp, we use a compilation of literature measurements to explore and construct relationships between Dp with other soil properties such as the clay and silt fractions.

However, many past laboratory studies used the wet sedimentation or wet sieving technique to measure the texture of the soil samples. Wet sieving effectively breaks down soil microaggregates into disaggregated particles and can dissolve soluble minerals (Chatenet et al., 1996), thereby disturbing the estimations of the in situ soil median particle sizes. In contrast, dry sieving causes a minimal disruption to soil microaggregates, and thus Chatenet et al. (1996) argued that the dry-sieved soil PSDs are more representative of the in situ, aggregated soil PSDs. Although the soil texture is a disaggregated soil property, Dp might depend on soil texture fa and other soil properties because the strength of interparticle forces is contingent upon soil texture (the clay and silt content), which governs the extent of soil aggregation. Here, we use measurements from past laboratory studies (see Table S1), which contain site-scale, dry-sieved soil PSDs, wet-sieved soil texture fa, and other soil properties to investigate their statistical relations and infer a new global distribution of Dp. All studies listed in Table S1 have dry-sieved soil PSD measurements and the wet-sieved sand, silt, and clay fractions. Many studies have recorded soil organic carbon (SOC; %) and other properties such as calcite (CaCO3; %), pH value, and bulk density (g cm−3). Figure 1a shows the locations of the measurements of the employed soil studies, and the colors show the aridity where the sites are located. Some studies obtained measurements over a relatively large spatial domain, and we plot only one symbol at the domain centroid representing multiple measurements. Many studies reported PSD measurements extending to diameters in excess of 6000 µm, but we used only PSD measurements in the diameter range of 0 and 2000 µm that is relevant to dust emission (Zender et al., 2003a). For each dry-soil PSD measurement, we obtain the aggregated Dp by calculating the 50th percentile of the dry-soil PSD.

3.1.2 Deriving a global soil median diameter map

We classify the datasets into arid and nonarid groups, since we are primarily interested in Dp over desert regions (although we also display the soil behaviors over nonarid regions). We follow past studies (Mahowald et al., 2006, 2010; Kok et al., 2014b) which defined arid (or dust emission) regions using the criterion of LAI smaller than a threshold LAIthr, which we take to be 1 (see Sect. 2.3). Section 3.2.2 also describes the MERRA-2 LAI we used in this study to identify the world's arid regions.

After dividing the data into median dry diameters for arid and nonarid soils, we examine the statistical relationships between Dp and the soil properties (see Figs. 1b and S2). Figure 1b shows a scatterplot of Dp versus the sum of the soil component that produces substantial cohesion, namely the silt and clay fractions (fsilt+clay=fsilt+fclay). The data exhibit distinctly different trends for nonarid versus arid soils: for nonarid soils, Dp increases from 100 µm to greater than 1000 µm with fsilt+clay (regression p value =7.3×10-6), likely due to increasing cohesion with increasing clay and silt content. In contrast, Dp for arid soils shows a small and statistically insignificant increasing trend with fsilt+clay (p value = 0.77) with a smaller Dp variability (50–250 µm). This flat trend indicates that fsilt+clay does not effectively explain the median diameter of aggregated soil particles in arid regions. We examined the relationships of Dp with the individual fractions of sand, silt, and clay, as well as with other soil properties including SOC, pH, and CaCO3 (Fig. S2), but these relationships are not statistically significant. We obtain a surprisingly simple finding from the available measurements that there is limited variability in the aggregated Dp over the arid regions across different soil textures. We thus use a constant Dp0 as an approximation for arid regions. From Fig. 1b, we summarize the relationship between Dp and fsilt+clay as

(14) D p = Ψ 0 + Ψ 1 f silt + clay D p 0 , for LAI > LAI thr for LAI LAI thr ,

where Ψ0=7.8±3.7µm, Ψ1=124±36µm, Dp0=127±47µm, and LAIthr=1 as specified in Eq. (11). This empirical formula suggests that some models' assumptions of the relationship between Dp and soil texture were inaccurate (e.g., Table 2 of Laurent et al. (2008) assumed Dp decreases with fsilt+clay), and this result could substantially simplify model parameterizations. Additionally, our diameter of 127 µm over arid regions is larger than Z03's assumption of a globally constant optimal diameter of 75 µm. This translates to a modest increase of u∗ft0 from 0.204 to 0.216 m s−1 (given ρa=1.225 kg m−3), which slightly decreases global dust emissions by 18 % (see Sect. 4.1). The uncertainty in Dp0=127±47µm translates to an uncertainty of u∗ft0 between 0.204 to 0.234 m s−1.

Figure 1Constructing a global map of the median diameter Dp of aggregated soil particles using soil particle size and texture data. (a) The locations of literature measurements, with symbols indicating the names of the studies and the color indicating the aridity of the locations; a site is classified as arid (red color) if its location has MERRA-2 LAI < 1 and otherwise nonarid (blue color). (b) Literature measurements of soil dry median diameter Dp versus silt + clay fraction (fsilt+clay). (c) The predicted soil median diameter Dp map (in µm) derived by projecting our derived Dpfsilt+clay relationship of Eq. (14) on the SoilGrids (Hengl et al., 2017) soil texture data (Fig. S3). The circles represent the locations of the sites the same as panel (a), and their colors show the measured median diameter Dp at those sites. (d) The predicted Dp using Eq. (14) versus measured Dp from past studies.

We then project our derived relation between Dp and fsilt+clay on the available soil texture and properties database. We employ global soil properties data from the SoilGrids database (Hengl et al., 2014, 2015, 2017), a global soil mapping project that used machine learning (random forest) to regress in situ measurements of soil variables (moisture, temperature, nutrients, etc.). SoilGrids provides global maps of soil texture and other soil properties with a horizontal resolution of 250 m and eight soil depths down to 200 cm (Hengl et al., 2017). We use SoilGrids instead of other available soil databases as it shows better performance against observed soil profiles than other soil databases (Dai et al., 2019). Figure S3 shows the SoilGrids relative fractions of sand, silt, and clay with a 0.1× 0.1 horizontal resolution for the topmost soil layer. Figure 1c shows our global 0.1× 0.1 soil median diameter Dp map. Following Eq. (14), the arid and semiarid regions are set to have a Dp0 of 127 µm, whereas for nonarid regions, Dp increases with fsilt+clay. Our derived Dp values are largely consistent with the site Dp measurements from past studies (overlaid points), showing a similar spatial distribution, with a fit-line slope of 0.98 (p value = 0.007) and an R2 of 81 % (Fig. 1d). Note that, since the predictions for arid regions (red points) are a constant without variability, the agreement between the predictions and observations is essentially dominated by the linear Dpfsilt+clay relation over nonarid regions (blue points). Figure 1d shows that Eq. (14) gives satisfactory agreement in predicting global Dp, but dust emission modeling will depend exclusively on the predicted Dp over arid regions. We anticipate that as more measurements emerge in the future, more statistical or machine learning modeling approaches can more robustly decipher the intricate relationships between Dp and various soil properties over arid regions.

Since nonarid regions of LAI > 1 will generate zero emissions (Eq. 11), we simplify Eq. (14) and Fig. 1c by imposing a globally constant Dp0=127µm.

3.2 A wind drag partition scheme for decreasing wind stress and erosion

We now present a methodology to account for the wind drag partition effect due to nonerodible roughness elements including vegetation and rocks that protect the bare soil by absorbing part of the surface wind stress. We calculate the rock drag partition feff,r using z0a since global z0a observations are available, and we calculate the vegetation drag partition feff,v using vegetation cover which is a proxy of λ (e.g., Shao et al., 1996; Okin, 2008), since gridded plant cover is often parameterized in GCMs (e.g., Wu et al., 2016; Foroutan et al., 2017; Meier et al., 2022). Here we use two separate drag partition schemes (Marticorena and Bergametti, 1995; Okin, 2008) to quantify the roughness effect of rocks (Sect. 3.2.1) and vegetation (Sect. 3.2.2), respectively. Then, we propose a unifying approach to combine the two effects into a hybrid factor Feff (Sect. 3.2.3).

3.2.1 Drag partition due to roughness of rocks

In this study, we use the aeolian roughness length z0a to quantify the drag partition effect due to rocks. Whereas the smooth z0s and the aerodynamic momentum z0m can be derived from pre-existing datasets, it is more challenging to quantify the aeolian z0a. Existing efforts employed satellite and field measurements to quantify the roughness over deserts (e.g., Greeley et al., 1997; Roujean et al., 1997; Marticorena et al., 2004; Laurent et al., 2005; Prigent et al., 2005; Marticorena et al., 2006; Prigent et al., 2012). For instance, Marticorena et al. (1997) and Callot et al. (2000) developed a 1× 1 z0a map over Africa and the Middle East by combining topographic data, geological information, aerial pictures, and in situ observations. Prigent et al. (2005) and Prigent et al. (2012) further used radar measurements to yield global maps of backscatter coefficient, which is a measure of surface roughness because rougher surfaces generally scatter more radar signals to different directions and reduce the backscattering. Comparisons between satellite backscattering signals and field measurements of z0a yielded an empirical formula for extrapolating a global dataset of backscattering signal to global z0a. We use here the global aeolian z0a dataset from Prigent et al. (2005) (hereafter Pr05), which contains the climatological monthly mean z0a (12 monthly values per grid) derived from the backscatter coefficient observed by the scatterometer at 5.3 GHz on board the European Remote Sensing (ERS) satellite. Since satellite z0a measurements could quantify the roughness of both rocks and vegetation, we take the minimum value out of the 12 months for all grids to obtain a static aeolian z0a map to eliminate as much as possible the vegetation effect on the inferred roughness. Furthermore, we apply this map over arid regions only (LAI < 1), where the backscatter signal is mostly generated by rocks with little contribution from vegetation roughness. The resulting 2-D map of z0a (in centimeters) thus mostly represents time-invariant rock roughness and is plotted in Fig. 2a.

Marticorena and Bergametti (1995) derived a parameterization to quantify the drag partition effect using both z0a and z0s. They assumed that this equation is valid for roughness elements that are not too closely spaced (small wake), i.e., z0a<1 cm (Darmenova et al., 2009). Here we use their semiempirical equation to quantify the drag partition due to rocks, feff,r (also see Eq. 8):

(15) f eff , r = 1 - ln z 0 a z 0 s ln b 1 X z 0 s b 2 ,

where X is the distance downstream the point of discontinuity in roughness, a length parameter that scales with the IBL height δ behind the obstacle in Eq. (8), i.e., δz0s=b1Xz0sb2 following Marticorena and Bergametti (1995), and b1=0.7 and b2=0.8 are empirical constants (King et al., 2005; Darmenova et al., 2009). X should be a function of land type and implicitly space and time, but thus far most dust modeling studies have used a global constant for X (e.g., Darmenova et al. (2009) used a globally constant X=0.1 m). We use a globally constant X=10 m in this study, which is different from what the past studies suggested, because the scale of the rocks and plants we focus on in deserts is larger and is of the order of 100–101 m. Some studies considered even larger roughness and used X∼122 m for vegetated deserts (MacKinnon et al., 2004). We then obtain z0s from our derived global dataset of Dp in Sect. 3.1 using Eq. (7). When nonerodible roughness elements are abundant over a surface, z0az0s and feff,r1, causing the sheltering of the bare soil from the wind; when there are few roughness elements, z0a is small and close to z0s, and thus feff,r approaches 1. In Fig. 2a, the red areas with very small z0a are the most susceptible regions for dust emission. Figure 2a also shows that most arid and semiarid regions have z0a<0.2 cm, such that Eq. (15) can follow the criterion (z0a<1 cm) in Darmenova et al. (2009) well. Figure 2b shows the global feff,r over arid regions, which is dominated by the spatial pattern of z0a in Fig. 2a given feff,r is governed purely by z0a.

Figure 2Global roughness length and rock drag partition factor maps at a horizontal resolution of 0.25× 0.25. (a) Global mesoscale aeolian roughness length z0a (in centimeters) derived by Prigent et al. (2005). (b) Global static rock drag partition factor feff,r derived by Eq. (15) following Marticorena and Bergametti (1995), derived over arid and semiarid regions defined as MERRA-2 LAI < 1 in this study. The color schemes are set such that the most erodible regions appear red.

3.2.2 Drag partition due to roughness of vegetation

Unlike the very static and slowly evolving rock roughness, vegetation changes temporally. To include the effect of these dynamic vegetation changes on the drag partition, we follow the approach of Okin (2008) (hereafter O08), which uses unvegetated gap size (the distance between neighboring plants) to characterize the variability of the reduced wind stress. O08 argued that his scheme represents an advancement over the classical R93 scheme, since R93 uses the roughness density (or lateral cover) λ, which only quantifies how much roughness is on a surface but not how that roughness is spatially distributed. O08 pointed out that, given the same λ, roughness elements divided into small blocks spread over the soil surface would be more effective than elements stacked up like a telephone pole in partitioning wind stress (see Fig. 3 in Okin, 2008). Okin argued that since Raupach's model neglects the spatial variability of λ, the resulting simulated emission flux using the R93 scheme in Okin's paper decreased rapidly with increasing λ and unrealistically reached zero at relatively low λ. To partially compensate for this error, R93 introduced a tuning parameter m (Eq. 9a), serving to reduce the effective λ and thereby reducing the rapid decrease in dust flux. However, m is a tuning parameter not derived from first principles, and it is not clear how m changes over different surface conditions. Therefore, we use the O08 model here to better characterize the spatial variability of wind stress and the resulting dust emissions.

Here we describe the O08 scheme and adapt it for use in LSMs and GCMs. O08 assumes u drops significantly when encountering a roughness element (plant) and gradually recovers at the lee (downwind region) of the plant as a function of distance x, following

(16a) u s ( x / h ) = u f 0 + ( 1 - f 0 ) 1 - e - x / h c ,

where x/h is the dimensionless downwind distance from an obstacle normalized by vegetation height h (m), f0=usu|x=0 is the friction velocity ratio immediately behind the obstacle, and c is the dimensionless e-folding distance (normalized by h) over which u∗s locally recovers to u. In this formulation, the local drag partition factor due to vegetation as a function of distance x is

(16b) f local x h = u s u = f 0 + ( 1 - f 0 ) ( 1 - e - x / h c ) .

Note that in the limit of x/h,usu. O08 used measurements from Bradley and Mulhearn (1983) and fitted f0=0.32 and c=4.8 (i.e., the e-folding distance of u∗s recovery to u is 4.8 times the plant height h) for semiarid regions.

In order to use Eq. (16b) to obtain the drag partition feff,v relevant to a regionally vegetated area that is more applicable to GCMs, one needs to calculate an integral for the averaged and aggregated effect of drag partitioning feff,v (see Eq. 20a) instead of a locally varying flocal (Eq. 16b). Therefore, Okin employed a probability distribution function as a function of distance x/h to indicate the importance (or weight) of flocal at any x/h to the averaging of feff,v (McGlynn and Okin, 2006; Okin, 2008). The PDF is an exponential decay such that the weight of flocal decreases with distance x/h, so flocal at the immediate lee of the obstacle (which is smaller and close to f0) has more weight than the flocal farther away (which is larger and tends to 1). From McGlynn and Okin (2006), the PDF is a function of normalized distance x/h:


where L (m) is the mean gap length between obstacles (plants), which is conceptually related to fv; and K is the normalized gap length, which is the gap length L scaled by the plant height h. Physically, Pd is the probability that there is not another obstacle present within a downwind distance x/h. This exponential decay implies that the farther away from a plant (larger x/h), the higher the likelihood that there is another plant present within the downwind distance x/h, with the normalized gap length K quantifying the e-folding distance of the probability. This PDF governs the spatial domain over which u∗s recovers.

For O08, the mean gap length between obstacles K is the only required input for calculating the drag partition, since f0 and c are assumed to be invariant to surface conditions and desert biome. K can be expressed as a function of fv using some simple assumptions. First, O08 argued that the vegetation cover fraction is simply fvWL+W, where L is the mean gap length and W is the mean width of the plants within that vegetated area. Rearranging gives

(18a) L = W ( 1 f v - 1 ) .

Then we assume plants in arid regions (e.g., shrubs) are approximately hemispheres with radius R. Then, the plant height h=R and width W=2R are related by W=2h, which can be substituted into Eq. (18a) to yield

(18b) K L h = 2 1 f v - 1 ,

and thus we related K to fv. fv could be measured at the local level, and thus O08 was frequently applied in field studies (e.g., Li et al., 2013; Pierre et al., 2014a). However, what is novel in our study is that we are the first to propose the implementation of O08 into LSMs, because Eq. (18b) shows us that O08 could be formulated as a function of fv, which is a grid-level parameter. Here we propose to follow the Mahowald et al. (2006) assumption in Eq. (11) and approximate vegetation cover fraction as fv=1-fbare=LAI/LAIthr. Equation (18b) becomes

(18c) K L h = 2 1 f v - 1 = 2 ( 1 LAI / LAI thr - 1 ) ,

where we assume LAIthr=1 (in Eq. 11). The assumption of fv∼LAI is valid if we reasonably assume that leaf areas over arid regions overlap relatively little with each other. We note that by using LAI to quantify fv in Eq. (18c), we are only accounting for the vegetation drag partitioning due to green (photosynthetic) vegetation and miss that due to brown (nonphotosynthetic) vegetation. In the future, it is warranted that Eq. (18b) includes other proxies of brown vegetation drag partitioning, such as the vegetation cover quantified by Guerschman et al. (2015), which was adopted by later dust modeling studies such as Klose et al. (2021) and Huang and Foroutan (2022).

To estimate the reduced emission flux, O08 uses an integration approach without quantifying feff,v. O08 calculates the reduced dust emission flux Fred (kg m−2 s−1) by locally integrating the emission Fd following the spatially varying u∗s over the normalized distance x/h:

(19) F red = x / h = 0 P d x h F d u s x h d x h ,

where Fd (kg m−2 s−1) is the local emission as a function of u∗s which is itself a function of x/h. In the integration, Fd needs to be weighted by Pd (which means Fd at large x/h has proportionally less importance) because as x/h increases, the likelihood of the presence of another obstacle gets larger and larger, which will hinder the recovery of u∗s to u. Integrating the emission flux Fd from zero to infinity gives a reduced emission flux Fred, which will be smaller than the emission flux without roughness elements, defined as Fbare=x/h=0PdxhFdudxh=Fdu, in which Fd is a constant in space since us=u is a constant without obstacles.

However, since we need to also combine the vegetation drag partition with the rock partition effects, we need to quantify feff,v in order to form a hybrid drag partition factor for LSMs. Instead of directly implementing Eq. (19a) into LSMs, we require an alternative approach of quantifying feff,v such that Fdfeff,vu=Fred. Quantifying feff,v for O08 can be useful for comparisons against feff,v from other schemes such as R93 and Klose et al. (2021). In addition, quantifying feff,v for O08 makes it possible to generate a high-resolution, diagnostic feff,v dataset for mechanistic models with different resolutions as a model input.

An approach of evaluating feff,v from O08 was proposed by Pierre et al. (2014a). They obtained the expected value of the shear stress ratio us/u (SSR in Okin, 2008) between obstacles by evaluating the integral of us/u weighted by Pd, which represents the averaged flocal (in Eq. 16b) across the vegetated area and perfectly fits our purposes for implementing feff,v into LSMs:

(20a) f eff , v = x / h = 0 P d x h u s x h u d x h = x / h = 0 P d x h f local x h d x h .

Substituting Eq. (17) for Pd into Eq. (20a) and analytically evaluating the integral gives a simple algebraic equation for feff,v (Pierre et al., 2014a), representing the aggregated vegetation drag partition effect at the grid level:

(20b) f eff , v = K + f 0 c K + c .

This elegant formula conveys a clear physical intuition: if the obstacle does not effectively dissipate momentum (f0→1), feff,v1; if land is densely covered by vegetation (gap length K→0), feff,vf0(=0.32), the shear stress ratio at the immediate lee of the obstacle. An advantage of this approach is that it can be easily adopted by gridded models since modelers only need to code an algebraic equation instead of an integral.

We calculate 0.5× 0.625 global hourly feff,v data for Okin's model, using Eqs. (18c) and (20b) with hourly MERRA-2 LAI. We note that MERRA-2 LAI is based on the Advanced Very High Resolution Radiometer (AVHRR) observations (Reichle et al., 2017). Figure 3a shows the annually averaged MERRA-2 LAI for the year 2006 over arid regions with LAI < 1 (seasonal LAI maps are also shown in Fig. S4), and Fig. 3b shows the corresponding mean feff,v for areas where LAI < 1. The LAI plot shows the most erodible regions on Earth.

Figure 3Vegetation drag partition factor feff,v derived from the Okin (2008) and Pierre et al. (2014) drag partition model for the year 2006 on a 0.5× 0.625 grid. (a) Annual mean MERRA-2 LAI, with color bar saturated at a value of 1. (b) Annually averaged feff,v derived using the Okin (2008) and Pierre et al. (2014) drag partition model. White areas indicate water body, ice/snow, or LAI > 1.

3.2.3 Combining drag partition factors of rocks and vegetation

After obtaining both the static feff,r map of rocks and the time-varying feff,v map of vegetation, we now propose a methodology to combine the two drag partition sources to capture and represent the total drag partition effect for dust emission. LSMs need a single drag partition factor capturing all roughness effects to estimate the total reduction of the surface winds. Thus, we compute a hybrid drag partition factor map Feff that can be used as input for dust modules in GCMs. To achieve this, we need to know the fractions of a grid that consists of areas dominated by rocks and areas dominated by plants, which can be obtained from several recent studies (Lawrence et al., 2016; ESA, 2017; Klein Goldewijk et al., 2017; Kobayashi et al., 2017). We obtained these data from the European Space Agency Climate Change Initiative (ESA CCI) dataset (, last access: 21 June 2022). The land cover product classifies the land cover of the whole globe into 37 categories (Li et al., 2018), with relevant land cover over arid regions such as shrub, herbaceous, sparse vegetation, cropland, grassland, and consolidated (gravels and rocks) and unconsolidated (soil) bare land. This dataset has a horizontal resolution of 300 m, making the dataset capable of counting the portion of the grid consisting of rocks and vegetation over a larger MERRA-2 0.5× 0.625 grid box (a MERRA-2 grid box consists of  35 000 grids of 300 m). This dataset gives a representation of the annually varying land covers, so the rock and vegetation area fractions we use are a function of space only within the simulation year of 2006. We describe our approach to synthesizing the ESA CCI land cover maps and drag partition datasets in the following.

We incorporate the drag partition effects by identifying two roughness regimes using the ESA CCI dataset. The first regime is the rock regime (Fig. 4a), for which we combine the consolidated (gravel and rocks) and unconsolidated (soil) bare land types (types 34–36). This regime is subject to the rock drag partition effect. The second regime is the vegetation regime (Fig. 4b), which includes different vegetation types such as shrubland and herbaceous (types 19–23, 28–29, 32), sparse vegetation (types 26–27), cropland (types 2–5), grassland (type 24), mixed vegetation (type 18), and other vegetation mosaic (types 6–7). Since O08 does not specify the differences in drag partition for different plant functional types (PFTs), here we assume all PFTs produce the same drag partition effect. The overall drag partition effect Feff for a grid is thus defined by the summation of emissions, with emission Fd,r over the rock regime with a fractional area of Ar, and emission Fd,v over the vegetation regime with another fractional area Av:

(21a) F d u F eff = A r F d , r + A v F d , v = A r F d ( u f eff , r ) + A v F d ( u f eff , v ) .

Given that dust emissions approximately scale with the cube of u∗s (Zender et al., 2003a; Kok et al., 2014b) and neglecting the effect of the dust emission threshold, Eq. (21a) can be simplified to

(21b) F eff 3 = A r f eff , r 3 + A v f eff , v 3 ,

such that Feff is simply the weighted mean of drag partition effects. The fractional areas are simply calculated by counting the total occupied area of the ESA CCI land cover corresponding to a certain regime and then dividing by the total area of the grid box. We use Eq. (21b) to obtain the spatiotemporally varying Feff(long,lat,t), given feff,v(long,lat,t) and Ar, Av, and feff,r as functions of (long,lat). We then apply the obtained Feff here to Eq. (6) to yield u∗s for the dust emission equation. We discuss in Sects. 5 and S6.2 the caveats and limitations of this hybrid drag partition scheme.

Figure 4a–b show the fractional areas of the two regimes. The rock regime (Fig. 4a) is located mostly over the Sahara, the Middle East, and the Asian deserts. The vegetation regime (Fig. 4b) is concentrated mostly over Australia, the United States, South America, and southern Africa. Figure 4c shows the resulting annually averaged Feff using Eq. (21b). The regions with the highest Feff are the Bodélé Depression, El Djouf, the Arabian Desert, and Taklamakan due to high feff,r. The Strzelecki–Sturt Stony deserts in Australia, the Kyzylkum, and Patagonia also have high Feff ( 0.7) due to high feff,v. Regions with both high rock and vegetation roughness are located in parts of the Middle East and North America with low Feff values.

Figure 4The 0.5× 0.625 hybrid drag partition factor Feff incorporated using the European Space Agency Climate Change Initiative (ESA CCI) dataset. (a–b) The fractional areas of (a) the rock regime (consolidated/unconsolidated land) and (b) the vegetation regime (shrubs, herbaceous plants, croplands, grassland, and sparse vegetation) over arid regions. (c) The hybrid drag partition factor Feff by a combination of rock drag partition feff,r and vegetation drag partition feff,v for the year 2006.

3.3 Parameterizing the dust emission intermittency

The above improvements enable a more accurate calculation of emission when wind speeds are sufficient to initiate dust emission. Next, we will improve the calculation of the resulting dust emission flux by accounting for the effects of boundary-layer turbulence on dust emission intermittency. Dust emission intermittency exists because saltation is driven by turbulent surface winds, which exhibit strong spatiotemporal fluctuations in speed and direction. Instantaneous winds can thus pass within short timescales across the emission thresholds for initiating or ceasing saltation (Martin and Kok, 2018). Consequently, saltation can be highly intermittent (Comola et al., 2019b), with pronounced variability on timescales of seconds to hours (Dupont et al., 2013). In contrast, existing dust emission parameterizations describe saltation as uniform in time and space and driven by a constant downward momentum flux within a model time step. The disconnect between the reality of intermittent dust emissions and uniform emissions in current theories is likely contributing to the poor performance of dust emission simulations (Barchyn et al., 2014; Todd et al., 2008). Comola et al. (2019b) argued that the intermittency effect is more prevalent for regions with low-intensity dust emissions when u∗s is regularly fluctuating around the threshold to turn on or shut off dust emissions. Neglecting intermittent dust emissions in current models thus likely degrades the accuracy of dust emission simulations for arid regions during low-wind periods and for marginal dust source regions such as semiarid areas (since soil cohesion increases u∗ft but does not affect u∗it), which dominate in much of the Southern Hemisphere (Ginoux et al., 2012; Ito and Kok, 2017).

Accounting for the intermittency effect on dust emission fluxes is complicated by the hysteresis of dust emission due to the existence of double thresholds for dust emission physics. The instantaneous wind at the saltation level ũs (we use the tilde to denote instantaneous quantities and take away the asterisk to denote winds at the saltation level of zsal∼0.1 m using Eq. (S4a) instead of a velocity scale) needs to exceed the fluid threshold uft (also defined at the saltation level) to initiate saltation, but it only needs to exceed a smaller impact threshold uit to sustain it (Kok et al., 2012; Martin and Kok, 2018; Comola et al., 2019b). When ũs at a moment lies between both thresholds (uit<ũs<uft), saltation is active if transport was more recently initiated (ũs>uft) and inactive if transport was more recently terminated (ũs<uit). This process is known as hysteresis (Kok, 2010; Martin and Kok, 2018; Comola et al., 2019b). As a result, if us (mean of ũs within a model time step) is between uit and uft, there will be fluctuating emission fluxes in reality, while models using a fluid threshold scheme would predict zero emission within a model time step, thereby underestimating the emissions. Meanwhile, models using an impact threshold scheme without considering turbulence will have uniform positive dust emission within the time interval. However, because in reality high-frequency winds can pass below uit and shut off dust emissions, using average us in an impact threshold scheme will overestimate dust emissions. It is thus important for GCMs to account for the effects of turbulence causing both intermittency and hysteresis of dust emission.

As GCMs have a relatively large time step and a coarse horizontal resolution (e.g.,  30 min for a 1 GCM), they are not designed to resolve turbulence and cannot capture high-frequency ( 0.1–5 min) turbulent wind speed fluctuations. As a result, models cannot directly simulate the dust emission intermittency. Therefore, accounting for intermittent dust emission requires a parameterization that links the low-frequency ( 30 min) variables of boundary-layer turbulence that are resolved in GCMs to the high-frequency intermittency dynamics. Comola et al. (2019b) formulated a parameterization (hereafter the C19 scheme) of intermittent saltation fluxes by quantifying wind fluctuations due to both shear-driven and buoyancy-driven turbulence in terms of resolved model parameters, including uit and the Monin–Obukhov length L. C19 showed that when a dust emission equation employs uit and accounts for the intermittency effect, it can successfully capture the magnitudes of small dust fluxes otherwise missed by models using uft (Fig. 3 of Comola et al., 2019b). The C19 scheme will thus moderate the temporal variability of modeled dust emissions due to diurnal wind cycles continuously crossing the thresholds. Additionally, it will also capture more lower-intensity emissions over marginal sources missed by many current models (Zhao et al., 2022).

In the C19 scheme, the dust emission flux Fd is calculated using u∗it instead of u∗ft. We update K14 (Eq. 13) with ut=uit as the threshold, giving


where Cd is still a function of u∗st, and ust=uftρa/ρa0 is the same standardized fluid threshold as in the default K14, and u∗it is computed using Eq. (5). Because u∗it <u∗ft, this modified equation allows more small dust fluxes over the marginal source regions that are otherwise missed by employing u∗ft as the threshold (see Figs. 7g–h).

Next, we account for the intermittency effect by introducing the intermittency factor η, which is the fraction of time that saltation is active in a model time step (e.g.,  30 min). η corrects the horizontal sand saltation flux, which scales with dust emission flux (Shao et al., 1993), thereby also representing the fraction of time that dust emission is active in a model time step. C19 accounts for the effect of intermittency by multiplying dust emission by η as follows:

(22b) F d , η = η F d ,

where η[0,1]. Note that C19 parameterizes η using wind information at the typical saltation height of zsal=0.1 m instead of the velocity scales. η is thus formulated as a function of the wind speeds us, uit, and uft at height zsal (see Sect. S3 Eqs. S3–S6) and the standard deviation σũs of the instantaneous ũs (Eq. 23):

(22c) η = η ( u s , σ u ̃ s , u it , u ft ) .

σũs is defined given that ũs can be described by a normal distribution (Chu et al., 1996), with its mean being the model time step mean (at 0.1 m) us and its standard deviation σũs. From Eq. (22c), η→1 when usuft (active emission for the whole time step) and η→0 when usuit (no emission for the time step). The further away us is from the thresholds uft and uit, the smaller the probability of the instantaneous ũs sweeping across the thresholds and the more dichotomous η behaves (either zero or 1). If us is very close to uft or uit, or is indeed between them (uit<ũs<uft), the frequency of crossing the threshold is determined by the magnitude of the turbulent fluctuation σũs. σũs is parameterized using the similarity theory (Panofsky et al., 1977):

(23) σ u ̃ s = u s 12 - 0.5 z i L 1 / 3 for 12 - 0.5 z i L 0 ,

where L is the Obukhov length, and zi is the modeled PBL height. Note that MERRA-2 does not provide L output, and in this study we computed L from the MERRA-2 outputs of u, ρa, sensible heat flux H, and temperature T for our simulations (see Sect. S3). In boundary-layer dynamics, turbulence is generated by mechanical shear and buoyancy (Stull, 1988). From Eq. (23), high-frequency wind fluctuations σũs increase with shear (us>0) and buoyancy (L<0). A larger σũs makes it easier for ũs to sweep across uit and shut off dust emission, leading to η<1. In a time step, if usuft+σũs, ũs will be unlikely to sweep across uit, and η will approach 1. If us is slightly larger than uft, the instantaneous ũs will be likely to sweep across uit, leading to η<1. In the hysteresis regime (uit<us<uft), η will be around 0.3–0.7, since ũs will sweep across both thresholds given σũs, leading to a reduced emission flux (meanwhile, other parameterizations predict a zero emission flux since they use uft only). When us<uit, η could also be greater than zero when σũs is large enough so that the instantaneous ũs sweeps across uit. However, C19 would not generate any emission according to Eq. (22a) (which is a technical flaw of C19; see a discussion in Sect. S6.3). We note that Eq. (23) is not the traditional Monin–Obukhov similarity theory, as the zonal fluctuation was shown to correlate poorly with z/L but relates much better with zi/L (Panofsky et al., 1977). We also note that Eq. (23) only applies to the convective PBL, but dust emission often occurs during the daytime within the convective boundary layer (Yu et al., 2021). With the complete C19 scheme in Eqs. (S3)–(S6), we can compute η to yield the dust emission with intermittency effect Fd,η as the final dust emission for the LSM. The full C19 intermittency scheme is described in Sect. S3 and also discussed in Comola et al. (2019b). See a discussion of the limitations of this scheme in Sects. 5 and S6.3.

Here we show some significant results from the intermittency scheme. Figure 5 shows the global dust emission thresholds in 2006 computed using MERRA-2 fields. Figure 5a shows u∗it, computed using a globally constant Dp0=127μm. Its spatial variability is purely a function of ρa (Eq. 2). u∗it is around  0.16–0.24 m s−1, and higher u∗it implies higher altitude. A lower u∗it leads to a smaller aerodynamic drag force from the airflows given the same wind speed; conversely, there will be a large u∗it for soils over low ρa regions. Figure 5b shows u∗ft (Eq. 1). It varies between 0.2 and 0.9 m s−1. Its spatial variability is dictated by the spatial variability of soil moisture w (see Fig. S1). Regions with the lowest u∗ft are the driest places in the world, which are all deserts. Regions with the highest u∗ft are wet soils covered by rainforests, boreal forests, tundras/permafrosts, and snow. Figure 5c shows the ratio of uft/uit, for which the spatial variability is again dictated by that of w. The magnitude of the ratio conveys not only the strength of the soil moisture effect on the threshold but also the width of the hysteresis regime. Deserts with uft/uit1/Bit have a narrow hysteresis regime (uit<us<uft) and smaller thresholds and thus tend to have more continuous dust emissions. Semiarid and nonarid regions with larger uft/uit tend to have a wide hysteresis regime, and thus dust emissions will be more intermittent.

Figure 5The dust emission thresholds using the Shao and Lu (2000) scheme for the year 2006 on a 0.5× 0.625 grid. (a) The impact threshold u∗it calculated using Dp0=127μm (Eq. 4), (b) the wet fluid threshold u∗ft (Eq. 1), and (c) the ratio between the wet fluid threshold and impact threshold which is fm/0.81, where fm is the moisture effect on u∗ft (Eq. 3). The larger this ratio, the wider the range of wind speeds for which hysteresis in dust emission occurs, and the more important it is to account for intermittency in dust emissions.

Figure 6a shows the 2006 annual mean intermittency effect over the Bodélé Depression as an example. Figure 6a shows hourly mean η as a function of the hourly mean us. It demonstrates the properties of η discussed above: e.g., when us>uft, η→1 and when us<uit, η→0. In both regimes, the behavior of the dust emission intermittency is asymptotic to dichotomous (0 or 1) activity which is the same as that of the conventional emission schemes. Near the intermittency or hysteresis regime uit<us<uft, η is intermediate between zero and 1, and thus a scheme using u∗it gives a small finite emission flux while conventional schemes using u∗ft give a prediction of zero. The color code shows the strength of convection -zi/L. -zi/L is positive (red) when buoyant convection is active (L<0) and is negative (blue) when the PBL is statically stable (L>0). The color shows that there is a modest correlation between -zi/L and us, but the correlation is not necessarily strong, and the strongest buoyancy (dark red) often happens when us (or shear u∗s) is moderate. The strongest buoyancy associates with moderate η values of  0.5 only, and for the highest η values -zi/L is mildly unstable (light red). This means that the turbulent fluctuation is primarily governed by shear u∗s instead of controlled by buoyancy -zi/L, and the intermittency behavior is dictated by shear-driven instead of buoyancy-driven turbulence. Equation (23) could essentially be simplified into σũs121/3us.

Figure 6b shows the global spatial distribution of the annual mean η for the year 2006, averaged across time steps during which saltation and dust emissions are occurring over the grid (and thus η during time steps when Fd=0 is not counted). Most marginal sources have small η<0.3 (red color), indicating dust emissions are fluctuating and intermittent. These emissions may not be existent in other LSMs employing u∗ft in the dust emission equation (but also dependent on their threshold tuning). Over these regions, because of the intermittent shut-offs of emissions, the emissions need to be scaled down by the fraction of time η which is also missed by other LSMs. Intermittency is thus critical in accounting for emissions over semiarid regions. Regions with high η (more continuous emission; light blue color) include El Djouf; the Bodélé Depression; the Libyan–Nubian Desert over Libya, Egypt, and Sudan; the Rub' al Khali Desert within the greater Arabian Desert; the Lut Desert in Iran; the Taklamakan Desert; and the Strzelecki Desert in Australia. During saltation, these regions tend to have wind episodes further away from the thresholds, leading to a high η. However, regions with high η are not necessarily regions with the highest dust emissions (Fig. 7a–b) because their saltation frequency could be low, or the strength of their emission fluxes is limited by other factors such as the fragmentation exponent and soil moisture. In Fig. S5, we further show the factor η averaged over all time steps of 2006, including periods of no emissions over the grid. Since η is close to 0 when there is no emission, Fig. S5 shows much smaller η compared to Fig. 6b.

Figure 6Simulated dust emission intermittency effect for the year 2006. (a) The hourly mean intermittency factor η versus hour mean wind speed us (m s−1) at 0.1 m height over the Bodélé Depression. The color indicates -zi/L, with red indicating -zi/L>0 (unstable), blue indicating -zi/L<0 (stable), and gray indicating -zi/L0 (neutral). The two vertical dotted black lines indicate the annual mean impact threshold (left) and the fluid threshold (right) wind speed (m s−1) at 0.1 m. (b) The annual mean intermittency factor η, or equivalently the fraction of time within a time step that emission is active, averaged over times when emission is active (Fd>0).

4 Results of our new dust emission scheme

In this section, we implement the three new parameterizations of key dust emissions processes with the K14 model into R to investigate the resulting spatial variability of dust emissions. The MERRA-2 data and ESA CCI land cover data are for the year 2006. Other inputs in Table 1 (SoilGrids data, Prigent's roughness, and source functions) are 2-D spatial datasets with no 2006 data available, but they are slowly time-varying variables and are generally used for present-day simulations in other years. In Sect. 4.1, we then compute the dust emissions and analyze their spatial characteristics. We examine the effects of each new modification on the simulated dust emissions by conducting different simulations with different individual parameterizations added. Then, in Sect. 4.2, we inspect the effects of the grid resolution of input data on simulating dust emissions and propose a simple method to calibrate the spatial variability of low-resolution dust emissions to match the spatial variability of high-resolution emissions.

4.1 Effects of different new physics on global dust emissions

In this subsection, we show the effects of different modifications on the resulting dust emissions (in kg m−2 yr−1) in Fig. 7. We demonstrate the effect of each modification by creating a suite of sensitivity experiments as follows: (I) first we simulate emissions using the default K14 scheme, (II) then we use the K14 scheme with our proposed globally constant soil diameter of Dp0127µm for simulation, (III) we further add in the hybrid drag partition physics on top of (II), (IV) we switch from using the fluid threshold to the impact threshold in (III), and finally (V) we include the intermittency effect on top of (IV) for simulation. We note that past studies derived the approximate magnitude of the global total emission (e.g., Tegen and Fung, 1995; Zender et al., 2004; Evan et al., 2014; Kok et al., 2021), but there are no global observations of dust emissions. In the field of dust modeling, there are currently no first principles that can derive the essential dust emission proportionality constants to constrain modeled emissions at a correct order of magnitude, which means scientists still have insufficient knowledge in aeolian physics to generate emissions predictions in the correct order of magnitude. Recognizing that the spatiotemporal characteristics of the predictions are more credible, it is very common for dust modelers to rescale the emissions according to the known constraints of observed atmospheric dust mass or the global DAOD. For instance, Li et al. (2022) scaled all their simulations to achieve a global mean DAOD of 0.03, based on Ridley et al. (2016). Thus, what matters the most is how each modification changes the spatial variability and the relative magnitudes of dust emitted from one region compared to the others, and the absolute magnitude changes are of secondary importance. In our experiments, we normalize all simulations to a global total of 5000 Tg yr−1, which is around the current constraint of global PM20 dust emission flux from past studies (Evan et al., 2014; Kok et al., 2021a, b). Figure 7 shows the normalized emissions of K14 and our scheme (Fig. 7a–b) and differences between the normalized emissions from one modification to another (Fig. 7c–j). The left panels (Fig. 7c, e, g, i) show the normalized emission differences, and the right panels (Fig. 7d, f, h, j) show the normalized emission ratios. The left columns show the regions with the greatest changes in absolute magnitude, which would mostly be dominated by hyper-arid regions and primary sources, whereas the right columns show the regions with the biggest percentage changes with respect to their own order of magnitude. Figure S6 shows the original, unnormalized emission maps for all experiments, and Fig. S7 shows the differences in unnormalized emissions between different experiments. Figure S8 shows the normalized emission maps for all experiments. Table S2 summarizes the global total changes in the normalized emissions done by all modifications.

Figure 7a shows the normalized emissions from the default K14 scheme (experiment I) using MERRA-2 inputs. The emission map shows a similar spatial pattern compared with the K14 simulations in Kok et al. (2014b) and Li et al. (2022) despite differences in the data used for different land-surface fields (e.g., for LAI and soil moisture) and a larger LAI limit of 1 used in this study. The most significant emissions are over the Bodélé Depression in Chad, the Nubian Desert in Sudan and Egypt, the whole Arabian Peninsula, most of Iran, the Taklamakan Desert in China, and the Strzelecki Desert in Australia. Over regions with emissions, the spatial variability of the emissions is dictated by the moisture w (see Figs. S1 and 5b), which fundamentally shapes the threshold u∗ft, the soil erodibility Cd, and the intermittency η. Smaller emissions occur over southern Africa, South America, and the western United States. Without normalization, the simulated emissions in Fig. S6a gives a global total of  29 300 Tg yr−1 using the inputs in Table 1.

Figure 7c–d show the effect of changing the global soil diameter from the current standard of 75 µm to our new constraint of Dp0=127µm (experiment II). As described in Sect. 3.1, a globally larger Dp leads to heavier particles, resulting in higher thresholds and lower emissions across the globe. Figure 7c shows that the normalized emission tends to redistribute from semiarid regions to hyper-arid regions, but the changes are small overall. The color bar in Fig. 7c shows that the effect of employing a new global Dp0 is relatively mild (of the largest order of 0.001 kg m−2 yr−1) compared to the drag partition effect and intermittency (Fig. 7e and g; of the largest order of  0.1 kg m−2 yr−1). Figure 7d shows the ratio map of the normalized emissions of experiment II to those of experiment I. Employing a new globally constant D¯p0 does not strongly impact the spatial pattern of the emissions, so the ratio map is of order 1 around the globe. While Fig. 7c more clearly shows that the largest emission changes in magnitude occur over the major sources, Fig. 7d shows that, after rescaling, there are some stronger emissions reductions in percentage (in blue) over marginal regions (e.g., the Arctic) compared to the minimal increases over the major sources (in white and light red, e.g., over the Bodélé Depression). The major sources are less affected, since in the high u∗s regime Fd becomes more sensitive to u∗s than to u∗ft. Thus, a uniform increase in u∗ft around the globe tends to eliminate small emissions more than large ones. The main effect of employing a larger Dp is therefore a very modest shift of emissions from the marginal regions toward the arid areas. Summing up the absolute magnitude of changes in Fig. 7c, out of 5000 Tg yr−1 there are 250 Tg yr−1 of dust redistributed within source regions. Figure S6b shows the unnormalized global emission flux of  23 900 Tg yr−1, which is 18 % less than the unnormalized emission flux of experiment I as a result of globally increased thresholds. We note that the difference maps of unnormalized emissions in Fig. S7a show that larger emissions reductions occur over the major sources just because the emission magnitudes are larger there.

Figure 7e–f show the effect of including the hybrid drag partition effect Feff (experiment III). The clear contrast between major and marginal sources is shown in Fig. 7e and f, which mirrors the spatial pattern of Feff in Fig. 4c. Compared with experiment II, experiment III has emissions redistributed from semiarid to hyper-arid regions, since the drag partition effect leads to stronger inhibitions of emissions over the nonarid regions than the arid regions. For example, normalized emissions increase (in red) over major sources such as the Bodélé Depression, El Djouf, and the Rub' al Khali Desert; normalized emissions over the western United States and western Australia are significantly reduced. El Djouf is a significant dust source over Africa (Yu et al., 2018), yet K14 fails to represent its high emissions because of the strong moisture effect (see Fig. 5c) compared to other major sources such as the Bodélé Depression. However, Feff highlights El Djouf as a highly erodible surface and helps mitigate the low emission issue over there. Similarly, the underrepresented emissions over the Taklamakan and the Arabian Desert by K14 are partially mitigated by accounting for the drag partitioning (light red in Fig. 7f). Summing up the absolute magnitude of changes in Fig. 7e, out of 5000 Tg yr−1 there is 3611 Tg yr−1 of dust redistributed within the source regions. This shows that the drag partition effect has a much bigger influence on the spatial variability of dust emissions than changing Dp in Fig. 7c. For unnormalized emissions, the global total emission in experiment III decreases drastically by 85 % to  2880 Tg yr−1 relative to that of experiment II. In Fig. S7b, many significant emissions reductions (in dark blue) occur over the Sahara where Feff<0.7, such as Egypt, Sudan, and the western Sahara. K14 struggles to distinguish major sources (e.g., the Bodélé) from less significant sources (e.g., Sudan) and predicts similar levels of emissions. Feff effectively introduces the effect of surface roughness on mitigating emissions over secondary sources, reducing the emissions by at least 1 order of magnitude compared to Fig. S6b.

Figure 7g–h show the effects of implementing the C19 intermittency scheme. It consists of employing u∗it in K14 (experiment IV) and further multiplying the dust emission flux by the intermittency factor η (experiment V). Figure 7g–h show their combined effects by comparing experiment V with experiment III. As seen from Fig. 7g and h, emission schemes employing u∗it will have much stronger emissions over marginal sources. This is because not only is uft/uit bigger but the fragmentation exponent κ (which scales with u∗ft) is also greater over marginal sources. As a result, the main feature in the spatial pattern of Fig. 7h is that the marginal sources (in red color) now have more emissions than experiments II and III. The most apparent emissions increases are over Patagonia, the Horn of Africa (HOA), and western US deserts. Another observable change is that there are many more high-latitude dust emissions, such as over the Arctic, Canada, and Alaska. Studies reported that many models greatly underestimate high-latitude dust emissions (Bullard et al., 2016; Meinander et al., 2022), and the use of u∗it will mitigate this issue. Figure 7g shows that although the changes over nonarid regions are high in ratios (e.g., >1000 times over Canada in Fig. 7h), the emission redistributions in magnitude are still small (e.g.,  10−5 kg m−2 yr−1 in Canada) compared to the major sources (e.g., >10-3 kg m−2 yr−1 in the Sahara) because the magnitudes in experiment III are too small over nonarid regions. The unnormalized total emissions (Fig. S6d) vastly increase to  13 200 Tg yr−1 when employing u∗it instead of u∗ft, more than 4 times that of experiment III.

On the other hand, the effect of multiplying the emission flux with the intermittency factor η is less dramatic than the effect of using u∗it. η tends to scale down small emissions (Fig. 5b), so fluxes from major sources are only moderately reduced. In contrast, fluxes from marginal source regions (e.g., high-latitude boreal forests) are typically reduced by  1–3 orders of magnitude (blue areas in Fig. S6e). However, experiment V has a global total of  11 700 Tg yr−1, which is only 11 % smaller than experiment IV, because those remote regions as such already have small emissions. All in all, accounting for both the impact threshold and the intermittency factor will increase the global total emission from  2880 Tg yr−1 in experiment III to  11 700 Tg yr−1 in experiment V, which is about a 4-fold increase. Figure 7h shows that C19 mainly increases marginal emissions; the overall effect of C19 is thus to move emissions from the hyper-arid regions to semiarid regions. Summing up the absolute magnitude of changes in Fig. 7g, out of 5000 Tg yr−1 there is 3163 Tg yr−1 of dust redistributed within the source regions, indicating that the intermittency scheme induces a similar magnitude of changes compared to employing Feff. Both the hybrid drag partition scheme and the intermittency scheme lead to >60 % of dust emissions redistributed, showing that both effects modify the modeled emission behavior much more strongly compared to the effect of changing the value of Dp (experiment II).

Figure 7b shows the final emission map of our new dust emission scheme with all new physics, and Fig. 7i–j shows the resulting emission changes and ratios from K14 (experiment I) to our scheme (experiment V). Figure 7i–j show that compared to K14, our scheme's emission fluxes over densely vegetated regions (e.g., equatorial Africa and northern Australia) are reduced due to the drag partition effect, while there are increases in marginal sources like the Arctic and mid-latitude boreal forests due to the intermittency effect. The figures show essentially a combination of drag partition (Fig. 7e–f) and intermittency (Fig. 7g–h) effects. Major sources are more affected by the drag partition effect (e.g., the Bodélé Depression and El Djouf), while marginal sources are more dominated by turbulence and intermittency (e.g., the Arctic). For regions where both effects take place, more vegetated semiarid regions are more affected by the drag partition (e.g., western United States), while less vegetated semiarid regions are more affected by the intermittency (e.g., Patagonia, the Great Plains of the United States, and southern Australia). For unnormalized emissions, our scheme's global total of 11 700 Tg yr−1 is  60 % smaller than the K14 emission. A notable feature is that the new mechanisms favor the emissions over the Horn of Africa (HOA) the most, with an emission increase of  2 kg m−2 yr−1 (as seen in Figs. 7g and S7f). This is because in C19, the HOA has low u∗it, high summertime u, low roughness element cover (and thus high Feff), and moderate soil moisture (and thus high dust emission coefficient Cd). This issue could be problematic, since it could introduce too much dust in the GCM over the HOA, which is further discussed in Sect. 5.

Figure 7The effects of the proposed improvements to the parameterization of dust emissions on the default Kok et al. (2014a, b) dust emission scheme. (a, b) Globally normalized dust emission fluxes simulated by (a) the default K14 scheme (expt. (experiment) I) and (b) our new scheme (expt. V). (c–j) Maps of normalized emission (c, e, g, i) differences and (d, f, h, j) ratios, with individual improvements added on top of the default K14 scheme. The individual improvements are respectively (c, d) changing the soil median diameter to 127 µm (expt. II), (e, f) including the drag partition effect (expt. III), and (g, h) employing the Comola et al. (2019) intermittency scheme (expt. V). (i, j) Maps of (i) normalized emission differences and (j) emission ratios for our new scheme and the K14 scheme. The color bars of the maps of differences are drawn to log10 scale.

4.2 The grid-scale dependence of our new dust emission scheme

The spatial resolution of a GCM strongly affects the budget and spatiotemporal variability of dust emissions because modeled emissions scale nonlinearly with input meteorological fields (Ridley et al., 2013; Meng et al., 2021). Many current GCMs have a horizontal resolution of  1–2 (Zhao et al., 2022). Most of the datasets employed in this paper, such as the 0.5 MERRA-2 fields and the even finer aeolian roughness length z0a, are datasets of higher horizontal resolutions. Thus, GCMs need to regrid the datasets to the model native grid resolution as model input. Since dust emission has nonlinear dependences on multiple variables such as u and w, using a simple, area-weighted spatial average of u to calculate dust emission would be inaccurate, as it is different from an area-weighted average of high-resolution dust emissions per se, i.e., for any n>1 usn<usn leads to Fdusn<Fd(usn) (Ridley et al., 2013) (here we use the bar as a spatial average from a finer to a coarser grid). This inequality applies to our model and datasets as well: simply taking the area-weighted mean of high-resolution u∗s or w in a model grid box will omit the locally high u∗s (or low w) values that can produce locally extremely high emissions, resulting in an underestimation of emissions relative to direct area-weighted averaging of the emissions. Additionally, the presence of thresholds in dust emission parameterizations further intensifies the scale dependence, because spatially averaging u∗s might cause us<uit which leads to zero emission over a coarse grid box, whereas fine emissions could be large than zero when us>uit in any fine grid. Dust emission flux will thus be strongly dependent on model resolution more than other linear processes (Zender et al., 2003a; Ridley et al., 2013; Meng et al., 2021), which is undesirable. There is a need to better upscale low-resolution dust emissions to match the variability of high-resolution emissions such that dust emissions tend to be less resolution dependent. To address the problem of missing emissions due to the smoothing of the subgrid wind maxima, a common approach is to employ a grid-by-grid Weibull distribution to the GCM winds to represent the subgrid wind maxima and thus obtain the subgrid emission peaks (Cakmur et al., 2004; Grini et al., 2005; Cowie et al., 2015; Zhang et al., 2016; Menut, 2018; Tai et al., 2021). However, the shape factor of the distribution needs to be empirically determined and thus might not capture the interannual variability and changes in climate. In this subsection, we examine the scale dependence of our dust emission scheme and then propose an alternative approach, which is to derive a simple spatial map and upscale the spatial variability of dust emissions from low-resolution ones to high-resolution ones. A discussion about our proposed approach versus the more common Weibull distribution approach is detailed in Sect. S6.4.

We first examine the scale dependence of our dust emission scheme. We achieve this by performing an area-weighted mean of all input meteorological and land-surface variables to various coarser resolutions. Starting from the native 0.5× 0.625 resolution of MERRA-2, we regrid fields to 0.9× 1.25, 1.9× 2.5, and 4× 5. Then we use these input fields with our new scheme to compute hourly dust emissions for 2006 across these four resolutions. We then compare the emission outputs across these resolutions.

We examine the global, regional, and grid-level scale dependence of our dust emission scheme in Fig. 8. At the regional level, Fig. 8a shows the unnormalized annual total emissions (in Tg yr−1) of nine major dust source regions. The source regions are defined following Kok et al. (2021a) (see Figs. S9 and 8c–d). The regional dust abundance in Fig. 8a is mainly consistent with the regional dust emissions in Kok et al. (2021b; see their Fig. 2). The highest emissions occur over the Middle East/Central Asia, followed by the Sahel and northern Africa. Smaller emissions are over East Asia, North and South America, Australia, and southern Africa. The dust emission scheme also shows scale dependence across different resolutions. Some regions may have a sharper decrease in emissions from 0.9× 1.25 to 1.9× 2.5 (e.g., northeastern Africa), and some regions may have a sharper decrease from 1.9× 2.5 to 4× 5 (e.g., Sahel). This difference is contingent upon the degree of smoothing of the input fields such as u∗s at a particular resolution. The emissions will drastically drop when the local extrema of u∗s are smoothed out and no longer can be represented at a particular resolution, and over different regions this cutoff may occur at different resolutions depending on the spatial heterogeneity of the local-scale meteorological fields. Nevertheless, in general, the coarser the resolution, the worse the model can represent the local variability of u∗s and other input fields and subsequently the emissions, and thus the magnitudes of the emissions decrease with resolution. Figure 8a shows that the relative differences in regional emissions can be different in different resolutions (e.g., Sahel can have larger emissions than the Middle East/Central Asia in the 4× 5 simulation), which will subsequently affect the spatial variability of other major dust cycle variables (such as DAOD or deposition). Therefore, it is always preferential for GCMs to simulate the dust cycle in high resolutions. Figure 8b shows the scale dependence of the global emissions, which also shows the same decrease in emissions from fine to coarse grid resolutions.

We also examine the scale dependence of dust emissions at the grid level. Figure 8c–d show the spatial distributions of 0.5× 0.625 and 1.9× 2.5 unnormalized emissions (with a global total of 11 700 and 5450 Tg yr−1, respectively). The 0.5× 0.625 simulation shows a more detailed local spatial variability compared to the coarser 1.9× 2.5 simulation. The 1.9× 2.5 simulation fails to capture some of the high emission regimes in the 0.5× 0.625 simulation such as over the Taklamakan, and it fails to simulate the local emission peaks over marginal sources such as the Chihuahuan Desert and Patagonia. Coarse-gridded simulations lose emissions in both major and marginal dust sources. We calculate in Fig. 8e the grid-by-grid ratio of unnormalized 0.5× 0.625 emissions to 1.9× 2.5 emissions to show the emissions missed by the 1.9× 2.5 simulation in each grid. Both the coarse (y axis) and fine (x axis) emissions are plotted in the log10 scale to show the emission ratios in terms of the order of magnitude only. It can be seen that most of the points are above the 1:1 line (dashed black line), meaning the coarse emissions are underestimating some local dust fluxes accounted for by the fine emissions. The low-resolution simulation can capture the largest emissions well (top righthand corner), as predicted by the high-resolution simulation. However, the smaller the emission, the more significant the difference in emissions, and for minimal emissions (<10-5 kg m−2 yr−1), the difference can be up to 3 orders of magnitude. There are very few grids with low-resolution emissions higher than high-resolution emissions. Those are exceptional cases due to the spatial variability of moisture w more heterogeneous than that of u∗s, leading to the smoothing of local maxima of u∗ft instead of u∗s and thus the overestimation of Fd in the coarser simulation. In conclusion, the coarse-resolution models omit many local input features and thus fail to represent the correct spatial variability of dust emissions.

Figure 8The dependence on horizontal resolution of dust emissions simulated with our new dust emission scheme. (a, b) Bar plots of dust emissions as a function of grid resolutions for (a) nine major dust emission regions and (b) the globe. (c) Unnormalized emissions of 0.5× 0.625 and (d) 1.9× 2.5. (e) A scatterplot of 0.5× 0.625 versus 1.9× 2.5 unnormalized dust emissions. The rectangular boxes show the nine source regions in (a) defined as in Fig. S9.

To mitigate the scale dependence of our dust emission simulations, here we propose a simple approach to upscale the simulated low-resolution emissions to match the high-resolution emissions. This approach assumes that the fine emissions have a more adequate magnitude and spatiotemporal variability than the coarse emissions so that the coarse emissions are calibrated to match the fine emissions. Our approach is that, by dividing the normalized fine-resolution emission map Fd,f by the coarse-resolution emission map Fd,c, we obtain a map of scaling factors K̃c to account for the differences in the spatial variability of dust emissions between high- and low-resolution simulations:

(24) K ̃ c long , lat = F d , f long , lat / F d , c long , lat .

We obtain a global map of correction factors K̃c and apply this map to all simulated coarse emissions from both historical and future simulations, correcting their spatial variability to match that of the high-resolution emissions. We show the seasonal variability of K̃c in Fig. S10. This map of scaling factors contains some temporal and seasonal variability, so it would be preferable to apply a seasonally varying dataset of K̃c. However, K̃c varies modestly across season because the subgrid variability of the meteorological and land-surface fields is partly determined by spatial structures such as orography or land use/land cover, which change across a relatively long enough timescale (decadal to multidecadal or longer). We note that theoretically this approach will fail to work if, over the remote regions, the low-resolution emission is zero throughout the entire simulation period while the high-resolution emission is a small positive definite, since K̃c will go to infinity. We also note that employing different input fields or different emission schemes will change the subgrid variability and thus the spatial representation of this correction map. For instance, one will obtain a slightly different correction map if one uses ERA-Interim meteorology instead of MERRA-2 or a moderately different map if one employs the Z03 or any other dust emission equations instead of K14 or our new scheme. Therefore, although here we present a standard correction map which is likely accurate and realistic, we suggest each model should make their own correction maps for their specific model configurations. We discuss more caveats of this approach in Sect. 5.4.

Figure 9 shows the ratio maps normalized fine emissions to coarse emissions. Figure 9a shows the ratio of normalized 0.5× 0.625 emissions (Fd,0.5) to constrained 1.9× 2.5 emissions (Fd,2), and Fig. 9b shows the ratio of constrained Fd,0.5 to constrained 0.9× 1.25 emissions (Fd,1). Since all emissions are constrained to match the global dust budget, the correction maps display the relative changes in spatial variability between two resolutions. From Fig. 9a–b, it can be seen that 0.5× 0.625 emissions tend to generate relatively fewer emissions (blue color) over the major sources, such as the Sahara, the Arabian Desert, and the Taklamakan Desert. Emissions of 0.5× 0.625 also tend to have relatively more emissions (red color) over the peripheries of the major sources, such as Algeria, Yemen, and the Taklamakan. This is in line with the above discussion, because the high-resolution simulations tend to be more capable of representing the local peaks of u∗s and therefore can more likely pass the thresholds and produce Fd; the low-resolution simulations would miss a lot of marginal emissions because the low-resolution wind will be smaller than the low-resolution threshold and yield a zero Fd. The ratios over the marginal sources can be up to 100 or even more because of the much smaller emissions in low-resolution simulations. The contrast is smaller (with ratios of 0.5–0.9) over major sources where coarser simulations are more capable of representing the large-scale emission fluxes. For this reason, the correction values in Fig. 9a for 1.9× 2.5 (K̃c,2) generally have bigger magnitudes than those in Fig. 9a for 0.9× 1.25 (K̃c,1). These maps indicate that in coarse-gridded simulations, dust emissions are overall underrepresented over marginal sources and overrepresented over major sources. Therefore, we propose implementing these maps into GCMs of  1 or coarser resolution to correct the dust emission spatial variability accordingly. Scaling all simulations across different spatial resolutions to the finest spatial resolution will move our dust emission scheme toward a scale-aware and grid-independent formulation.

Figure 9Gridded maps of a 0.9× 1.25, 1.9× 2.5, and 4× 5 scaling factor to rescale the coarse dust emission simulations to match the spatial variability of high-resolution dust emissions. This is achieved by calculating the ratios of (a) 0.5× 0.625 to 1.9× 2.5 emissions, (b) 0.5× 0.625 to 0.9× 1.25 emissions, and (c) 0.5× 0.625 to 4× 5 emissions. Emissions of all resolutions are constrained and normalized to have a global total of 5000 Tg yr−1 before calculating the ratios.

4.3 Comparison of our dust emission scheme against other dust estimates

To validate the emissions produced by our new scheme, this subsection focuses on comparing the resulting 0.5× 0.625 emissions from our new scheme against other existing emission data. Since there are no globally gridded observations of dust emissions, GCMs and ESMs mostly evaluate their schemes using observable atmospheric dust products such as satellite and ground-based DAOD data, as well as dust surface concentrations and dust deposition flux measurements (Ridley et al., 2012; Kok et al., 2014b; Pu and Ginoux, 2018; Parajuli et al., 2019; Klose et al., 2021). Since this study focuses on simulating dust emissions, not their subsequent transport and deposition, we compare our results against constraints on the fraction of annual dust emission contributed by nine major source regions (Table 1 in Kok et al., 2021b). These constraints on regional emission fluxes were obtained from the Dust Constraints from joint Observational-Modelling-experiMental analysis (DustCOMM) dataset (Kok et al., 2021a, b), which used inverse modeling to combine an ensemble of model simulations of the global dust cycle with constraints on the regional DAOD near major dust source regions (Ridley et al., 2016), the dust size distribution (Adebiyi and Kok, 2020), and dust extinction efficiency (Kok et al., 2017). The DustCOMM constraints on regional dust emissions include error estimates that account for the spread in model results and the uncertainties in the constraints on dust properties and abundance. Comparisons against independent measurements of dust surface concentrations and deposition fluxes indicated that the DustCOMM product is more accurate than GCM simulations and the MERRA-2 dust reanalysis and that uncertainties are realistic (Kok et al., 2021a, b). The total emissions from all nine source regions obtained by the DustCOMM dataset were 4.7 (3.4–9.1) × 103 Tg yr−1 for dust PM20. The global total of  4700 Tg yr−1 is close to the 5000 Tg yr−1 we adopted in this study for normalization, and we again normalized the DustCOMM global emissions to 5000 Tg yr−1. Note that Kok et al. (2021a, b) only constrained the emissions and other dust variables for each broad region, but its subregional spatial distribution of dust is a multimodel mean and thus unconstrained. Therefore, in Fig. 10b, we sum up the emissions of DustCOMM to regional total emissions, which are constrained by the regional DAOD from Ridley et al. (2016). We also sum up the emissions of all other schemes to regional levels to evaluate each scheme's regional spatial variability against DustCOMM. Tables S3 and S4 summarize the global total emissions and regional emissions of DustCOMM and all other schemes.

We first compare the global emission maps between DustCOMM and different schemes. Figure 10a shows the gridded global spatial distribution of DustCOMM dust emissions (Kok et al., 2021a). Here we compare the gridded simulations of K14 (Fig. 7a) and our scheme (Fig. 7b) against the gridded K21 DustCOMM emissions. Our new scheme's simulation successfully captures most of the major peaks in DustCOMM emissions, except that there are more northern US and high-latitude emissions in our new scheme which were not represented and constrained in DustCOMM's inverse analysis. Our scheme and DustCOMM emissions have a gridded spatial correlation coefficient of r=0.71, showing the resemblance of the two emission maps and our scheme's ability in physically capturing the emission peaks. On the other hand, K14 emissions also share a similar spatial distribution with DustCOMM emissions but show more emissions over central Africa, central India, and northern Australia. Its gridded spatial correlation coefficient with DustCOMM is r=0.61, indicating it does not match DustCOMM emissions in spatial variability as well as our scheme does.

We also conducted simulations using the Z03 scheme (Eq. 10) for comparison with our scheme's simulation. The Z03 scheme requires a source function S, and in this study we adopted two popular source functions: one is the Zender et al. (2003b) geomorphological source function (e.g., used in CESM; Oleson et al., 2013) and the other one is the Ginoux et al. (2001) source function (e.g., used in GEOS-Chem; Fairlie et al., 2007). Both source functions are plotted in Fig. 2 in Kok et al. (2014b). Figure S11a shows the simulations of the Z03 scheme with Ginoux et al. (2001) S (Z03–G), which shows almost an identical pattern with the Z03–G scheme in the GEOS-Chem simulations (e.g., top panel of Fig. 1 in Meng et al., 2021). The Z03–G scheme is mostly consistent with the DustCOMM multimodel emissions (r=0.57) but captures more African emissions in the northwest Africa such as Algeria, Morocco, and the western Sahara. On the other hand, the Z03–Z scheme employs a geomorphic S that possesses a spatial distribution of the upstream area where surface runoff is collected, from Zender et al. (2003b). Figure S11b shows the simulations of the Z03–Z scheme, which shows a highly similar pattern with the Z03–Z scheme in the CESM simulations (e.g., Fig. 2a in Li et al., 2022). It captures the emission peaks across the globe but is quite spatially heterogenous, yielding a relatively low r=0.35 with DustCOMM.

We further obtained MERRA-2-simulated dust emissions for comparison. MERRA-2 simulates the dust cycle using the global ozone chemistry aerosol radiation and transport (GOCART) model, which employs the Ginoux et al. (2001) scheme (hereafter G01). We thus use MERRA-2-simulated dust emissions here for a comparison between G01, Z03, K14, our scheme, and DustCOMM. However, since DustCOMM employs GOCART dust as one of the six models for inverse analysis and assimilation (see Table 1 of Kok et al., 2021a), and because both DustCOMM and MERRA-2 use remotely sensed aerosol optical depth, DustCOMM and MERRA-2 emissions show relatively similar spatial distributions. Figure S12a shows MERRA-2 annual mean emissions, which show very similar features of grid-by-grid variability especially over Australia, East Asia, and the Middle East. For the same reason, MERRA-2 emissions have a high correlation of r=0.86 with DustCOMM.

Next, we aggregate the emission maps to regional total emissions to examine their regional variability. In Fig. 10b, we compare the simulated dust emissions summed over each of nine source regions against the regional DustCOMM emissions. We also summed the emissions outside all rectangular boxes as the high-latitude emissions, to yield the 10 data points in Fig. 10b. High-latitude emissions from our scheme (Fig. 7b) mainly include emissions from Alaska, Canada, Greenland, and Iceland, and there are no emissions from Antarctica because of the lack of necessary input data there (e.g., soil texture and roughness due to rocks). However, since K21 does not provide emission estimates outside of the nine defined source regions, we compare against the estimate of high-latitude dust emissions from Bullard et al. (2016) (hereafter B16) obtained from GCM results. Their definition of high-latitude emission not only includes emissions from the abovementioned regions but also from Patagonia, so we add Patagonia (39–56 S of S. America) emissions as part of the high-latitude emissions in Fig. 10b. B16 estimated that high-latitude emissions (without error estimates) accounted for 4 %–5 % of their assumed 2000 Tg yr−1 of global total emissions. We normalized their estimate to match our global total of 5000 Tg yr−1, yielding a high-latitude emission range of 200–250 Tg yr−1. We take the average, which is 225 Tg yr−1. For all schemes and datasets discussed here, Table S4 provides a list of regional emission contributions to the global total emission.

Figure 10b shows that our scheme's emissions are in overall better agreement with DustCOMM emissions than the K14 scheme. Some of the most significant differences in emissions between our scheme and K14 are over regions including Australia, North America, and southern Africa, where the vegetation drag partitioning causes strong reductions in winds and emission fluxes (Fig. 7f) from K14 to our scheme. Our scheme's East Asian emission is significantly higher than K14's (also shown in Fig. 7j), primarily due to the switch from using u∗ft to u∗it in the dust emission equation (Figs. 7h and S8d). Emissions over South America, the Middle East, and the three regions of northern Africa have relatively small and negligible differences between K14 and our scheme. This occurs because both the drag partitioning and intermittency effects create only minimal changes to the emissions over these regions (Fig. 7j). The results with our scheme (blue color) better match the DustCOMM regional emissions than results with the K14 scheme, lying substantially closer to the 1:1 (black) line over most regions including Africa, Asia, and Australia. There are two notable exceptions where our scheme has less agreement than K14 with DustCOMM, namely North America and the high-latitude emissions. Our scheme generates fewer dust emissions over the Mojave–Sonoran–Chihuahuan deserts over the United States–Mexico border compared to the K14 emissions (Fig. 7a), because of the high LAI (annual mean >0.4) over the western United States that leads to the strong wind drag partitioning. Meanwhile, our scheme generates significant high-latitude emissions over the Arctic, which were not captured by K14 emissions due to the very high u∗ft. Using the emission intermittency parameterization, our scheme represents one of the earliest attempts to successfully capture significant levels of high-latitude emissions. Our high-latitude emissions account for 262 Tg yr−1 (without Patagonia) and 308 Tg yr−1 (with Patagonia), in total accounting for 5 %–6 % of a global sum of 5000 Tg yr−1, which is very close to the percentage B16 suggested. Our scheme has an R2 of 89 % and a root mean squared error (RMSE) of 141 Tg yr−1. We note in Fig. 10b the normalized RMSE (NRMSE) of 28 %, which is the RMSE divided by the mean of the DustCOMM emissions (5000 Tg yr−1/10 data points = 500 Tg yr−1). Our scheme's performance is substantially better than K14's performance with an R2 of 65 % and an RMSE = 259 Tg yr−1 (NRMSE = 52 %).

On the other hand, the Z03–Z scheme (Fig. 10c in green) has a similar level of performance compared with K14, with a higher RMSE of 317 Tg yr−1 (NRMSE = 63 %) and an R2 of 64 %. The Z03–G simulation (Fig. 10c in violet) has a higher R2 of 83 % and also a smaller RMSE of 237 Tg yr−1 (NRMSE = 43 %) against DustCOMM compared with K14 and Z03–Z. Meanwhile, MERRA-2 (Fig. 10d) has a high regional correlation of R2 of 88 % and RMSE of 187 Tg yr−1 (NRMSE = 37 %) against DustCOMM regional variability. In conclusion, our scheme outperforms all the aforementioned simulations in matching against the DustCOMM estimates of regional dust emissions.

To evaluate our simulations of the dust emission thresholds, we also compare our simulations of dust emission thresholds against observationally based threshold estimates from Pu et al. (2020). They compared reanalyzed wind speed distributions against observationally derived DAOD distributions to obtain a threshold wind speed for each grid box that corresponds to a threshold DAOD value (e.g., 0.5 over arid regions and 0.05 over semiarid regions), above which is defined as a dust emission event. We show that our simulations of dust emission threshold overall match their derived threshold wind speed in magnitude and spatial variability (see Sect. S5 and Fig. S13).

Figure 10Dust emissions simulated using different schemes compared against DustCOMM (K21) and Bullard et al. (2016) constraints on regional dust emissions. (a) Globally gridded DustCOMM emissions (kg m−2 yr−1) based on emissions from six different models that were adjusted using inverse modeling to match constraints on particle size distribution, extinction efficiency, and regional dust aerosol optical depth. The rectangular regions specify the nine source regions defined by Kok et al. (2021a) as also shown in Fig. 8c–d. (b) DustCOMM regional emissions (based on fractional emissions reported in the fifth column of Table 1 in K21b scaled to a global total of 5000 Tg yr−1) versus the regional emissions computed by the K14 scheme and our new scheme. (c, d) Comparison of the DustCOMM regional emissions with the regional emissions computed by the (c) Z03–Z and Z03–G schemes and with (d) MERRA-2 dust emissions. The regional emissions of all simulations are obtained from DustCOMM following the nine source regions in panel (a), with one extra data point representing the high-latitude emissions estimated by Bullard et al. (2016). The error bars show 1 standard error, except that the B16 high-latitude emission does not include an error estimate. The black line shows the 1:1 line.

5 Discussion of the caveats and limitations of the new parameterization

Because of the complexities of simulating the fine-scale process of dust emission in a large-scale gridded model, the parameterizations of dust emission processes presented in Sect. 4 and the dust schemes employed from all past studies necessarily made a number of assumptions. Below, we highlight one important caveat or limitation for each modification made in this paper. We provide discussions on further limitations for each modification in Sect. S6.

For the new dry-soil median diameter Dp representation, we used a globally constant Dp0=127±47µm for calculating dust emission thresholds. In theory, Dp should be a function of soil properties (Hillel, 1980) and therefore implicitly a function of space and time, but we obtained a simple relationship for Dp over arid regions because (1) there were no statistically significant correlations between Dp and standard soil properties like soil texture (Fig. S2), (2) data are insufficient for a detailed statistical analysis against a wider range of soil properties, (3) the uncertainties of soil data are moderately large, and (4) most Dp measurements over arid regions found Dp within 40–250 µm (Fig. 1c), which limited u∗ft0 to vary within the relatively small range of  0.204–0.268 m s−1 (from Eq. 2 assuming ρa=1.225 kg m−3). Our analysis of a compilation of Dp measurements also suggests that the spatial variability of Dp over deserts is relatively small compared to Dp over nonarid regions. Our results thus surprisingly suggest that the Dp parameterization can reasonably be much simplified by using a constant value over arid regions. We also compared our approach in deriving Dp against some other derived Dp maps (e.g., Fig. S14) from previous dust modeling studies in Sect. S6.1.

For our hybrid drag partition scheme (Sect. 3.2), we have proposed to combine the rock and vegetation drag partition effects using a weighted mean approach (Eq. 21). This approach assumed that there is a rock regime and a vegetated regime of certain fractions for every grid box, which has the advantage of smoothing the transitions of roughness from very exposed regions (e.g., the Sahara) to semiarid regions with more vegetation (e.g., the Sahel). By using the weighted mean approach, we separate the treatments of rock and vegetation drag partition effects and avoid dealing with the need of adding z0a of rocks and z0a of plants, which is challenging because roughness length is not an additive quantity. For simplicity, Eq. (21) neglects regimes where the roughness of rocks and plants both contribute substantially to the total aeolian roughness. We also note that ESA CCI vegetation cover fractions are annual values and do not exhibit seasonality. The seasonal variability of the simulated Feff in this study is caused by the temporal variability of the LAI input.

For the intermittency scheme (C19 in Sect. 3.3), we note that it has exponential dependences on u∗s, σũs, u∗it, and u∗ft (see Sect. S2) and is thus very sensitive to the accuracy of the GCM simulations of these four variables. For instance, if the thresholds are overestimated by the employed threshold schemes, not only will emissions be underestimated but η from C19 will also be close to zero and further worsen the low bias of the simulated dust emissions. Therefore, a prerequisite of employing the C19 scheme is that the wind friction velocity u∗s and the thresholds u∗it and u∗ft should be adequately simulated and have reasonable ranges of variability throughout the year. We also note that C19 was designed to be used in climate models that run in the Reynolds-averaged Navier–Stokes (RANS) mode. If climate models are resolving turbulence (i.e., in the LES mode), there is minimal need to use C19 to account for intermittent emissions.

For the dust emission scaling method (Sect. 4.3), we note that the scaling factors neglect seasonal variability, which Fig. S10 indicates is moderate. However, employing an annual scaling map like Fig. 9 will already address a large part of the scale-dependence problem. Although we suggested the use of an annual scaling map, ESMs and CTMs that focus on present-day simulations can also perform multiyear simulations in both high and native (coarser) grid resolutions to obtain their own monthly climatological maps of scaling factors. Afterwards, ESMs only need to read in the climatological monthly scaling maps to rescale the native grid dust emissions every month before passing the dust emissions to the atmospheric model component. We further discussed our scaling approach versus some other existent approaches (e.g., using a Weibull distribution for the winds) in Sect. S6.4.

Finally, all emission maps produced in this paper depend on the accuracy of the representations of the input meteorological fields and land-surface variables in various datasets. Our results are particularly sensitive to the soil moisture simulated by MERRA-2, mainly because dust schemes are very sensitive to soil moisture. As a result, although the final simulation of our scheme outperforms other dust emission schemes employed in this paper (Fig. 10), some features in the dust emissions map are unrealistic. First, for instance, the Australian dust emissions are of a comparable order of magnitude to the East Asian emissions (and even larger than East Asian dust emissions in coarser resolutions, per Fig. 8a), which might be because the soil moisture over Australia is slightly underestimated by MERRA-2. Second, northeastern China has larger emissions than northwestern Chinese deserts (Fig. 7b), and similar spatial variability is also seen in other past studies (e.g, Kok et al., 2014b), which might be due to the stronger friction velocity over northeastern China (annual mean > 0.3 m s−1 in MERRA-2) than over northwestern China (annual mean  0.2 m s−1 in MERRA-2) in the GCMs or input MERRA-2 fields (see Fig. S15 for the spatial distribution of MERRA-2 u over northeastern vs. northwestern China). Third, our scheme generates very high summertime emissions over Sudan and the Horn of Africa because of the very high MERRA-2 u ( 1 m s−1 in the summer), low soil moisture, and aeolian roughness (see Fig. 7i). However, this emission peak is not consistent with dust aerosol optical depth (DAOD) observations. There might be several reasons for the HOA emissions peak, including that (i) the input fields are biased over the HOA, (ii) some unknown mechanisms are responsible for suppressing the HOA dust, and (iii) normalized emissions outside of the HOA are overly suppressed relative to the HOA by the drag partitioning and the intermittency effect.

6 Conclusions and significance of our new parameterization

This study presented a new desert dust emission scheme for GCMs and CTMs. The major advances of our scheme compared with existing schemes are the following: (1) we obtained improved parameterizations for several key aspects of dust emission, (2) these improved parameterizations were informed by multiple observations that constrained critical parameters, and (3) we proposed a method to reduce the grid-resolution dependence of the emission scheme that is a common problem to many other existing schemes.

To achieve these advances, we have implemented the following modifications to the existing dust emission scheme of Kok et al. (2014a, b): our first improvement involved the use of soil particle size distributions from multiple past studies to estimate and constrain the soil median diameter Dp as a critical parameter that determines the dust emission thresholds u∗it (impact) and u∗ft (fluid). We found that over the arid desert regions (LAI < 1), Dp can be approximated as a global constant of 127 ± 47 µm, and over nonarid regions Dp increases linearly with silt and clay content. This finding indicates that past dust modeling approaches which parameterized Dp as a function of soil types can be simplified.

Second, we presented a parameterization of the effects of surface rocks and vegetation on the wind drag partition effects, which is not included by many of the current GCMs and CTMs. In particular, a major advance of our drag partition scheme is that we propose a novel method to combine the effects of rocks and vegetation by getting a weighted mean of both effects according to the globally gridded rock and vegetation land-cover area fractions from land-cover datasets (e.g., Klein Goldewijk et al., 2017; ESA, 2017; Kobayashi et al., 2017). Many dust modeling studies only attempted to include the drag partition effect of either one of these roughness elements, and this study represents one of the earliest attempts along with a few other papers (e.g., Darmenova et al., 2009; Foroutan et al., 2017; Klose et al., 2021) to combine and unify the effects on the wind partition of both kinds of roughness elements. Future work should also account for the time-varying vegetation drag partition effect to further enhance the realistic coupling of dust emissions to vegetation dynamics and variability.

Third, we incorporated the boundary-layer turbulence effects on dust emission intermittency by coupling the intermittency scheme formulated by Comola et al. (2019b) to our Kok et al. (2014b) dust emission schemes. The C19 scheme is formulated based on field measurements of simultaneous high-frequency measurements of sand transport and the turbulent wind. This is one of the first studies to have included the turbulence effects on dust emissions, amongst others which focused more on the turbulence effects on convective dust emissions (e.g., Klose et al., 2014; Li et al., 2014). Including the turbulence-driven intermittency effect is important for marginal dust emission sources where the wind speed is normally below the fluid threshold, e.g., dust emissions from high-latitude regions. The C19 scheme also allows dust emission physics to couple better with boundary-layer dynamics and variability, such that simulated dust will have a day-to-day and seasonal variability that is physically linked to the characteristics of the turbulent boundary layer.

Fourth, we proposed a simple scaling method to reduce the inconsistencies in the spatial distributions of the high-resolution and low-resolution dust emission simulations within an LSM. We propose to rescale the low-resolution dust emissions to match the spatial variability of the high-resolution emissions by comparing the spatial distributions of the high- and low-resolution dust emission maps, thereby obtaining a climatological map of scaling factors. The correction maps can thus be applied to other simulations of similar settings, e.g., experiments with the same meteorological and land-surface inputs but different sea/land ice, ocean, stratospheric, or plant physiological forcings. This approach can alleviate the long-standing problem of grid-resolution dependence and spatial distribution inconsistencies of dust emissions across grid sizes among a GCM (e.g., Ridley et al., 2013). Although grid-scale dependence exists in most physical variables in the GCMs, dust emission is exceptionally vulnerable to the grid-resolution-dependence problem because of the very strong nonlinearity (a power of 3–5 or more) of dust emissions to the meteorological fields.

These new approaches act synergistically to improve dust emission modeling. Our new scheme's dust emission simulation, driven by the MERRA-2 meteorological and land-surface fields, shows higher consistency with the Kok et al. (2021a, b) DustCOMM multimodel mean emissions (Fig. 10), which were observationally constrained by an inverse modeling approach and thus contain a realistic regional distribution of dust emissions. Our scheme shows the best agreement against the multimodel mean dust emissions in terms of regional characteristics with R2=89 %, meanwhile other schemes, such as Kok et al. (2014a, b) and Zender et al. (2003a, b) respectively yielded R2=65 % and R2=64 %. This indicates that adding the missing aeolian physics to the existing emission schemes is critical to correctly capturing the dust emission spatial variability and that our scheme displayed almost identical regional characteristics as the inverted multimodel emissions. Our emission map also shows more distinctively the major dust sources including the Bodélé Depression, El Djouf, the Arabian Desert, the Australian Desert, and Patagonia. In our companion paper (Leung et al., 2023b), we will examine the dust cycle simulations of this newly proposed dust emission scheme in CESM and evaluate its performance with other dust cycle variables such as dust PM concentration, dust AOD, lifetime, and deposition flux.

Finally, we note that although our scheme employs a specific emission threshold scheme (i.e., Shao and Lu, 2000) and a specific dust emission equation from Kok et al. (2014b), the modifications we proposed could be applied to different dust emission schemes. For instance, one could use Iversen and White (1982) as a threshold scheme with our newly proposed soil median diameter Dp. One could also use Ginoux et al. (2001), Tegen et al. (2002), or any other dust emission equation to combine with the Comola et al. (2019b) intermittency scheme and our hybrid drag partition scheme. Therefore, our formulation in this paper is highly versatile and adaptable to most of the existing GCMs and CTMs. As such, the new dust emission parameterization presented here can improve the global dust cycle in most GCMs, ESMs, RCMs, and CTMs.

Code availability

The source code of the dust emission model (in R) is available at (Leung, 2023a).

Data availability

MERRA-2 meteorological fields, land-surface fields, and dust emissions are available at (Earthdata, 2023; last access: 12 August 2022). SoilGrids soil texture and properties can be obtained from (Soilgrids, 2023; last access: 20 December 2021). ESA CCI land cover can be obtained from (ESA, 2023; last access: 21 June 2022). DustCOMM inverse analysis results are available at (UCLA, 2023; last access: 21 June 2022). The satellite-derived aeolian roughness data are available upon contacting Catherine Prigent.


The supplement related to this article is available online at:

Author contributions

JFK and DMLe conceptualized the study. DMLe performed the model development, conducted the simulations, analyzed the simulation results, and conducted the evaluations. DMLe wrote the original paper and plotted all figures under JFK's supervision. LL, GSO, CPG-P, MK, LM, DMLa, and NMM assisted with the conceptualization and model development. CP provided the roughness length satellite retrievals. MC assisted with the implementation of the intermittency scheme. All authors contributed to the paper preparation, discussion, and writing.

Competing interests

The contact author has declared that none of the authors has any competing interests.


Publisher’s note: Copernicus Publications remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.


Danny M. Leung thanks Francesco Comola for the helpful comments on the implementation of the dust emission intermittency scheme. Computing and data storage resources, including the Cheyenne supercomputer (, were provided by the Computational and Information Systems Laboratory (CISL) at the National Center for Atmospheric Research (NCAR).

Financial support

Danny M. Leung and Jasper F. Kok are funded by the National Science Foundation (NSF) (grant nos. 1552519 and 1856389). Longlei Li and Natalie M. Mahowald acknowledge support from the Department of Energy (DOE) DE-SC0021302 and the Earth Surface Mineral Dust Source Investigation (EMIT), a NASA Earth Ventures-Instrument (EVI-4) Mission. Gregory S. Okin is funded by the Army Research Office (award no. W911NF-21-1-0070). Martina Klose has received funding through the Helmholtz Association's Initiative and Networking Fund (grant agreement no. VH-NG-1533). Carlos Pérez García-Pando acknowledges support from the European Research Council under the European Union's Horizon 2020 research and innovation programme (grant no. 773051; FRAGMENT) and the AXA Research Fund (AXA Chair on Sand and Dust Storms based at the Barcelona Supercomputing Center). Marcelo Chamecki is supported by the National Science Foundation (grant no. 1856389).

Review statement

This paper was edited by Susannah Burrows and reviewed by two anonymous referees.


Adebiyi, A. A. and Kok, J. F.: Climate models miss most of the coarse dust in the atmosphere, Sci. Adv., 6, eaaz9507,, 2020. 

Albani, S., Mahowald, N. M., Perry, A. T., Scanza, R. A., Zender, C. S., Heavens, N. G., Maggi, V., Kok, J. F., and Otto-Bliesner, B. L.: Improved dust representation in the Community Atmosphere Model, J. Adv. Model. Earth Syst., 6, 541–570,, 2014. 

Alfaro, S. C. and Gomes, L.: Modeling mineral aerosol production by wind erosion: Emission intensities and aerosol size distributions in source areas, J. Geophys. Res.-Atmos., 106, 18075–18084,, 2001. 

Anderson, R. S.: Saltation of sand: a qualitative review with biological analogy, P. Roy. Soc. Edinb. B, 96, 149–165,, 1989. 

Andreotti, B., Claudin, P., and Pouliquen, O.: Measurements of the aeolian sand transport saturation length, Geomorphology, 123, 343–348,, 2010. 

Arya, S. P. S.: A drag partition theory for determining the large-scale roughness parameter and wind stress on the Arctic pack ice, J. Geophys. Res. (1896–1977), 80, 3447–3454,, 1975. 

Bagnold, R. A.: The Transport of Sand by Wind, The Geogr. J., 89, 409–438,, 1937. 

Bagnold, R. A.: The Physics of Blown Sand and Desert Dunes, 1st Edn., Springer Netherlands, New York, 265 pp.,, 1941. 

Barchyn, T. E., Martin, R. L., Kok, J. F., and Hugenholtz, C. H.: Fundamental mismatches between measurements and models in aeolian sediment transport prediction: The role of small-scale variability, Aeolian Res., 15, 245–251,, 2014. 

Barth, M. C., Rasch, P. J., Kiehl, J. T., Benkovitz, C. M., and Schwartz, S. E.: Sulfur chemistry in the National Center for Atmospheric Research Community Climate Model: Description, evaluation, features, and sensitivity to aqueous chemistry, J. Geophys. Res.-Atmos., 105, 1387–1415,, 2000. 

Bradley, E. F. and Mulhearn, P. J.: Development of velocity and shear stress distribution in the wake of a porous shelter fence, J. Wind Eng. Ind. Aerod., 15, 145–156,, 1983. 

Bullard, J. E., Baddock, M., Bradwell, T., Crusius, J., Darlington, E., Gaiero, D., Gassó, S., Gisladottir, G., Hodgkins, R., McCulloch, R., McKenna-Neuman, C., Mockford, T., Stewart, H., and Thorsteinsson, T.: High-latitude dust in the Earth system, Rev. Geophys., 54, 447–485,, 2016. 

Cakmur, R. V., Miller, R. L., and Torres, O.: Incorporating the effect of small-scale circulations upon dust emission in an atmospheric general circulation model, J. Geophys. Res.-Atmos., 109, D7,, 2004. 

Callot, Y., Marticorena, B., and Bergametti, G.: Geomorphologic approach for modelling the surface features of arid environments in a model of dust emissions: application to the Sahara desert, Geodin. Acta, 13, 245–270,, 2000. 

Chappell, A. and Webb, N. P.: Using albedo to reform wind erosion modelling, mapping and monitoring, Aeolian Res., 23, 63–78,, 2016. 

Chappell, A., Webb, N., Hennen, M., Zender, C., Ciais, P., Schepanski, K., Edwards, B., Ziegler, N., Jones, S., Balkanski, Y., Tong, D., Leys, J., Heidenreich, S., Hynes, R., Fuchs, D., Zeng, Z., Ekström, M., Baddock, M., Lee, J., and Kandakji, T.: Weaknesses in dust emission modelling hidden by tuning to dust in the atmosphere, Geosci. Model Dev. Discuss. [preprint],, 2021. 

Chatenet, B., Marticorena, B., Gomes, L., and Bergametti, G.: Assessing the microped size distributions of desert soils erodible by wind, Sedimentology, 43, 901–911,, 1996. 

Chu, C. R., Parlange, M. B., Katul, G. G., and Albertson, J. D.: Probability density functions of turbulent velocity and temperature in the atmospheric surface layer, Water Resour. Res., 32, 1681–1688,, 1996. 

Claudin, P. and Andreotti, B.: A scaling law for aeolian dunes on Mars, Venus, Earth, and for subaqueous ripples, Earth Planet. Sci. Lett., 252, 30–44,, 2006. 

Comola, F., Gaume, J., Kok, J. F., and Lehning, M.: Cohesion-Induced Enhancement of Aeolian Saltation, Geophys. Res. Lett., 46, 5566–5574,, 2019a. 

Comola, F., Kok, J. F., Chamecki, M., and Martin, R. L.: The Intermittency of Wind-Driven Sand Transport, Geophys. Res. Lett., 46, 13430–13440,, 2019b. 

Cowie, S. M., Marsham, J. H., and Knippertz, P.: The importance of rare, high-wind events for dust uplift in northern Africa, Geophys. Res. Lett., 42, 8208–8215,, 2015. 

Dai, Y., Shangguan, W., Wei, N., Xin, Q., Yuan, H., Zhang, S., Liu, S., Lu, X., Wang, D., and Yan, F.: A review of the global soil property maps for Earth system models, SOIL, 5, 137–158,, 2019. 

Darmenova, K., Sokolik, I. N., Shao, Y., Marticorena, B., and Bergametti, G.: Development of a physically based dust emission module within the Weather Research and Forecasting (WRF) model: Assessment of dust emission parameterizations and input parameters for source regions in Central and East Asia, J. Geophys. Res.-Atmos., 114, D14,, 2009. 

Di Biagio, C., Balkanski, Y., Albani, S., Boucher, O., and Formenti, P.: Direct Radiative Effect by Mineral Dust Aerosols Constrained by New Microphysical and Spectral Optical Data, Geophys. Res. Lett., 47, e2019GL086186,, 2020. 

Dunne, J. P., Horowitz, L. W., Adcroft, A. J., Ginoux, P., Held, I. M., John, J. G., Krasting, J. P., Malyshev, S., Naik, V., Paulot, F., Shevliakova, E., Stock, C. A., Zadeh, N., Balaji, V., Blanton, C., Dunne, K. A., Dupuis, C., Durachta, J., Dussin, R., Gauthier, P. P. G., Griffies, S. M., Guo, H., Hallberg, R. W., Harrison, M., He, J., Hurlin, W., McHugh, C., Menzel, R., Milly, P. C. D., Nikonov, S., Paynter, D. J., Ploshay, J., Radhakrishnan, A., Rand, K., Reichl, B. G., Robinson, T., Schwarzkopf, D. M., Sentman, L. T., Underwood, S., Vahlenkamp, H., Winton, M., Wittenberg, A. T., Wyman, B., Zeng, Y., and Zhao, M.: The GFDL Earth System Model Version 4.1 (GFDL-ESM 4.1): Overall Coupled Model Description and Simulation Characteristics, J. Adv. Model. Earth Syst., 12, e2019MS002015,, 2020. 

Dupont, S., Bergametti, G., Marticorena, B., and Simoëns, S.: Modeling saltation intermittency, J. Geophys. Res.-Atmos., 118, 7109–7128,, 2013. 

Durán, O., Claudin, P., and Andreotti, B.: On aeolian transport: Grain-scale interactions, dynamical mechanisms and scaling laws, Aeolian Res., 3, 243–270,, 2011. 

Earthdata: Data Collections, Ges Dic [data set], (last access: 28 April 2023), 2023. 

Elbelrhiti, H., Claudin, P., and Andreotti, B.: Field evidence for surface-wave-induced instability of sand dunes, Nature, 437, 720–723,, 2005. 

ESA: Land Cover CCI Product User Guide Version 2 Technical Report, European Space Agency, (last access: 21 June 2022), 2017. 

ESA: Download CCI LC Products, ESA [data set], (last access: 21 June 2022), 2023. 

Evan, A. T.: Surface Winds and Dust Biases in Climate Models, Geophys. Res. Lett., 45, 1079–1085,, 2018. 

Evan, A. T., Flamant, C., Fiedler, S., and Doherty, O.: An analysis of aeolian dust in climate models, Geophys. Res. Lett., 41, 5996–6001,, 2014. 

Evan, A. T., Fiedler, S., Zhao, C., Menut, L., Schepanski, K., Flamant, C., and Doherty, O.: Derivation of an observation-based map of North African dust emission, Aeolian Res., 16, 153–162,, 2015. 

Evans, S., Ginoux, P., Malyshev, S., and Shevliakova, E.: Climate-vegetation interaction and amplification of Australian dust variability, Geophys. Res. Lett., 43, 11823–11830,, 2016. 

Fairlie, T. D., Jacob, D. J., and Park, R. J.: The impact of transpacific transport of mineral dust in the United States, Atmos. Environ., 41, 1251–1266,, 2007. 

AO/IIASA/ISRIC/ISS-CAS/JRC: Harmonized World Soil Database (version 1.2), FAO, Rome, Italy and IIASA, Laxenburg, Austria, 2012. 

Fécan, F., Marticorena, B., and Bergametti, G.: Parametrization of the increase of the aeolian erosion threshold wind friction velocity due to soil moisture for arid and semi-arid areas, Ann. Geophys., 17, 149–157,, 1999. 

Feng, Y., Wang, H., Rasch, P. J., Zhang, K., Lin, W., Tang, Q., Xie, S., Hamilton, D. S., Mahowald, N., and Yu, H.: Global Dust Cycle and Direct Radiative Effect in E3SM Version 1: Impact of Increasing Model Resolution, J. Adv. Model. Earth Syst., 14, e2021MS002909,, 2022. 

Foroutan, H., Young, J., Napelenok, S., Ran, L., Appel, K. W., Gilliam, R. C., and Pleim, J. E.: Development and evaluation of a physics-based windblown dust emission scheme implemented in the CMAQ modeling system, J. Adv. Model. Earth Syst., 9, 585–608,, 2017. 

Froyd, K. D., Yu, P., Schill, G. P., Brock, C. A., Kupc, A., Williamson, C. J., Jensen, E. J., Ray, E., Rosenlof, K. H., Bian, H., Darmenov, A. S., Colarco, P. R., Diskin, G. S., Bui, T., and Murphy, D. M.: Dominant role of mineral dust in cirrus cloud formation revealed by global-scale measurements, Nat. Geosci., 15, 177–183,, 2022. 

Gelaro, R., McCarty, W., Suárez, M. J., Todling, R., Molod, A., Takacs, L., Randles, C. A., Darmenov, A., Bosilovich, M. G., Reichle, R., Wargan, K., Coy, L., Cullather, R., Draper, C., Akella, S., Buchard, V., Conaty, A., Silva, A. M. da, Gu, W., Kim, G.-K., Koster, R., Lucchesi, R., Merkova, D., Nielsen, J. E., Partyka, G., Pawson, S., Putman, W., Rienecker, M., Schubert, S. D., Sienkiewicz, M., and Zhao, B.: The Modern-Era Retrospective Analysis for Research and Applications, Version 2 (MERRA-2), J. Climate, 30, 5419–5454,, 2017. 

Gillette, D. A. and Passi, R.: Modeling dust emission caused by wind erosion, J. Geophys. Res.-Atmos., 93, 14233–14242,, 1988. 

Ginoux, P., Chin, M., Tegen, I., Prospero, J. M., Holben, B., Dubovik, O., and Lin, S.-J.: Sources and distributions of dust aerosols simulated with the GOCART model, J. Geophys. Res.-Atmos., 106, 20255–20273,, 2001. 

Ginoux, P., Prospero, J. M., Gill, T. E., Hsu, N. C., and Zhao, M.: Global-scale attribution of anthropogenic and natural dust sources and their emission rates based on MODIS Deep Blue aerosol products, Rev. Geophys., 50, 3,, 2012. 

Greeley, R. and Iversen, J. D.: Wind as a Geological Process: On Earth, Mars, Venus and Titan, CUP Archive, 356 pp.,, 1985. 

Greeley, R., Blumberg, D. G., McHone, J. F., Dobrovolskis, A., Iversen, J. D., Lancaster, N., Rasmussen, K. R., Wall, S. D., and White, B. R.: Applications of spaceborne radar laboratory data to the study of aeolian processes, J. Geophys. Res.-Planets, 102, 10971–10983,, 1997. 

Grini, A., Myhre, G., Zender, C. S., and Isaksen, I. S. A.: Model simulations of dust sources and transport in the global atmosphere: Effects of soil erodibility and wind speed variability, J. Geophys. Res.-Atmos., 110, D2,, 2005. 

Guerschman, J. P., Scarth, P. F., McVicar, T. R., Renzullo, L. J., Malthus, T. J., Stewart, J. B., Rickards, J. E., and Trevithick, R.: Assessing the effects of site heterogeneity and soil properties when unmixing photosynthetic vegetation, non-photosynthetic vegetation and bare soil fractions from Landsat and MODIS data, Remote Sens. Environ., 161, 12–26,, 2015. 

Gutman, G. and Ignatov, A.: The derivation of the green vegetation fraction from NOAA/AVHRR data for use in numerical weather prediction models, Int. J. Remote Sens., 19, 1533–1543,, 1998. 

Hamilton, D. S., Moore, J. K., Arneth, A., Bond, T. C., Carslaw, K. S., Hantson, S., Ito, A., Kaplan, J. O., Lindsay, K., Nieradzik, L., Rathod, S. D., Scanza, R. A., and Mahowald, N. M.: Impact of Changes to the Atmospheric Soluble Iron Deposition Flux on Ocean Biogeochemical Cycles in the Anthropocene, Global Biogeochem. Cy., 34, e2019GB006448,, 2020. 

Hengl, T., Jesus, J. M. de, MacMillan, R. A., Batjes, N. H., Heuvelink, G. B. M., Ribeiro, E., Samuel-Rosa, A., Kempen, B., Leenaars, J. G. B., Walsh, M. G., and Gonzalez, M. R.: SoilGrids1km – Global Soil Information Based on Automated Mapping, PLOS ONE, 9, e105992,, 2014. 

Hengl, T., Heuvelink, G. B. M., Kempen, B., Leenaars, J. G. B., Walsh, M. G., Shepherd, K. D., Sila, A., MacMillan, R. A., Jesus, J. M. de, Tamene, L., and Tondoh, J. E.: Mapping Soil Properties of Africa at 250 m Resolution: Random Forests Significantly Improve Current Predictions, PLOS ONE, 10, e0125814,, 2015. 

Hengl, T., Jesus, J. M. de, Heuvelink, G. B. M., Gonzalez, M. R., Kilibarda, M., Blagotiæ, A., Shangguan, W., Wright, M. N., Geng, X., Bauer-Marschallinger, B., Guevara, M. A., Vargas, R., MacMillan, R. A., Batjes, N. H., Leenaars, J. G. B., Ribeiro, E., Wheeler, I., Mantel, S., and Kempen, B.: SoilGrids250m: Global gridded soil information based on machine learning, PLOS ONE, 12, e0169748,, 2017. 

Hillel, D.: Fundamentals of Soil Physics, Elsevier,, 1980. 

Huang, X. and Foroutan, H.: Effects of Non-Photosynthetic Vegetation on Dust Emissions, J. Geophys. Res.-Atmos., 127, e2021JD035243,, 2022. 

Huneeus, N., Schulz, M., Balkanski, Y., Griesfeller, J., Prospero, J., Kinne, S., Bauer, S., Boucher, O., Chin, M., Dentener, F., Diehl, T., Easter, R., Fillmore, D., Ghan, S., Ginoux, P., Grini, A., Horowitz, L., Koch, D., Krol, M. C., Landing, W., Liu, X., Mahowald, N., Miller, R., Morcrette, J.-J., Myhre, G., Penner, J., Perlwitz, J., Stier, P., Takemura, T., and Zender, C. S.: Global dust model intercomparison in AeroCom phase I, Atmos. Chem. Phys., 11, 7781–7816,, 2011. 

Ito, A. and Kok, J. F.: Do dust emissions from sparsely vegetated regions dominate atmospheric iron supply to the Southern Ocean?, J. Geophys. Res.-Atmos., 122, 3987–4002,, 2017. 

Iversen, J. D. and White, B. R.: Saltation threshold on Earth, Mars and Venus, Sedimentology, 29, 111–119,, 1982. 

Jickells, T. D., An, Z. S., Andersen, K. K., Baker, A. R., Bergametti, G., Brooks, N., Cao, J. J., Boyd, P. W., Duce, R. A., Hunter, K. A., Kawahata, H., Kubilay, N., laRoche, J., Liss, P. S., Mahowald, N., Prospero, J. M., Ridgwell, A. J., Tegen, I., and Torres, R.: Global Iron Connections Between Desert Dust, Ocean Biogeochemistry, and Climate, Science, 308, 67–71,, 2005. 

Jin, Q., Wei, J., Lau, W. K. M., Pu, B., and Wang, C.: Interactions of Asian mineral dust with Indian summer monsoon: Recent advances and challenges, Earth-Sci. Rev., 215, 103562,, 2021. 

Kawai, K., Matsui, H., Kimura, R., and Shinoda, M.: High Sensitivity of Asian Dust Emission, Transport, and Climate Impacts to Threshold Friction Velocity, Sola, 17, 239–245,, 2021. 

King, J., Nickling, W. G., and Gillies, J. A.: Representation of vegetation and other nonerodible elements in aeolian shear stress partitioning models for predicting transport threshold, J. Geophys. Res.-Earth Surf., 110, F4,, 2005. 

Kinne, S., Schulz, M., Textor, C., Guibert, S., Balkanski, Y., Bauer, S. E., Berntsen, T., Berglen, T. F., Boucher, O., Chin, M., Collins, W., Dentener, F., Diehl, T., Easter, R., Feichter, J., Fillmore, D., Ghan, S., Ginoux, P., Gong, S., Grini, A., Hendricks, J., Herzog, M., Horowitz, L., Isaksen, I., Iversen, T., Kirkevåg, A., Kloster, S., Koch, D., Kristjansson, J. E., Krol, M., Lauer, A., Lamarque, J. F., Lesins, G., Liu, X., Lohmann, U., Montanaro, V., Myhre, G., Penner, J., Pitari, G., Reddy, S., Seland, O., Stier, P., Takemura, T., and Tie, X.: An AeroCom initial assessment – optical properties in aerosol component modules of global models, Atmos. Chem. Phys., 6, 1815–1834,, 2006. 

Klein Goldewijk, K., Beusen, A., Doelman, J., and Stehfest, E.: Anthropogenic land use estimates for the Holocene – HYDE 3.2, Earth Syst. Sci. Data, 9, 927–953,, 2017. 

Klose, M. and Shao, Y.: Stochastic parameterization of dust emission and application to convective atmospheric conditions, Atmos. Chem. Phys., 12, 7309–7320,, 2012. 

Klose, M. and Shao, Y.: Large-eddy simulation of turbulent dust emission, Aeolian Res., 8, 49–58,, 2013. 

Klose, M., Shao, Y., Li, X., Zhang, H., Ishizuka, M., Mikami, M., and Leys, J. F.: Further development of a parameterization for convective turbulent dust emission and evaluation based on field observations, J. Geophys. Res.-Atmos., 119, 10441–10457,, 2014. 

Klose, M., Jorba, O., Gonçalves Ageitos, M., Escribano, J., Dawson, M. L., Obiso, V., Di Tomaso, E., Basart, S., Montané Pinto, G., Macchia, F., Ginoux, P., Guerschman, J., Prigent, C., Huang, Y., Kok, J. F., Miller, R. L., and Pérez García-Pando, C.: Mineral dust cycle in the Multiscale Online Nonhydrostatic AtmospheRe CHemistry model (MONARCH) Version 2.0, Geosci. Model Dev., 14, 6403–6444,, 2021. 

Kobayashi, T., Tateishi, R., Alsaaideh, B., Sharma, R., Wakaizumi, T., Miyamoto, D., Bai, X., Long, B., Gegentana, G., Maitiniyazi, A., Cahyana, D., Haireti, A., Morifuji, Y., Abake, G., Pratama, R., Zhang, N., Alifu, Z., Shirahata, T., Mi, L., Iizuka, K., Yusupujiang, A., Rinawan, F., Bhattarai, R., and Phong, D.: Production of Global Land Cover Data – GLCNMO2013, J. Geogr. Geol., 9, p1,, 2017. 

Kok, J. F.: An improved parameterization of wind-blown sand flux on Mars that includes the effect of hysteresis, Geophys. Res. Lett., 37, 12,, 2010. 

Kok, J. F., Parteli, E. J. R., Michaels, T. I., and Karam, D. B.: The physics of wind-blown sand and dust, Rep. Prog. Phys., 75, 106901,, 2012. 

Kok, J. F., Mahowald, N. M., Fratini, G., Gillies, J. A., Ishizuka, M., Leys, J. F., Mikami, M., Park, M.-S., Park, S.-U., Van Pelt, R. S., and Zobeck, T. M.: An improved dust emission model – Part 1: Model description and comparison against measurements, Atmos. Chem. Phys., 14, 13023–13041,, 2014a. 

Kok, J. F., Albani, S., Mahowald, N. M., and Ward, D. S.: An improved dust emission model – Part 2: Evaluation in the Community Earth System Model, with implications for the use of dust source functions, Atmos. Chem. Phys., 14, 13043–13061,, 2014b. 

Kok, J. F., Ridley, D. A., Zhou, Q., Miller, R. L., Zhao, C., Heald, C. L., Ward, D. S., Albani, S., and Haustein, K.: Smaller desert dust cooling effect estimated from analysis of dust size and abundance, Nat. Geosci., 10, 274–278,, 2017. 

Kok, J. F., Ward, D. S., Mahowald, N. M., and Evan, A. T.: Global and regional importance of the direct dust-climate feedback, Nat. Commun., 9, 241,, 2018. 

Kok, J. F., Adebiyi, A. A., Albani, S., Balkanski, Y., Checa-Garcia, R., Chin, M., Colarco, P. R., Hamilton, D. S., Huang, Y., Ito, A., Klose, M., Leung, D. M., Li, L., Mahowald, N. M., Miller, R. L., Obiso, V., Pérez García-Pando, C., Rocha-Lima, A., Wan, J. S., and Whicker, C. A.: Improved representation of the global dust cycle using observational constraints on dust properties and abundance, Atmos. Chem. Phys., 21, 8127–8167,, 2021a. 

Kok, J. F., Adebiyi, A. A., Albani, S., Balkanski, Y., Checa-Garcia, R., Chin, M., Colarco, P. R., Hamilton, D. S., Huang, Y., Ito, A., Klose, M., Li, L., Mahowald, N. M., Miller, R. L., Obiso, V., Pérez García-Pando, C., Rocha-Lima, A., and Wan, J. S.: Contribution of the world's main dust source regions to the global cycle of desert dust, Atmos. Chem. Phys., 21, 8169–8193,, 2021b. 

Kok, J. F., Storelvmo, T., Karydis, V. A., Adebiyi, A. A., Mahowald, N. M., Evan, A. T., He, C., and Leung, D. M.: Mineral dust aerosol impacts on global climate and climate change, Nat. Rev. Earth Environ., 4, 1–16,, 2023. 

Koven, C. D. and Fung, I.: Identifying global dust source areas using high-resolution land surface form, J. Geophys. Res.-Atmos., 113, D22,, 2008. 

Laurent, B., Marticorena, B., Bergametti, G., Chazette, P., Maignan, F., and Schmechtig, C.: Simulation of the mineral dust emission frequencies from desert areas of China and Mongolia using an aerodynamic roughness length map derived from the POLDER/ADEOS 1 surface products, J. Geophys. Res.-Atmos., 110, D18,, 2005. 

Lawrence, D. M., Hurtt, G. C., Arneth, A., Brovkin, V., Calvin, K. V., Jones, A. D., Jones, C. D., Lawrence, P. J., de Noblet-Ducoudré, N., Pongratz, J., Seneviratne, S. I., and Shevliakova, E.: The Land Use Model Intercomparison Project (LUMIP) contribution to CMIP6: rationale and experimental design, Geosci. Model Dev., 9, 2973–2998,, 2016. 

LeGrand, S. L., Polashenski, C., Letcher, T. W., Creighton, G. A., Peckham, S. E., and Cetola, J. D.: The AFWA dust emission scheme for the GOCART aerosol model in WRF-Chem v3.8.1, Geosci. Model Dev., 12, 131–166,, 2019. 

LeGrand, S. L., Letcher, T. W., Okin, G. S., Webb, N. P., Gallagher, A. R., Dhital, S., Hodgdon, T. S., Ziegler, N. P., and Michaels, M. L.: Application of a satellite-retrieved sheltering parameterization (v1.0) for dust event simulation with WRF-Chem v4.1, Geosci. Model Dev., 16, 1009–1038,, 2023. 

Leung, D. M., Kok, J. F., Li, L., Okin, G. S., Prigent, C., Klose, M., Garcia-Pando, C. P., Menut, L., Mahowald, N. M., Lawrence, D. M., and Chamecki, M.: Source code for “ new process-based and scale-respecting desert dust emission scheme for globa l climate models - Part I: description and evaluation against inverse modeling emissions” (1.0), Zenodo [code],, 2023a. 

Leung, D. M., Kok, J. F., Li, L., Mahowald, N. M., Lawrence, D. M., Tilmes, S., Kluzek, E., Klose, M., and Pérez García-Pando, C.: A new process-based and scale-aware desert dust emission scheme for global climate models – Part II: evaluation in the Community Earth System Model (CESM2), EGUsphere [preprint],, 2023. 

Li, J., Okin, G. S., Herrick, J. E., Belnap, J., Miller, M. E., Vest, K., and Draut, A. E.: Evaluation of a new model of aeolian transport in the presence of vegetation, J. Geophys. Res.-Earth Surf., 118, 288–306,, 2013. 

Li, L., Mahowald, N. M., Miller, R. L., Pérez García-Pando, C., Klose, M., Hamilton, D. S., Gonçalves Ageitos, M., Ginoux, P., Balkanski, Y., Green, R. O., Kalashnikova, O., Kok, J. F., Obiso, V., Paynter, D., and Thompson, D. R.: Quantifying the range of the dust direct radiative effect due to source mineralogy uncertainty, Atmos. Chem. Phys., 21, 3973–4005,, 2021. 

Li, L., Mahowald, N. M., Kok, J. F., Liu, X., Wu, M., Leung, D. M., Hamilton, D. S., Emmons, L. K., Huang, Y., Sexton, N., Meng, J., and Wan, J.: Importance of different parameterization changes for the updated dust cycle modeling in the Community Atmosphere Model (version 6.1), Geosci. Model Dev., 15, 8181–8219,, 2022. 

Li, W., MacBean, N., Ciais, P., Defourny, P., Lamarche, C., Bontemps, S., Houghton, R. A., and Peng, S.: Gross and net land cover changes in the main plant functional types derived from the annual ESA CCI land cover maps (1992–2015), Earth Syst. Sci. Data, 10, 219–234,, 2018. 

Li, X., Feng, G., Sharratt, B. S., Zheng, Z., Pi, H., and Gao, F.: Soil Wind Erodibility Based on Dry Aggregate-Size Distribution in the Tarim Basin, Soil Sci. Soc. Am. J., 78, 2009–2016,, 2014. 

Lin, S.-J.: A “Vertically Lagrangian” Finite-Volume Dynamical Core for Global Models, Mon. Weather Rev., 132, 2293–2307,<2293:AVLFDC>2.0.CO;2, 2004. 

Linke, C., Möhler, O., Veres, A., Mohácsi, Á., Bozóki, Z., Szabó, G., and Schnaiter, M.: Optical properties and mineralogical composition of different Saharan mineral dust samples: a laboratory study, Atmos. Chem. Phys., 6, 3315–3323,, 2006. 

Liu, H., Jacob, D. J., Bey, I., and Yantosca, R. M.: Constraints from 210Pb and 7Be on wet deposition and transport in a global three-dimensional chemical tracer model driven by assimilated meteorological fields, J. Geophys. Res.-Atmos., 106, 12109–12128,, 2001. 

MacKinnon, D. J., Clow, G. D., Tigges, R. K., Reynolds, R. L., and Chavez, P. S.: Comparison of aerodynamically and model-derived roughness lengths (zo) over diverse surfaces, central Mojave Desert, California, USA, Geomorphology, 63, 103–113,, 2004. 

Mahowald, N. M., Baker, A. R., Bergametti, G., Brooks, N., Duce, R. A., Jickells, T. D., Kubilay, N., Prospero, J. M., and Tegen, I.: Atmospheric global dust cycle and iron inputs to the ocean, Global Biogeochem. Cy., 19,, 2005. 

Mahowald, N. M., Muhs, D. R., Levis, S., Rasch, P. J., Yoshioka, M., Zender, C. S., and Luo, C.: Change in atmospheric mineral aerosols in response to climate: Last glacial period, preindustrial, modern, and doubled carbon dioxide climates, J. Geophys. Res.-Atmos., 111, D10, 2006. 

Mahowald, N. M., Kloster, S., Engelstaedter, S., Moore, J. K., Mukhopadhyay, S., McConnell, J. R., Albani, S., Doney, S. C., Bhattacharya, A., Curran, M. A. J., Flanner, M. G., Hoffman, F. M., Lawrence, D. M., Lindsay, K., Mayewski, P. A., Neff, J., Rothenberg, D., Thomas, E., Thornton, P. E., and Zender, C. S.: Observed 20th century desert dust variability: impact on climate and biogeochemistry, Atmos. Chem. Phys., 10, 10875–10893,, 2010. 

Mailler, S., Menut, L., Khvorostyanov, D., Valari, M., Couvidat, F., Siour, G., Turquety, S., Briant, R., Tuccella, P., Bessagnet, B., Colette, A., Létinois, L., Markakis, K., and Meleux, F.: CHIMERE-2017: from urban to hemispheric chemistry-transport modeling, Geosci. Model Dev., 10, 2397–2423,, 2017. 

Marticorena, B. and Bergametti, G.: Modeling the atmospheric dust cycle: 1. Design of a soil-derived dust emission scheme, J. Geophys. Res.-Atmos., 100, 16415–16430,, 1995. 

Marticorena, B., Bergametti, G., Aumont, B., Callot, Y., N'Doumé, C., and Legrand, M.: Modeling the atmospheric dust cycle: 2. Simulation of Saharan dust sources, J. Geophys. Res.-Atmos., 102, 4387–4404,, 1997. 

Marticorena, B., Chazette, P., Bergametti, G., Dulac, F., and Legrand, M.: Mapping the aerodynamic roughness length of desert surfaces from the POLDER/ADEOS bi-directional reflectance product, Int. J. Remote Sens., 25, 603–626,, 2004. 

Marticorena, B., Kardous, M., Bergametti, G., Callot, Y., Chazette, P., Khatteli, H., Hégarat-Mascle, S. L., Maillé, M., Rajot, J.-L., Vidal-Madjar, D., and Zribi, M.: Surface and aerodynamic roughness in arid and semiarid areas and their relation to radar backscatter coefficient, J. Geophys. Res.-Earth Surf., 111, F3,, 2006. 

Martin, R. L. and Kok, J. F.: Distinct Thresholds for the Initiation and Cessation of Aeolian Saltation From Field Measurements, J. Geophys. Res.-Earth Surf., 123, 1546–1565,, 2018. 

Martin, R. L. and Kok, J. F.: Size-Independent Susceptibility to Transport in Aeolian Saltation, J. Geophys. Res.-Earth Surf., 124, 1658–1674,, 2019. 

McGlynn, I. O. and Okin, G. S.: Characterization of shrub distribution using high spatial resolution remote sensing: Ecosystem implications for a former Chihuahuan Desert grassland, Remote Sens. Environ., 101, 554–566,, 2006. 

McGraw, Z., Storelvmo, T., David, R. O., and Sagoo, N.: Global Radiative Impacts of Mineral Dust Perturbations Through Stratiform Clouds, J. Geophys. Res.-Atmos., 125, e2019JD031807,, 2020. 

McKenna Neuman, C. and Sanderson, S.: Humidity control of particle emissions in aeolian systems, J. Geophys. Res.-Earth Surf., 113, F2,, 2008. 

Meier, R., Davin, E. L., Bonan, G. B., Lawrence, D. M., Hu, X., Duveiller, G., Prigent, C., and Seneviratne, S. I.: Impacts of a revised surface roughness parameterization in the Community Land Model 5.1, Geosci. Model Dev., 15, 2365–2393,, 2022. 

Meinander, O., Dagsson-Waldhauserova, P., Amosov, P., Aseyeva, E., Atkins, C., Baklanov, A., Baldo, C., Barr, S. L., Barzycka, B., Benning, L. G., Cvetkovic, B., Enchilik, P., Frolov, D., Gassó, S., Kandler, K., Kasimov, N., Kavan, J., King, J., Koroleva, T., Krupskaya, V., Kulmala, M., Kusiak, M., Lappalainen, H. K., Laska, M., Lasne, J., Lewandowski, M., Luks, B., McQuaid, J. B., Moroni, B., Murray, B., Möhler, O., Nawrot, A., Nickovic, S., O’Neill, N. T., Pejanovic, G., Popovicheva, O., Ranjbar, K., Romanias, M., Samonova, O., Sanchez-Marroquin, A., Schepanski, K., Semenkov, I., Sharapova, A., Shevnina, E., Shi, Z., Sofiev, M., Thevenet, F., Thorsteinsson, T., Timofeev, M., Umo, N. S., Uppstu, A., Urupina, D., Varga, G., Werner, T., Arnalds, O., and Vukovic Vimic, A.: Newly identified climatically and environmentally significant high-latitude dust sources, Atmos. Chem. Phys., 22, 11889–11930,, 2022. 

Meng, J., Martin, R. V., Ginoux, P., Hammer, M., Sulprizio, M. P., Ridley, D. A., and van Donkelaar, A.: Grid-independent high-resolution dust emissions (v1.0) for chemical transport models: application to GEOS-Chem (12.5.0), Geosci. Model Dev., 14, 4249–4260,, 2021. 

Menut, L.: Modeling of Mineral Dust Emissions with a Weibull Wind Speed Distribution Including Subgrid-Scale Orography Variance, J. Atmos. Ocean. Technol., 35, 1221–1236,, 2018. 

Menut, L., Pérez, C., Haustein, K., Bessagnet, B., Prigent, C., and Alfaro, S.: Impact of surface roughness and soil texture on mineral dust emission fluxes modeling, J. Geophys. Res.-Atmos., 118, 6505–6520,, 2013. 

Miller, R. L. and Tegen, I.: Climate Response to Soil Dust Aerosols, J. Climate, 11, 3247–3267,<3247:CRTSDA>2.0.CO;2, 1998. 

Miller, R. L., Tegen, I., and Perlwitz, J.: Surface radiative forcing by soil dust aerosols and the hydrologic cycle, J. Geophys. Res.-Atmos., 109, D4,, 2004. 

Minvielle, F., Marticorena, B., Gillette, D. A., Lawson, R. E., Thompson, R., and Bergametti, G.: Relationship between the Aerodynamic Roughness Length and the Roughness Density in Cases of Low Roughness Density*, Environm. Fluid Mecha., 3, 249–267,, 2003. 

Mokhtari, M., Gomes, L., Tulet, P., and Rezoug, T.: Importance of the surface size distribution of erodible material: an improvement on the Dust Entrainment And Deposition (DEAD) Model, Geosci. Model Dev., 5, 581–598,, 2012. 

Nikuradse, J.: Laws of flow in rough pipes, Tech. Rep. NACATM-1292, National Advisory Committee for Aeronautics, (last access: 20 February 2023), 1950. 

Okin, G. S.: A new model of wind erosion in the presence of vegetation, J. Geophys. Res.-Earth Surf., 113, F2,, 2008. 

Oleson, K., Lawrence, D., Bonan, G., Drewniak, B., Huang, M., Koven, C., Levis, S., Li, F., Riley, W., Subin, Z., Swenson, S., Thornton, P., Bozbiyik, A., Fisher, R., Heald, C., Kluzek, E., Lamarque, J.-F., Lawrence, P., Leung, L., Lipscomb, W., Muszala, S., Ricciuto, D., Sacks, W., Sun, Y., Tang, J., and Yang, Z.-L.: Technical description of version 4.5 of the Community Land Model (CLM), No. NCAR/TN503+STR),, 2013. 

Pähtz, T., Clark, A. H., Valyrakis, M., and Durán, O.: The Physics of Sediment Transport Initiation, Cessation, and Entrainment Across Aeolian and Fluvial Environments, Rev. Geophys., 58, e2019RG000679,, 2020. 

Panofsky, H. A., Tennekes, H., Lenschow, D. H., and Wyngaard, J. C.: The characteristics of turbulent velocity components in the surface layer under convective conditions, Bound.-Lay. Meteorol., 11, 355–361,, 1977. 

Parajuli, S. P., Stenchikov, G. L., Ukhov, A., and Kim, H.: Dust Emission Modeling Using a New High-Resolution Dust Source Function in WRF-Chem With Implications for Air Quality, J. Geophys. Res.-Atmos., 124, 10109–10133,, 2019. 

Parajuli, S. P., Stenchikov, G. L., Ukhov, A., Mostamandi, S., Kucera, P. A., Axisa, D., Gustafson Jr., W. I., and Zhu, Y.: Effect of dust on rainfall over the Red Sea coast based on WRF-Chem model simulations, Atmos. Chem. Phys., 22, 8659–8682,, 2022. 

Petroff, A. and Zhang, L.: Development and validation of a size-resolved particle dry deposition scheme for application in aerosol transport models, Geosci. Model Dev., 3, 753–769,, 2010. 

Pierre, C., Bergametti, G., Marticorena, B., Kergoat, L., Mougin, E., and Hiernaux, P.: Comparing drag partition schemes over a herbaceous Sahelian rangeland, J. Geophys. Res.-Earth Surf., 119, 2291–2313,, 2014a. 

Pierre, C., Bergametti, G., Marticorena, B., AbdourhamaneTouré, A., Rajot, J.-L., and Kergoat, L.: Modeling wind erosion flux and its seasonality from a cultivated sahelian surface: A case study in Niger, CATENA, 122, 61–71,, 2014b. 

Prigent, C., Tegen, I., Aires, F., Marticorena, B., and Zribi, M.: Estimation of the aerodynamic roughness length in arid and semi-arid regions over the globe with the ERS scatterometer, J. Geophys. Res.-Atmos., 110, D9,, 2005. 

Prigent, C., Jiménez, C., and Catherinot, J.: Comparison of satellite microwave backscattering (ASCAT) and visible/near-infrared reflectances (PARASOL) for the estimation of aeolian aerodynamic roughness length in arid and semi-arid regions, Atmos. Meas. Tech., 5, 2703–2712,, 2012. 

Prospero, J. M.: Long-range transport of mineral dust in the global atmosphere: Impact of African dust on the environment of the southeastern United States, P. Natl. Acad. Sci. USA, 96, 3396–3403,, 1999. 

Pu, B. and Ginoux, P.: How reliable are CMIP5 models in simulating dust optical depth?, Atmos. Chem. Phys., 18, 12491–12510,, 2018. 

Pu, B., Ginoux, P., Guo, H., Hsu, N. C., Kimball, J., Marticorena, B., Malyshev, S., Naik, V., O'Neill, N. T., Pérez García-Pando, C., Paireau, J., Prospero, J. M., Shevliakova, E., and Zhao, M.: Retrieving the global distribution of the threshold of wind erosion from satellite data and implementing it into the Geophysical Fluid Dynamics Laboratory land–atmosphere model (GFDL AM4.0/LM4.0), Atmos. Chem. Phys., 20, 55–81,, 2020. 

Rahimi, S. R., Wu, C., Liu, X., and Brown, H.: Exploring a Variable-Resolution Approach for Simulating Regional Climate Over the Tibetan Plateau Using VR-CESM, J. Geophys. Res.-Atmos., 124, 4490–4513,, 2019. 

Ralaiarisoa, V., Dupont, P., Moctar, A. O. E., Naaim-Bouvet, F., Oger, L., and Valance, A.: Particle impact on a cohesive granular media, Phys. Rev. E, 105, 054902,, 2022. 

Raupach, M. R.: Drag and drag partition on rough surfaces, Bound.-Lay. Meteorol., 60, 375–395,, 1992. 

Raupach, M. R., Gillette, D. A., and Leys, J. F.: The effect of roughness elements on wind erosion threshold, J. Geophys. Res.-Atmos., 98, 3023–3029,, 1993. 

Reichle, R. H., Draper, C. S., Liu, Q., Girotto, M., Mahanama, S. P. P., Koster, R. D., and Lannoy, G. J. M. D.: Assessment of MERRA-2 Land Surface Hydrology Estimates, J. Climate, 30, 2937–2960,, 2017. 

Ridley, D. A., Heald, C. L., and Ford, B.: North African dust export and deposition: A satellite and model perspective, J. Geophys. Res.-Atmos., 117, D2,, 2012. 

Ridley, D. A., Heald, C. L., Pierce, J. R., and Evans, M. J.: Toward resolution-independent dust emissions in global models: Impacts on the seasonal and spatial distribution of dust, Geophys. Res. Lett., 40, 2873–2877,, 2013. 

Ridley, D. A., Heald, C. L., and Prospero, J. M.: What controls the recent changes in African mineral dust aerosol across the Atlantic?, Atmos. Chem. Phys., 14, 5735–5747,, 2014. 

Rosenfeld, D., Rudich, Y., and Lahav, R.: Desert dust suppressing precipitation: A possible desertification feedback loop, P. Natl. Acad. Sci. USA, 98, 5975–5980,, 2001. 

Roujean, J.-L., Tanré, D., Bréon, F.-M., and Deuzé, J.-L.: Retrieval of land surface parameters from airborne POLDER bidirectional reflectance distribution function during HAPEX-Sahel, J. Geophys. Res.-Atmos., 102, 11201–11218,, 1997. 

Shangguan, W., Dai, Y., Duan, Q., Liu, B., and Yuan, H.: A global soil data set for earth system modeling, J. Adv. Model. Earth Syst., 6, 249–263,, 2014. 

Shao, Y.: A model for mineral dust emission, J. Geophys. Res.-Atmos., 106, 20239–20254,, 2001. 

Shao, Y.: Simplification of a dust emission scheme and comparison with data, J. Geophys. Re.-Atmos., 109, D4,, 2004. 

Shao, Y.: Physics and Modelling of Wind Erosion, Springer Science & Business Media, 459 pp.,, 2008. 

Shao, Y. and Lu, H.: A simple expression for wind erosion threshold friction velocity, J. Geophys. Res.-Atmos., 105, 22437–22443,, 2000. 

Shao, Y. and Yang, Y.: A scheme for drag partition over rough surfaces, Atmos. Environ., 39, 7351–7361,, 2005. 

Shao, Y., Raupach, M. R., and Findlater, P. A.: Effect of saltation bombardment on the entrainment of dust by wind, J. Geophys. Res.-Atmos., 98, 12719–12726,, 1993. 

Shao, Y., Wyrwoll, K.-H., Chappell, A., Huang, J., Lin, Z., McTainsh, G. H., Mikami, M., Tanaka, T. Y., Wang, X., and Yoon, S.: Dust cycle: An emerging core theme in Earth system science, Aeolian Res., 2, 181–204,, 2011a. 

Shao, Y., Ishizuka, M., Mikami, M., and Leys, J. F.: Parameterization of size-resolved dust emission and validation with measurements, J. Geophys. Res.-Atmos., 116, D8,, 2011b. 

Shao, Y. P., Raupach, M. R., and Leys, J. F.: A model for predicting aeolian sand drift and dust entrainment on scales from paddock to region, Soil Res., 34, 309–342,, 1996. 

Sherman, D. J.: An equilibrium relationship for shear velocity and apparent roughness lenght in aeolian saltation, Geomorphology, 5, 419–431,, 1992. 

Shi, Y. and Liu, X.: Dust Radiative Effects on Climate by Glaciating Mixed-Phase Clouds, Geophys. Res. Lett., 46, 6128–6137,, 2019. 

Smith, M. B., Mahowald, N. M., Albani, S., Perry, A., Losno, R., Qu, Z., Marticorena, B., Ridley, D. A., and Heald, C. L.: Sensitivity of the interannual variability of mineral aerosol simulations to meteorological forcing dataset, Atmos. Chem. Phys., 17, 3253–3278,, 2017. 

Sokolik, I. and Golitsyn, G.: Investigation of optical and radiative properties of atmospheric dust aerosols, Atmos. Environ. Part A., General Topics, 27, 2509–2517,, 1993. 

Sokolik, I. N. and Toon, O. B.: Direct radiative forcing by anthropogenic airborne mineral aerosols, Nature, 381, 681–683,, 1996. 

Stull, R.: An Introduction to Boundary Layer Meteorology, Springer Dordrecht,, 670 pp., 1988. 

Tai, A. P. K., Ma, P. H. L., Chan, Y.-C., Chow, M.-K., Ridley, D. A., and Kok, J. F.: Impacts of climate and land cover variability and trends on springtime East Asian dust emission over 1982–2010: A modeling study, Atmos. Environ., 254, 118348,, 2021. 

Tegen, I. and Fung, I.: Contribution to the atmospheric mineral aerosol load from land surface modification, J. Geophys. Res.-Atmos., 100, 18707–18726,, 1995. 

Tegen, I., Harrison, S. P., Kohfeld, K., Prentice, I. C., Coe, M., and Heimann, M.: Impact of vegetation and preferential source areas on global dust aerosol: Results from a model study, J. Geophys. Res.-Atmos., 107, AAC14-1–AAC14-27,, 2002. 

Todd, M. C., Karam, D. B., Cavazos, C., Bouet, C., Heinold, B., Baldasano, J. M., Cautenet, G., Koren, I., Perez, C., Solmon, F., Tegen, I., Tulet, P., Washington, R., and Zakey, A.: Quantifying uncertainty in estimates of mineral dust flux: An intercomparison of model performance over the Bodélé Depression, northern Chad, J. Geophys. Res.-Atmos., 113, D24,, 2008. 

UCLA: The DustCOMM data set: Dust Constraints from joint Observational-Modelling-Experimental analysis, UCLA [data set], (last access: 30 March 2023), 2023. 

Van der Does, M., Knippertz, P., Zschenderlein, P., Harrison, R. G., and Stuut, J.-B. W.: The mysterious long-range transport of giant mineral dust particles, Sci. Adv., 4, eaau2768,, 2018. 

Webb, N. P., Galloza, M. S., Zobeck, T. M., and Herrick, J. E.: Threshold wind velocity dynamics as a driver of aeolian sediment mass flux, Aeolian Res., 20, 45–58,, 2016. 

Webb, N. P., Chappell, A., LeGrand, S. L., Ziegler, N. P., and Edwards, B. L.: A note on the use of drag partition in aeolian transport models, Aeolian Res., 42, 100560,, 2020. 

Wu, C., Lin, Z., He, J., Zhang, M., Liu, X., Zhang, R., and Brown, H.: A process-oriented evaluation of dust emission parameterizations in CESM: Simulation of a typical severe dust storm in East Asia, J. Adv. Model. Earth Syst., 8, 1432–1452,, 2016. 

Wu, C., Lin, Z., and Liu, X.: The global dust cycle and uncertainty in CMIP5 (Coupled Model Intercomparison Project phase 5) models, Atmos. Chem. Phys., 20, 10401–10425,, 2020. 

Wu, M., Liu, X., Yang, K., Luo, T., Wang, Z., Wu, C., Zhang, K., Yu, H., and Darmenov, A.: Modeling Dust in East Asia by CESM and Sources of Biases, J. Geophys. Res.-Atmos., 124, 8043–8064,, 2019. 

Xi, X. and Sokolik, I. N.: Seasonal dynamics of threshold friction velocity and dust emission in Central Asia, J. Geophys. Res.-Atmos., 120, 1536–1564,, 2015. 

Yu, Y., Kalashnikova, O. V., Garay, M. J., Lee, H., and Notaro, M.: Identification and Characterization of Dust Source Regions Across North Africa and the Middle East Using MISR Satellite Observations, Geophys. Res. Lett., 45, 6690–6701,, 2018. 

Yu, Y., Kalashnikova, O. V., Garay, M. J., Lee, H., Choi, M., Okin, G. S., Yorks, J. E., Campbell, J. R., and Marquis, J.: A global analysis of diurnal variability in dust and dust mixture using CATS observations, Atmos. Chem. Phys., 21, 1427–1447,, 2021. 

Zender, C. S., Bian, H., and Newman, D.: Mineral Dust Entrainment and Deposition (DEAD) model: Description and 1990s dust climatology, J. Geophys. Res.-Atmos., 108, D14,, 2003a. 

Zender, C. S., Newman, D., and Torres, O.: Spatial heterogeneity in aeolian erodibility: Uniform, topographic, geomorphic, and hydrologic hypotheses, J. Geophys. Res.-Atmos., 108, D17,, 2003b. 

Zender, C. S., Miller, R. L. R. L., and Tegen, I.: Quantifying mineral dust mass budgets:Terminology, constraints, and current estimates, Eos, Transactions American Geophysical Union, 85, 509–512,, 2004. 

Zhang, K., Zhao, C., Wan, H., Qian, Y., Easter, R. C., Ghan, S. J., Sakaguchi, K., and Liu, X.: Quantifying the impact of sub-grid surface wind variability on sea salt and dust emissions in CAM5, Geosci. Model Dev., 9, 607–632,, 2016.  

Zhang, L., Gong, S., Padro, J., and Barrie, L.: A size-segregated particle dry deposition scheme for an atmospheric aerosol module, Atmos. Environ., 35, 549–560,, 2001. 

Zhao, A., Ryder, C. L., and Wilcox, L. J.: How well do the CMIP6 models simulate dust aerosols?, Atmos. Chem. Phys., 22, 2095–2119,, 2022. 

Short summary
Desert dust modeling is important for understanding climate change, as dust regulates the atmosphere's greenhouse effect and radiation. This study formulates and proposes a more physical and realistic desert dust emission scheme for global and regional climate models. By considering more aeolian processes in our emission scheme, our simulations match better against dust observations than existing schemes. We believe this work is vital in improving dust representation in climate models.
Final-revised paper