# 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, et al.

► To cite this version:

Danny M. Leung, Jasper F. Kok, Longlei Li, Gregory S. Okin, Catherine Prigent, et al. A new process-based and scale-aware desert dust emission scheme for global climate models - Part I: Description and evaluation against inverse modeling emissions. Atmospheric Chemistry and Physics, 2023, 23, pp.6487-6523. 10.5194/acp-23-6487-2023. insu-04195512

HAL Id: insu-04195512

https://insu.hal.science/insu-04195512

Submitted on 4 Sep 2023

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
<figure><img src="image_2.png" /></figure>

<figure><img src="image_3.png" /></figure>

# 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<fnref n="1" />, Jasper F. Kok<fnref n="1" />, Longlei Li<fnref n="2" />, Gregory S. Okin<fnref n="3" />, Catherine Prigent<fnref n="4" />, Martina Klose<fnref n="5" />,
Carlos Pérez García-Pando⁶,<fnref n="7" />, Laurent Menut<fnref n="8" />, Natalie M. Mahowald<fnref n="2" />, David M. Lawrence⁹, and
Marcelo Chamecki<fnref n="1" />

**Correspondence:** Danny M. Leung (dannymleung@ucla.edu)

Received: 17 October 2022 – Discussion started: 24 October 2022

Revised: 7 February 2023 – Accepted: 30 March 2023 – Published: 14 June 2023

**Abstract.** 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.

<footnote n="1">¹Department of Atmospheric and Oceanic Sciences, University of California – Los Angeles, Los Angeles, California, USA</footnote>

<footnote n="2">²Department of Earth and Atmospheric Sciences, Cornell University, Ithaca, New York, USA</footnote>

<footnote n="3">³Department of Geography, University of California – Los Angeles, Los Angeles, California, USA</footnote>

<footnote n="4">⁴Sorbonne Université, Observatoire de Paris, Université PSL, CNRS, LERMA, Paris, France</footnote>

<footnote n="5">⁵Department Troposphere Research, Institute of Meteorology and Climate Research (IMK-TRO), Karlsruhe Institute of Technology (KIT), Karlsruhe, Germany
⁶Barcelona Supercomputing Center (BSC), Barcelona, Spain</footnote>

<footnote n="7">⁷Catalan Institution for Research and Advanced Studies (ICREA), Barcelona, Spain</footnote>

<footnote n="8">⁸Laboratoire de Météorologie Dynamique, École Polytechnique, Institut Polytechnique de Paris, Ecole Normale Supérieure, Sorbonne Université, CNRS, Palaiseau, France
⁹Climate and Global Dynamics Laboratory, National Center for Atmospheric Research, Boulder, Colorado, USA</footnote>
# 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 (Raubach 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., Darmen-ova 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 re-solve these issues, we propose a new approach that combines the drag partition effects of both elements, leveraging satel-lite 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 re-solve 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 fluctu-ate around the high emission threshold, as models do not ac-count 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 de-viation of the instantaneous wind to the boundary-layer dy-namical variables (Comola et al., 2019b). In this study, we will account for turbulent dust emissions by following Co-mola 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 mod-eling 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 emis-sions 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 reso-lutions, 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 oc-curs for soil moisture for instance, with its maxima smoothed out leading to an overestimation of dust emissions. Thus, al-though some GCMs and ESMs recently implemented more physical schemes (Zhao et al., 2022), their inability to re-solve 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 $u_{*it}/u_{*ft}$ 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 $D_p$ on both thresholds. These two thresholds are mainly governed by soil particle diameter $D_p$, air density $\rho_a$, and soil moisture $w$ (Greeley et al., 1997; Shao and Lu, 2000). Although there are multiple data sources of globally gridded products for $\rho_a$ and $w$, there are relatively few efforts on obtaining globally gridded $D_p$, since there are no methods for satellites to observe and derive surface $D_p$ 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 $\bar{D}_p$ 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):

$$u_{*ft} = u_{*ft0}(D_p, \rho_a)f_m(w), \quad (1)$$

where $u_{*ft0}$ is the “dry” fluid threshold friction velocity (in m s⁻¹) on a smooth and bare surface as a function of air density $\rho_a$(long, lat, t) (longitude, latitude, time) and $D_p$(long, lat), which in this study will be the median diameter $\bar{D}_p$ of a polydispersed, mixed soil PSD. $f_m$ is the correction factor for the presence of soil moisture $w$(long, lat, t); $f_m \ge 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 tun-
nel measurements (e.g., Greeley and Iversen, 1985), they obtained the following equation with fitting parameters:

$$u_{*ft0} = \sqrt{A(\rho_p g D_p + \gamma / D_p)} \rho_a^{-0.5}, \quad (2)$$

where $g = 9.81 \text{ m s}^{-2}$ is the gravitational acceleration, $\rho_p = 2650 \text{ kg m}^{-3}$ is the typical soil particle density, $\rho_a$ is in kg m$^{-3}$, and $A = 0.0123$ and $\gamma = 1.65 \times 10^{-4} \text{ kg s}^{-2}$ are empirical constants accounting for the aerodynamic forces and interparticle forces, respectively. Assuming an air density $\rho_a = 1.225 \text{ kg m}^{-3}$, Eq. (2) will yield the smallest $u_{*ft0}$ of $0.21 \text{ m s}^{-1}$ at $D_p = 80 \text{ µ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 $Re_p$. I&W82 calculates $u_{*ft0} = u_{*ft0}(\rho_p, \rho_a, g, D_p, Re_p)$ in a similar form to S&L00 (see the detailed solution in I&W82 or Kok et al., 2012), with

$$Re_p = u_{*ft0} D_p / \nu, \quad (3)$$

where $\nu$ is the kinematic viscosity of air. Since $u_{*ft0}$ is a function of $Re_p$, 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 $\rho_a$ and $D_p$ in calculation.

Regardless of whether S&L00 or I&W82 is used, different models make different assumptions for soil particle sizes $D_p$, 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 $D_p = 75 \text{ µm}$ for the threshold calculation following Zender et al. (2003a) (Olesen et al., 2013; also see the latest version of CLM5.0 technical note on https://escomp.github.io/ctsm-docs/versions/release-clm5.0/html/tech note/index.html, last access: 3 October 2022). Thus $u_{*ft0}$ becomes solely a function of $\rho_a$ in CESM (e.g., $u_{*ft0} = 0.21 \text{ m s}^{-1}$ for $\rho_a = 1.225 \text{ kg m}^{-3}$ at $D_p = 75 \text{ µm}$).

Most models parameterize the effect of soil moisture $f_m$ on $u_{*ft}$ following Fécan et al. (1999):

$$f_m = \sqrt{1 + 1.21(w - w_t)^{0.68}} \quad \text{for } w > w_t \quad (4a)$$

$$\begin{aligned} w_t &= a(17 f_{\text{clay}} + 14 f_{\text{clay}}^2) \\ &= a(0.17(\%\text{clay}) + 0.0014(\%\text{clay})^2), \end{aligned} \quad (4b)$$

where $w(\text{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); $w_t(\text{long, lat})$ is the threshold gravimetric water content above which $u_{*ft}$ increases; $f_{\text{clay}}(\text{long, lat})$ is the

fraction of clay content in the topmost layer of soil between zero and one; $\%\text{clay} = 100 f_{\text{clay}}$ 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). $w_t$ 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):

$$u_{*it} = B_{it} u_{*ft0}, \quad (5)$$

where $B_{it} = 0.82$ is approximately constant with soil properties and particle size (Bagnold, 1937; Kok et al., 2012). Eqs. (2)–(5) imply that $u_{*ft} \ge u_{*ft0} > u_{*it}$ 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 $f_m$ 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, $F_{\text{eff}}$, 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)

$$u_{*ft} = u_{*ft0}(D_p, \rho_a) f_m(w) / F_{\text{eff}}, \quad (6a)$$

where $F_{\text{eff}} < 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 $F_{eff}$ to $u_*$:

$$u_{*s} = u_* F_{eff}, \qquad (6b)$$

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 $z_0$ to quantify roughness. Because of the large differences in the length scales between mountain-s/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 $z_{0m}$ mainly represents roughness due to large-scale orography, forests, and/or urbanization (with sizes of $10^{-3}$ to $10^3$ m) with values ranging from ~1 cm to 1 m (Menut et al., 2013). Second, the aeolian roughness length $z_{0a}$ 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). $z_{0a}$ 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 $z_{0s}$ quantifies the roughness of a bed of fine soil particles in the absence of roughness elements. $z_{0s}$ characterizes the roughness of mobile, erodible soil particles over an exposed surface. $z_{0s}$ is directly related to the particle diameter $D_p$ by (Nikuradse, 1950; Sherman, 1992; Pierre et al., 2014b)

$$z_{0s} = 2D_p/30. \qquad (7)$$

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 $z_{0s}$. 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 $z_{0a}$. By arguing that the two wind speeds must be equal at the IBL height $\delta$, they derived $F_{eff}$ as a function of $z_{0a}$ and $z_{0s}$:

$$F_{eff} = 1 - \frac{\ln\left(\frac{z_{0a}}{z_{0s}}\right)}{\ln\left(\frac{\delta}{z_{0s}}\right)}. \qquad (8)$$

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 $\lambda$ to quantify roughness. $\lambda$ is defined as the total frontal area of roughness elements divided by the area of land $A: \lambda = \frac{n\lambda}{A}$, where $n$ and $\lambda$ are the obstacle height and width, and $\lambda$ 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 $\tau'_s$ is related to the total drag force $\tau$, given $\lambda$, the roughness-element basal area-to-frontal area ratio $\sigma$, and the ratio of the roughness element-to-surface drag coefficient $\beta$:

$$\frac{\tau'_s}{\tau} = \frac{1}{(1 - m\sigma\lambda)(1 + m\beta\lambda)}, \qquad (9a)$$

where $m$ is a geometric parameter to account for the spatial variability of $\tau'_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 $\sigma$, $m$, and $\beta$ (Darmenova et al., 2009; Xi and Sokolik, 2015). $\lambda$, however, is related to the abundance of obstacles and is thus spatially variable, and thus far there are no globally gridded datasets of $\lambda$ available. Most studies thus related grid-scale $\lambda$ to other grid-scale properties; for instance, Shao et al. (1996) linked $\lambda$ to the vegetation cover fraction $f_v$ using in situ observations:

$$\lambda = c_{\lambda} \ln(1 - f_v), \qquad (9b)$$

where $c_\lambda$ is a proportionality constant. Gridded $\lambda$ 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 $F_{eff}$ and $\lambda$ 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 $\lambda$ due to rocks. Similarly, vegetation roughness is historically mostly represented by $\lambda$ rather than $z_{0a}$, so there is no globally gridded $z_{0a}$ 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 $z_{0a}$ by relating plant $\lambda$ with $z_{0a}$ (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 $\lambda$ with $z_{0a}$, 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 $\lambda$ to $z_{0a}$ 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 $z_{0a}$ and $\lambda$ from both schools of drag partition parameterizations, quantifying the drag partitions of rocks and plants into one hybrid drag partition factor $F_{eff}$.

## 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:

$$F_d = \text{STC}_{\text{MB}} \varphi f_{\text{bare}} \frac{\rho_a}{g} u_{*s}^3 \left(1 - \frac{u_{*t}^2}{u_{*s}^2}\right) \left(1 + \frac{u_{*t}}{u_{*s}}\right) \quad (10)$$

where $F_d$ is the emission flux (in kg m$^2$ s$^{-1}$); $C_{MB}$ is a proportionality constant for bridging the gap between local-scale and large-scale dust fluxes; $\varphi = 10^{13.4 f_{\text{clay}}^{-6}$ is the sandblasting efficiency, the vertical dust emission flux produced per unit of horizontal saltation flux as a function of soil clay fraction $f_{\text{clay}}$; $u_{*t}$ is the dust emission threshold (in m s$^{-1}$; Z03 used $u_{*ft}$ as $u_{*t}$); and $f_{\text{bare}}$ 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 $f_{\text{bare}}$, Mahowald et al. (2006) used a simple parameterization in which $f_{\text{bare}}$ is a pure function of LAI (neglecting the effects of other objects such as snow, rocks, and buildings):

$$f_{\text{bare}} = 1 - \text{LAI}/\text{LAI}_{\text{thr}} \quad \text{for LAI} \le \text{LAI}_{\text{thr}} \quad (11a)$$

$$f_{\text{bare}} = 0 \quad \text{for LAI} > \text{LAI}_{\text{thr}}. \quad (11b)$$

While Mahowald et al. (2006) took $LAI_{thr} = 0.3$, we take $LAI_{thr} = 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 $LAI_{thr}$, emission ($F_d \propto f_{\text{bare}}$) drops more rapidly with LAI such that the drag partition effect is also incorporated in the $f_{\text{bare}}$ term. However, since we are considering $F_{eff}$ in this study, we can set a more realistic $LAI_{thr}$ value such that $f_{\text{bare}}$ 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 $C_d$ in K14) can be physically parameterized using the standardized fluid threshold $u_{*st} = u_{*ft} \sqrt{\rho_a / \rho_{a0}}$, which is $u_{*ft}$ scaled to the standard air density of $\rho_{a0} = 1.22 \text{ kg m}^{-3}$:

$$C_d = C_{d0} \exp \left( -C_e \frac{u_{*st} - u_{*st0}}{u_{*st0}} \right), \quad (12)$$

where $C_d(\text{long, lat, t})$ is the time-varying dust emission coefficient or soil erodibility coefficient, $C_{d0} = (4.4 \pm 0.5) \times 10^{-5}$, $C_e = 2.0 \pm 0.3$, and $u_{*st0} = 0.16 \text{ m s}^{-1}$. Furthermore, K14 derived a new dust emission equation for $F_d$ (kg m$^{-2}$ s$^{-1}$):

$$F_d = C_{\text{tune}} C_d f_{\text{bare}} f_{\text{clay}} \frac{\rho_a (u_{*s}^2 - u_{*t}^2)}{u_{*st}} \left(\frac{u_{*s}}{u_{*t}}\right)^{\kappa} \quad (13a)$$

$$\kappa = C_{\kappa} \frac{(u_{*st} - u_{*st0})}{u_{*st0}}, \quad (13b)$$

where $C_{\kappa} = 2.7 \pm 1.0$, $C_{\text{tune}} = 0.05$ is the proportionality constant, $f_{\text{bare}}$ is modeled by Eq. (11), $u_{*s} = u_{*t}$ in K14, and $u_{*t}$ is again the emission threshold (K14 assumed for simplicity that $u_{*t} = u_{*it} = u_{*ft}$). $\kappa$ is the fragmentation exponent which quantifies the sensitivity of $F_d$ to $u_{*s}$. Here we limit the value of $\kappa$ to 3 in order to prevent excessive sensitivity of the model to wind speeds, which can be problematic
around topography. From Eq. (12), $C_d$ 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^\circ \times 0.625^\circ$ 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.

# 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 $\bar{D}_p$ of the soil PSD. Thus, we ideally need a global gridded map of $\bar{D}_p$ 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 $\bar{D}_p$ 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 $\bar{D}_p$, we use a compilation of literature measurements to explore and construct relationships between $\bar{D}_p$ 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, $\bar{D}_p$ might depend on soil texture $f_a$ 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 $f_a$, and other soil properties to investigate their statistical relations and infer a new global distribution of $\bar{D}_p$. 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 ($CaCO_3$; %), 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 $\bar{D}_p$ by calculating the 50th percentile of the dry-soil PSD.
<table><caption>**Table 1.** Input meteorological and land-surface variables employed for running the standalone dust emission model in this study.</caption><thead><tr><th>Variable</th><th>Meteorological parameter (SI unit)</th></tr></thead><tbody><tr><td>u<sub>*</sub></td><td>MERRA-2 friction velocity (m s<sup>-1</sup>)</td></tr><tr><td>θ</td><td>MERRA-2 volumetric soil moisture (m<sup>3</sup> water m<sup>-3</sup> soil)</td></tr><tr><td>LAI</td><td>MERRA-2 leaf area index (m<sup>2</sup> leaf m<sup>-2</sup> land)</td></tr><tr><td>ρ<sub>a</sub></td><td>MERRA-2 air density (kg m<sup>-3</sup>)</td></tr><tr><td>z<sub>i</sub></td><td>MERRA-2 planetary-boundary-layer height (m)</td></tr><tr><td>T</td><td>MERRA-2 air temperature (K)</td></tr><tr><td>H</td><td>MERRA-2 sensible heat flux over land (W m<sup>-2</sup>)</td></tr><tr><td>φ</td><td>MERRA-2 porosity (dimensionless)</td></tr><tr><td>f<sub>clay</sub></td><td>SoilGrids clay fraction (fraction)</td></tr><tr><td>f<sub>silt</sub></td><td>SoilGrids silt fraction (fraction)</td></tr><tr><td>z<sub>0a</sub></td><td>Prigent et al. (2005) aeolian roughness length (m)</td></tr><tr><td>A<sub>r</sub></td><td>Rock and bare soil land cover derived from the European Space Agency land cover dataset (fraction)</td></tr><tr><td>A<sub>v</sub></td><td>Vegetation land cover derived from the European Space Agency land cover dataset (fraction)</td></tr><tr><td>S</td><td>Ginoux et al. (2001) or Zender et al. (2003b) source function</td></tr></tbody></table>

### 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 $\overline{D}_p$ 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 LAI<sub>thr</sub>, 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 $\overline{D}_p$ and the soil properties (see Figs. 1b and S2). Figure 1b shows a scatterplot of $\overline{D}_p$ versus the sum of the soil component that produces substantial cohesion, namely the silt and clay fractions ($f_{\text{silt+clay}} = f_{\text{silt}} + f_{\text{clay}}$). The data exhibit distinctly different trends for nonarid versus arid soils: for nonarid soils, $\overline{D}_p$ increases from 100 µm to greater than 1000 µm with $f_{\text{silt+clay}}$ (regression $p$ value = 7.3 × 10<sup>-6</sup>), likely due to increasing cohesion with increasing clay and silt content. In contrast, $\overline{D}_p$ for arid soils shows a small and statistically insignificant increasing trend with $f_{\text{silt+clay}}$ ($p$ value = 0.77) with a smaller $\overline{D}_p$ variability (50–250 µm). This flat trend indicates that $f_{\text{silt+clay}}$ does not effectively explain the median diameter of aggregated soil particles in arid regions. We examined the relationships of $\overline{D}_p$ with the individual fractions of sand, silt, and clay, as well as with other soil properties including SOC, pH, and $\text{CaCO}_3$ (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 $\overline{D}_p$ over the arid regions across different soil textures. We thus use a constant $\overline{D}_{p0}$ as an approximation for arid regions. From Fig. 1b, we summarize the relationship between $\overline{D}_p$

and $f_{\text{silt+clay}}$ as

$$ \overline{D}_p = \begin{cases} \Psi_0 + \Psi_1 f_{\text{silt+clay}} & \text{, for LAI} > \text{LAI}_{\text{thr}} \\ \overline{D}_{p0} & \text{, for LAI} \le \text{LAI}_{\text{thr}} \end{cases}, \quad (14) $$

where $\Psi_0 = 7.8 \pm 3.7$ µm, $\Psi_1 = 124 \pm 36$ µm, $\overline{D}_{p0} = 127 \pm 47$ µm, and LAI<sub>thr</sub> = 1 as specified in Eq. (11). This empirical formula suggests that some models' assumptions of the relationship between $\overline{D}_p$ and soil texture were inaccurate (e.g., Table 2 of Laurent et al. (2008) assumed $\overline{D}_p$ decreases with $f_{\text{silt+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_{*f0}$ from 0.204 to 0.216 m s<sup>-1</sup> (given $\rho_a = 1.225$ kg m<sup>-3</sup>), which slightly decreases global dust emissions by 18% (see Sect. 4.1). The uncertainty in $\overline{D}_{p0} = 127 \pm 47$ µm translates to an uncertainty of $u_{*f0}$ between 0.204 to 0.234 m s<sup>-1</sup>.

We then project our derived relation between $\overline{D}_p$ and $f_{\text{silt+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 $\overline{D}_p$ map. Following Eq. (14), the arid and semiarid regions are set to have a $\overline{D}_{p0}$ of 127 µm, whereas for nonarid
<figure><img src="image_4.png" /><figcaption>**Figure 1.** Constructing a global map of the median diameter $D_p$ 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 $D_p$ versus silt + clay fraction ($f_{silt+clay}$). (c) The predicted soil median diameter $\bar{D}_p$ map (in µm) derived by projecting our derived $\bar{D}_p-f_{silt+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 $\bar{D}_p$ at those sites. (d) The predicted $\bar{D}_p$ using Eq. (14) versus measured $\bar{D}_p$ from past studies.</figcaption></figure>

regions, $\bar{D}_p$ increases with $f_{silt+clay}$. Our derived $\bar{D}_p$ values are largely consistent with the site $\bar{D}_p$ 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 $R^2$ 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 $\bar{D}_p-f_{silt+clay}$ relation over nonarid regions (blue points). Figure 1d shows that Eq. (14) gives satisfactory agreement in predicting global $\bar{D}_p$, but dust emission modeling will depend exclusively on the predicted $\bar{D}_p$ 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 $\bar{D}_p$ 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 $\bar{D}_{p0} = 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 nonerodable 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 $f_{\text{eff,r}}$ using $z_{0a}$ since global $z_{0a}$ observations are available, and we calculate the vegetation drag partition $f_{\text{eff,v}}$ using vegetation cover which is a proxy of $\lambda$ (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 $F_{\text{eff}}$ (Sect. 3.2.3).

### 3.2.1 Drag partition due to roughness of rocks

In this study, we use the aeolian roughness length $z_{0a}$ to quantify the drag partition effect due to rocks. Whereas the smooth $z_{0s}$ and the aerodynamic momentum $z_{0m}$ can be derived from pre-existing datasets, it is more challenging to quantify the aeolian $z_{0a}$. 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, $f_{\text{eff,r}}$ (also see Eq. 8):

$$f_{\text{eff,r}} = 1 - \frac{\ln\left(\frac{z_{0a}}{z_{0s}}\right)}{\ln\left[b_1\left(\frac{X}{z_{0s}}\right)^{b_2}\right]}, \quad (15)$$

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., $\frac{\delta}{z_{0s}} = b_1 \left(\frac{X}{z_{0s}}\right)^{b_2}$ following Marticorena and Bergametti (1995), and $b_1 = 0.7$ and $b_2 = 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 10$^0$–10$^1$ 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 $\overline{D}_p$ in Sect. 3.1 using Eq. (7). When nonerodible roughness elements are abundant over a surface, z0a ≫ z0s and $f_{\text{eff,r}} \ll 1$, 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 $f_{\text{eff,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 $f_{\text{eff,r}}$ over arid regions, which is dominated by the spatial pattern of z0a in Fig. 2a given $f_{\text{eff,r}}$ is governed purely by z0a.

### 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 un-vegetated 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

$$u_{*s}(x/h) = u_* \left[ f_0 + (1 - f_0) \left(1 - e^{-x/h/\epsilon}\right) \right], \quad (16a)$$

where $x/h$ is the dimensionless downwind distance from an obstacle normalized by vegetation height h (m), $f_0 = \frac{u_{*s}}{u_*}|_{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
<figure><img src="image_5.png" /><figcaption>**Figure 2.** Global roughness length and rock drag partition factor maps at a horizontal resolution of 0.25° × 0.25°. (a) Global mesoscale aeolian roughness length $z_{0a}$ (in centimeters) derived by Prigent et al. (2005). (b) Global static rock drag partition factor $f_{eff,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.</figcaption></figure>

of distance $x$ is

$$f_{\text{local}}\left(\frac{x}{h}\right) = \frac{u_{*s}}{u_{*}} = f_0 + (1 - f_0)(1 - e^{-\frac{x/h}{c}}). \quad (16b)$$

Note that in the limit of $x/h \to \infty$, $u_{*s} \to u_{*}$. O08 used measurements from Bradley and Mulhearn (1983) and fitted $f_0 = 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 $f_{eff,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 $f_{eff,v}$ (see Eq. 20a) instead of a locally varying $f_{local}$ (Eq. 16b). Therefore, Okin employed a probability distribution function as a function of distance $x/h$ to indicate the importance (or weight) of $f_{local}$ at any $x/h$ to the averaging of $f_{eff,v}$ (McGlynn and Okin, 2006; Okin, 2008). The PDF is an exponential decay such that the weight of $f_{local}$ decreases with distance $x/h$, so $f_{local}$ at the immediate lee of the obstacle (which is smaller and close to $f_0$) has more weight than the $f_{local}$ farther away (which is larger and tends to 1). From McGlynn and Okin (2006), the PDF is a function of normalized distance $x/h$:

$$P_d(x/h) = \frac{1}{K}e^{-\frac{x/h}{K}} \quad (17a)$$

$$K = L/h, \quad (17b)$$

where $L$ (m) is the mean gap length between obstacles (plants), which is conceptually related to $f_v$; and $K$ is the normalized gap length, which is the gap length $L$ scaled by the plant height $h$. Physically, $P_d$ 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 $f_0$ and $c$ are assumed to be invariant to surface conditions and desert biome. $K$ can be expressed as a function of $f_v$ using some simple assumptions. First, O08 argued that the vegetation cover fraction is simply $f_v \equiv \frac{W}{L+W}$, where $L$ is the mean gap length and $W$ is the mean width of the plants within that vegetated area. Rearranging gives

$$L = W\left(\frac{1}{f_v} - 1\right). \quad (18a)$$

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

$$K = \frac{L}{h} = 2\left(\frac{1}{f_v} - 1\right), \quad (18b)$$

and thus we related $K$ to $f_v$. $f_v$ 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 $f_v$, 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 $f_v = 1 - f_{bare} = LAI/LAI_{thr}$. Equation (18b) becomes

$$K = \frac{L}{h} = 2\left(\frac{1}{f_v} - 1\right) = 2\left(\frac{1}{LAI/LAI_{thr}} - 1\right), \quad (18c)$$

where we assume $LAI_{thr} = 1$ (in Eq. 11). The assumption of $f_v \sim 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 $f_v$ 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 $f_{eff,v}$. O08 calculates the reduced dust emission flux $F_{red}$ (kg m⁻² s⁻¹) by locally integrating the emission $F_d$ following the spatially varying $u_{*s}$ over the normalized distance $x/h$:

$$F_{red} = \int_{x/h=0}^{\infty} P_d\left(\frac{x}{h}\right) F_d\left[u_{*s}\left(\frac{x}{h}\right)\right] d\left(\frac{x}{h}\right), \quad (19)$$

where $F_d$ (kg m⁻² s⁻¹) is the local emission as a function of $u_{*s}$ which is itself a function of $x/h$. In the integration, $F_d$ needs to be weighted by $P_d$ (which means $F_d$ 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 $F_d$ from zero to infinity gives a reduced emission flux $F_{red}$, which will be smaller than the emission flux without roughness elements, defined as $F_{bare} = \int_{x/h=0}^{\infty} P_d\left(\frac{x}{h}\right) F_d(u_*) d\left(\frac{x}{h}\right) = F_d(u_*)$, in which $F_d$ is a constant in space since $u_{*s} = 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 $f_{eff,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 $f_{eff,v}$ such that $F_d(u_{*s}, u_*) = F_{red}$. Quantifying $f_{eff,v}$ for O08 can be useful for comparisons against $f_{eff,v}$ from other schemes such as R93 and Klose et al. (2021). In addition, quantifying $f_{eff,v}$ for O08 makes it possible to generate a high-resolution, diagnostic $f_{eff,v}$ dataset for mechanistic models with different resolutions as a model input.

An approach of evaluating $f_{eff,v}$ from O08 was proposed by Pierre et al. (2014a). They obtained the expected value of the shear stress ratio $u_{*s}/u_*$ (SSR in Okin, 2008) between obstacles by evaluating the integral of $u_{*s}/u_*$ weighted by $P_d$, which represents the averaged $f_{local}$ (in Eq. 16b) across the vegetated area and perfectly fits our purposes for implementing $f_{eff,v}$ into LSMs:

$$\begin{aligned} f_{eff,v} &= \int_{x/h=0}^{\infty} P_d\left(\frac{x}{h}\right) \left[ \frac{u_{*s}\left(\frac{x}{h}\right)}{u_*} \right] d\left(\frac{x}{h}\right) \\ &= \int_{x/h=0}^{\infty} P_d\left(\frac{x}{h}\right) f_{local}\left(\frac{x}{h}\right) d\left(\frac{x}{h}\right). \end{aligned} \quad (20a)$$

Substituting Eq. (17) for $P_d$ into Eq. (20a) and analytically evaluating the integral gives a simple algebraic equation for

$f_{eff,v}$ (Pierre et al., 2014a), representing the aggregated vegetation drag partition effect at the grid level:

$$f_{eff,v} = \frac{K + f_{0c}}{K + c}. \quad (20b)$$

This elegant formula conveys a clear physical intuition: if the obstacle does not effectively dissipate momentum ($f_0 \to 1$), $f_{eff,v} \to 1$; if land is densely covered by vegetation (gap length $K \to 0$), $f_{eff,v} \to f_0(=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 $f_{eff,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 $f_{eff,v}$ for areas where LAI < 1. The LAI plot shows the most erodible regions on Earth.

### 3.2.3 Combining drag partition factors of rocks and vegetation

After obtaining both the static $f_{eff,r}$ map of rocks and the time-varying $f_{eff,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 $F_{eff}$ 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 (https://www.esa-landcover-cci.org/?q=node/164, 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
<figure><img src="image_6.png" /><figcaption>**Figure 3.** Vegetation drag partition factor $f_{eff,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 $f_{eff,v}$ derived using the Okin (2008) and Pierre et al. (2014) drag partition model. White areas indicate water body, ice/snow, or LAI > 1.</figcaption></figure>

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 $F_{eff}$ for a grid is thus defined by the summation of emissions, with emission $F_{d,r}$ over the rock regime with a fractional area of $A_r$, and emission $F_{d,v}$ over the vegetation regime with another fractional area $A_v$:

$$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}) \quad (21a)$$

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

$$F_{eff}^3 = A_r f_{eff,r}^3 + A_v f_{eff,v}^3 \quad (21b)$$

such that $F_{eff}$ 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 $F_{eff}$ (long, lat, t), given $f_{eff,v}$ (long, lat, t) and $A_r$, $A_v$, and $f_{eff,r}$ as functions of (long, lat). We then apply the obtained $F_{eff}$ 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 $F_{eff}$ using Eq. (21b). The regions with the highest $F_{eff}$ are the Bodélé Depresion, El Djouf, the Arabian Desert, and Taklamakan due to high $f_{eff,r}$. The Strzelecki–Sturt Stony deserts in Australia, the Kyzylkum, and Patagonia also have high $F_{eff}$ (~0.7) due to high $f_{eff,v}$. Regions with both high rock and vegetation roughness are located in parts of the Middle East and North America with low $F_{eff}$ values.

## 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
<figure><img src="image_7.png" /><figcaption>**Figure 4.** The 0.5° × 0.625° hybrid drag partition factor $F_{\text{eff}}$ 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 $F_{\text{eff}}$ by a combination of rock drag partition $f_{\text{eff,r}}$ and vegetation drag partition $f_{\text{eff,v}}$ for the year 2006.</figcaption></figure>

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 $\tilde{u}_s$ (we use the tilde to denote instantaneous quantities and take away the asterisk to denote winds at the saltation level of $z_{\text{sal}} \sim 0.1$ m using Eq. (S4a) instead of a velocity scale) needs to exceed the fluid threshold $u_{\text{ft}}$ (also defined at the saltation level) to initiate saltation, but it only needs to exceed a smaller impact threshold $u_{\text{it}}$ to sustain it (Kok et al., 2012; Martin and Kok, 2018; Comola et al., 2019b). When $\tilde{u}_s$ at a moment lies between both thresholds ($u_{\text{it}} < \tilde{u}_s < u_{\text{ft}}$), saltation is active if transport was more recently initiated ($\tilde{u}_s > u_{\text{ft}}$) and inactive if transport was more recently terminated ($\tilde{u}_s < u_{\text{it}}$). This process is known as hysteresis (Kok, 2010; Martin and Kok, 2018; Comola et al., 2019b). As a result, if $u_s$ (mean of $\tilde{u}_s$ within a model time step) is between $u_{\text{it}}$ and $u_{\text{ft}}$, 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.

How- ever, because in reality high-frequency winds can pass below $u_{\text{it}}$ and shut off dust emissions, using average $u_s$ 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 $u_{\text{it}}$ and the Monin–Obukhov length $L$. C19 showed that when a dust emission equation employs $u_{\text{it}}$ and accounts for the intermittency effect, it can successfully capture the magnitudes of small dust fluxes otherwise missed by models using $u_{\text{ft}}$ (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 $F_d$ is calculated using $u_{*it}$ instead of $u_{*ft}$. We update K14 (Eq. 13) with $u_{*t} = u_{*it}$ as the threshold, giving

$$F_d = C_{\text{tune}} C_d f_{\text{bare}} f_{\text{clay}} \frac{\rho_a (u_{*s}^2 - u_{*it}^2)}{u_{*it}} \left( \frac{u_{*s}}{u_{*it}} \right)^{\kappa} \quad (22a)$$
for $u_{*s} > u_{*it}$,

where $C_d$ is still a function of $u_{*st}$, and $u_{*st} = u_{*ft} \sqrt{\rho_a / \rho_{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 $\eta$, which is the fraction of time that saltation is active in a model time step (e.g., ~30 min). $\eta$ 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 $\eta$ as follows:

$$F_{d,\eta} = \eta F_d, \quad (22b)$$

where $\eta \in [0, 1]$. Note that C19 parameterizes $\eta$ using wind information at the typical saltation height of $z_{sal} = 0.1$ m instead of the velocity scales. $\eta$ is thus formulated as a function of the wind speeds $u_s$, $u_{it}$, and $u_{ft}$ at height $z_{sal}$ (see Sect. S3 Eqs. S3–S6) and the standard deviation $\sigma_{\tilde{u}_s}$ of the instantaneous $\tilde{u}_s$ (Eq. 23):

$$\eta = \eta(u_s, \sigma_{\tilde{u}_s}, u_{it}, u_{ft}). \quad (22c)$$

$\sigma_{\tilde{u}_s}$ is defined given that $\tilde{u}_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) $u_s$ and its standard deviation $\sigma_{\tilde{u}_s}$. From Eq. (22c), $\eta \to 1$ when $u_s \gg u_{ft}$ (active emission for the whole time step) and $\eta \to 0$ when $u_s \ll u_{it}$ (no emission for the time step). The further away $u_s$ is from the thresholds $u_{ft}$ and $u_{it}$, the smaller the probability of the instantaneous $\tilde{u}_s$ sweeping across the thresholds and the more dichotomous $\eta$ behaves (either zero or 1). If $u_s$ is very close to $u_{ft}$ or $u_{it}$, or is indeed between them ($u_{it} < \tilde{u}_s < u_{ft}$), the frequency of crossing the threshold is determined by the magnitude of the turbulent fluctuation $\sigma_{\tilde{u}_s}$. $\sigma_{\tilde{u}_s}$ is parameterized using the similarity theory (Panofsky et al., 1977):

$$\sigma_{\tilde{u}_s} = u_{*s} \left( 12 - 0.5 \frac{z_i}{L} \right)^{1/3} \quad \text{for } 12 - 0.5 \frac{z_i}{L} \ge 0, \quad (23)$$

where $L$ is the Obukhov length, and $z_i$ 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_*$, $\rho_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 $\sigma_{\tilde{u}_s}$ increase with shear ($u_{*s} > 0$) and buoyancy ($L < 0$). A larger $\sigma_{\tilde{u}_s}$ makes it easier for $\tilde{u}_s$ to sweep across $u_{it}$ and shut off dust emission, leading to $\eta < 1$. In a time step, if $u_s \gg u_{ft} + \sigma_{\tilde{u}_s}$, $\tilde{u}_s$ will be unlikely to sweep across $u_{it}$, and $\eta$ will approach 1. If $u_s$ is slightly larger than $u_{ft}$, the instantaneous $\tilde{u}_s$ will be likely to sweep across $u_{it}$, leading to $\eta < 1$. In the hysteresis regime ($u_{it} < u_s < u_{ft}$), $\eta$ will be around 0.3–0.7, since $\tilde{u}_s$ will sweep across both thresholds given $\sigma_{\tilde{u}_s}$, leading to a reduced emission flux (meanwhile, other parameterizations predict a zero emission flux since they use $u_{ft}$ only). When $u_s < u_{it}$, $\eta$ could also be greater than zero when $\sigma_{\tilde{u}_s}$ is large enough so that the instantaneous $\tilde{u}_s$ sweeps across $u_{it}$. 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 $z_i/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 $\eta$ to yield the dust emission with intermittency effect $F_{d,\eta}$ 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 $\overline{D}_{p0} = 127 \ \mu\text{m}$. Its spatial variability is purely a function of $\rho_a$ (Eq. 2). $u_{*it}$ is around ~0.16–0.24 m s⁻¹, 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 $\rho_a$ regions. Figure 5b shows $u_{*ft}$ (Eq. 1). It varies between 0.2 and 0.9 m s⁻¹. 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 $u_{*ft}/u_{*it}$, 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 $u_{*ft}/u_{*it} \sim 1/B_{it}$ have a narrow hysteresis regime ($u_{*it} < u_{*s} < u_{*ft}$) and smaller thresholds and thus tend to have more continuous dust emissions. Semi-arid and nonarid regions with larger $u_{*ft}/u_{*it}$ tend to have a wide hysteresis regime, and thus dust emissions will be more intermittent.

Figure 6a shows the 2006 annual mean intermittency effect over the Bodélé Depression as an example. Figure 6a
<figure><img src="image_8.png" /><figcaption>**Figure 5.** The 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 $\bar{D}_{p0}$ = 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 $f_m/0.81$, where $f_m$ 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.</figcaption></figure>

shows hourly mean $\eta$ as a function of the hourly mean $u_s$. It demonstrates the properties of $\eta$ discussed above: e.g., when $u_s > u_{ft}$, $\eta \to 1$ and when $u_s < u_{it}$, $\eta \to 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 $u_{it} < u_s < u_{ft}$, $\eta$ 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 $-z_i/L$. $-z_i/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 $-z_i/L$ and $u_s$, but the correlation is not necessarily strong, and the strongest buoyancy (dark red) often happens when $u_s$ (or shear $u_{*s}$) is moderate. The strongest buoyancy associates with moderate $\eta$ values of ~0.5 only, and for the highest $\eta$ values $-z_i/L$ is mildly unstable (light red). This means that the turbulent fluctuation is primarily governed by shear $u_{*s}$ instead of controlled by buoyancy $-z_i/L$, and the intermittency behavior is dictated by shear-driven instead of buoyancy-driven turbulence. Equation (23) could essentially be simplified into $\sigma_{\bar{u}_s} \approx 12^{1/3} u_{*s}$.

Figure 6b shows the global spatial distribution of the annual mean $\eta$ for the year 2006, averaged across time steps during which saltation and dust emissions are occurring over the grid (and thus $\eta$ during time steps when $F_d = 0$ is not counted). Most marginal sources have small $\eta < 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 $\eta$ which is also missed by other LSMs. Intermittency is thus critical in accounting for emissions over semiarid regions. Regions with high $\eta$ (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 $\eta$. However, regions with high $\eta$ 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 $\eta$ averaged over all time steps of 2006, including periods of no emissions over the grid. Since $\eta$ is close to 0 when there is no emission, Fig. S5 shows much smaller $\eta$ compared to Fig. 6b.

## 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
<figure><img src="image_9.png" /><figcaption>**Figure 6.** Simulated dust emission intermittency effect for the year 2006. (a) The hourly mean intermittency factor $\eta$ versus hour mean wind speed $u_s$ (m s$^{-1}$) at 0.1 m height over the Bodélé Depression. The color indicates $-z_i/L$, with red indicating $-z_i/L > 0$ (unstable), blue indicating $-z_i/L < 0$ (stable), and gray indicating $-z_i/L \sim 0$ (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 $\eta$, or equivalently the fraction of time within a time step that emission is active, averaged over times when emission is active ($F_d > 0$).</figcaption></figure>

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 (Soil-Grids 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 $\bar{D}_{p0} \sim 127$ µ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 PM$_{2.5}$ 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 $C_d$, and the intermittency $\eta$. 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 $\bar{D}_{p0} = 127$ µm (experiment II). As described in Sect. 3.1, a globally larger $\bar{D}_p$ 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 $\bar{D}_{p0}$ is relatively mild (of the largest order of $0.001 \text{ kg m}^{-2} \text{ 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 $\bar{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 $F_d$ 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 $\bar{D}_p$ 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 $F_{\text{eff}}$ (experiment III). The clear contrast between major and marginal sources is shown in Fig. 7e and f, which mirrors the spatial pattern of $F_{\text{eff}}$ 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, $F_{\text{eff}}$ 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 $\bar{D}_p$ 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 $F_{\text{eff}} < 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. $F_{\text{eff}}$ 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 $\eta$ (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 $u_{*ft}/u_{*it}$ bigger but the fragmentation exponent $\kappa$ (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.
ies 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. Fig-
ure 7g shows that although the changes over nonarid re-
gions are high in ratios (e.g., > 1000 times over Canada in
Fig. 7h), the emission redistributions in magnitude are still
small (e.g., ~10⁻⁵ kg m⁻² yr⁻¹ in Canada) compared to the
major sources (e.g., > 10⁻³ kg m⁻² yr⁻¹ in the Sahara) be-
cause the magnitudes in experiment III are too small over
nonarid regions. The unnormalized total emissions (Fig. S6d)
vastly increase to ~13 200 Tg yr⁻¹ when employing $u_{*it}$ in-
stead 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⁻¹, which
is only 11% smaller than experiment IV, because those re-
mote regions as such already have small emissions. All in
all, accounting for both the impact threshold and the inter-
mittency factor will increase the global total emission from
~2880 Tg yr⁻¹ in experiment III to ~11 700 Tg yr⁻¹ in ex-
periment V, which is about a 4-fold increase. Figure 7h shows
that C19 mainly increases marginal emissions; the overall ef-
fect of C19 is thus to move emissions from the hyper-arid
regions to semiarid regions. Summing up the absolute mag-
nitude of changes in Fig. 7g, out of 5000 Tg yr⁻¹ there is
3163 Tg yr⁻¹ of dust redistributed within the source regions,
indicating that the intermittency scheme induces a similar
magnitude of changes compared to employing $F_{\text{eff}}$. 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 $\bar{D}_p$ (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 (exper-
iment 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 north-
ern 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 ef-
fect. The figures show essentially a combination of drag par-
tition (Fig. 7e-f) and intermittency (Fig. 7g-h) effects. Ma-
jor 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⁻¹ is ~60% smaller than the K14 emis-
sion. 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⁻² yr⁻¹ (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
$F_{\text{eff}}$), and moderate soil moisture (and thus high dust emis-
sion coefficient $C_d$). 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.

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

The spatial resolution of a GCM strongly affects the bud-
get and spatiotemporal variability of dust emissions because
modeled emissions scale nonlinearly with input meteorolog-
ical 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 aeo-
lian roughness length $z_{0a}$, are datasets of higher horizontal
resolutions. Thus, GCMs need to regrid the datasets to the
model native grid resolution as model input. Since dust emis-
sion 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$ $\overline{u}_{*s}^n < \overline{u}_{*s}^{n_1}$ leads to
$F_d(\overline{u}_{*s}^n) < \overline{F}_d(\overline{u}_{*s}^n)$ (Ridley et al., 2013) (here we use the bar
as a spatial average from a finer to a coarser grid). This in-
equality 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 di-
rect area-weighted averaging of the emissions. Additionally,
the presence of thresholds in dust emission parameterizations
further intensifies the scale dependence, because spatially av-
eraging $u_{*s}$ might cause $\overline{u}_{*s} < \overline{u}_{*it}$ which leads to zero emis-
sion over a coarse grid box, whereas fine emissions could be
large than zero when $u_{*s} > u_{*it}$ in any fine grid. Dust emis-
sion 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 ad-
dress the problem of missing emissions due to the smoothing
of the subgrid wind maxima, a common approach is to em-
ploy a grid-by-grid Weibull distribution to the GCM winds to
represent the subgrid wind maxima and thus obtain the sub-
<figure><img src="image_10.png" /></figure>

**Figure 7.** The 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 $log_{10}$ scale.
grid 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⁻¹) 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⁻¹, 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 Chinese 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 log₁₀ 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⁻⁵ kg m⁻² yr⁻¹), 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_{*f}$ instead of $u_{*s}$ and thus the overestimation of $F_d$ 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.

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 $F_{d,f}$ by the coarse-resolution emission map $F_{d,c}$, we obtain a map of scaling factors $\tilde{K}_c$ to account for the differences in the spatial variability of dust emissions between high- and low-resolution simulations:

$$ \tilde{K}_c(\text{long, lat}) = F_{d,f}(\text{long, lat}) / F_{d,c}(\text{long, lat}). \quad (24) $$

We obtain a global map of correction factors $\tilde{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
<figure><img src="image_11.png" /><figcaption>**Figure 8.** The 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.</figcaption></figure>

seasonal variability of $\tilde{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 $\tilde{K}_c$. However, $\tilde{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 $\tilde{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 ($F_{d,0.5}$) to constrained 1.9° × 2.5° emissions ($F_{d,2}$), and Fig. 9b shows the ratio of constrained $F_{d,0.5}$ to constrained 0.9° × 1.25° emissions ($F_{d,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 $F_d$; 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 $F_d$. 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° ($\tilde{K}_{c,2}$) generally have bigger magnitudes than those in Fig. 9a for 0.9° × 1.25° ($\tilde{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 dif-
ferent spatial resolutions to the finest spatial resolution will move our dust emission scheme toward a scale-aware and grid-independent formulation.

## 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) × 10³ Tg yr⁻¹ for dust PM₂₀. The global total of ~4700 Tg yr⁻¹ is close to the 5000 Tg yr⁻¹ we adopted in this study for normalization, and we again normalized the DustCOMM global emissions to 5000 Tg yr⁻¹. 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
<figure><img src="image_12.png" /><figcaption>**Figure 9.** Gridded 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⁻¹ before calculating the ratios.</figcaption></figure>

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⁻¹ of global total emissions. We normalized their estimate to

match our global total of 5000 Tg yr⁻¹, yielding a high-latitude emission range of 200–250 Tg yr⁻¹. We take the average, which is 225 Tg yr⁻¹. 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 $R^2$ 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 $R^2$ 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 $R^2$ of 64%. The Z03–G simulation (Fig. 10c in violet) has a higher $R^2$ 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 $R^2$ 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).

# 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 $\bar{D}_p$ representation, we used a globally constant $\bar{D}_{p0} = 127 \pm 47$ µm for calculating dust emission thresholds. In theory, $\bar{D}_p$ 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 $\bar{D}_p$ over arid regions because (1) there were no statistically significant correlations between $\bar{D}_p$ 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 $\bar{D}_p$ measurements over arid regions found $\bar{D}_p$ 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 $\rho_a = 1.225$ kg m$^{-3}$). Our analysis of a compilation of $\bar{D}_p$ measurements also suggests that the spatial variability of $\bar{D}_p$ over deserts is relatively small compared to $\bar{D}_p$ over nonarid regions. Our results thus surprisingly suggest that the $\bar{D}_p$ parameterization can reasonably be much simplified by using a constant value over arid regions. We also compared our approach in deriving $\bar{D}_p$ against some other derived $\bar{D}_p$ 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 $z_{0a}$ of rocks and $z_{0a}$ 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 aero-lian roughness. We also note that ESA CCI vegetation cover fractions are annual values and do not exhibit seasonality. The seasonal variability of the simulated $F_{eff}$ 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}$, $\sigma_{\tilde{u}_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 $\eta$ 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
<figure><img src="image_13.png" /></figure>

**Figure 10.** Dust 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⁻² yr⁻¹) 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⁻¹) 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.

Navier–Stokes (RANS) mode. If climate models are resolv-
ing 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 an-
nual 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 sim-
ulations 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 climato-
logical 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 meteoro-

logical 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 fi-
nal simulation of our scheme outperforms other dust emis-
sion schemes employed in this paper (Fig. 10), some fea-
tures in the dust emissions map are unrealistic. First, for in-
stance, 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 northeast-
ern China (annual mean > 0.3 m s⁻¹ in MERRA-2) than over
northwestern China (annual mean ~ 0.2 m s⁻¹ in MERRA-
2) in the GCMs or input MERRA-2 fields (see Fig. S15
for the spatial distribution of MERRA-2 u* over northeast-
ern 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_*(\sim 1 \text{ 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 $\bar{D}_p$ as a critical parameter that determines the dust emission thresholds $u_{*\text{it}}$ (impact) and $u_{*\text{ft}}$ (fluid). We found that over the arid desert regions ($\text{LAI} < 1$), $\bar{D}_p$ can be approximated as a global constant of $127 \pm 47 \text{ µm}$, and over nonarid regions $\bar{D}_p$ increases linearly with silt and clay content. This finding indicates that past dust modeling approaches which parameterized $\bar{D}_p$ 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 $R^2 = 89\%$, meanwhile other schemes, such as Kok et al. (2014a, b) and Zender et al. (2003a, b) respectively yielded $R^2 = 65\%$ and $R^2 = 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 $\bar{D}_p$. 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 https://doi.org/10.5281/zenodo.7811426 (Leung, 2023a).

**Data availability.** MERRA-2 meteorological fields, land-surface fields, and dust emissions are available at https://disc.gsfc.nasa.gov/datasets?project=MERRA-2 (Earthdata, 2023; last access: 12 August 2022). SoilGrids soil texture and properties can be obtained from https://soilgrids.org/ (Soilgrids, 2023; last access: 20 December 2021). ESA CCI land cover can be obtained from https://www.esa-landcover-cci.org/?q=node/164 (ESA, 2023; last access: 21 June 2022). DustCOMM inverse analysis results are available at https://dustcomm.atmos.ucla.edu (UCLA, 2023; last access: 21 June 2022). The satellite-derived aeolian roughness data are available upon contacting Catherine Prigent.

**Supplement.** The supplement related to this article is available online at: https://doi.org/10.5194/acp-23-6487-2023-supplement.

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

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

**Acknowledgements.** 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 (https://doi.org/10.5065/D6RX99HX), 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.

## References

Adebiyi, A. A. and Kok, J. F.: Climate models miss most of the coarse dust in the atmosphere, Sci. Adv., 6, eaaz9507, https://doi.org/10.1126/sciadv.aaz9507, 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, https://doi.org/10.1002/2013MS000279, 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, https://doi.org/10.1029/2000JD900339, 2001.
Anderson, R. S.: Saltation of sand: a qualitative review with biological analogy, P. Roy. Soc. Edinb. B, 96, 149–165, https://doi.org/10.1017/S0269727000010903, 1989.

Andreotti, B., Claudin, P., and Pouliquen, O.: Measurements of the aeolian sand transport saturation length, Geomorphology, 123, 343–348, https://doi.org/10.1016/j geomorph.2010.08.002, 2010.

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

Bagnold, R. A.: The Transport of Sand by Wind, The Geogr. J., 89, 409–438, https://doi.org/10.2307/1786411, 1937.

Bagnold, R. A.: The Physics of Blown Sand and Desert Dunes, 1st Edn., Springer Netherlands, New York, 265 pp., https://doi.org/10.1007/978-94-009-5682-7, 1941.

Barchyn, T. E., Martin, R. L., Kok, J. F., and Hugen- holtz, C. H.: Fundamental mismatches between measurements and models in aeolian sediment transport prediction: The role of small-scale variability, Aeolian Res., 15, 245–251, https://doi.org/10.1016/j.aeolia.2014.07.002, 2014.

Barth, M. C., Rasch, P. J., Kiehl, J. T., Benkovitz, C. M., and Schwartz, S. E.: Sulfur chemistry in the National Cen- ter for Atmospheric Research Community Climate Model: Description, evaluation, features, and sensitivity to aque- ous chemistry, J. Geophys. Res.-Atmos., 105, 1387–1415, https://doi.org/10.1029/1999JD900773, 2000.

Bradley, E. F. and Mulhearn, P. J.: Development of veloc- ity and shear stress distribution in the wake of a porous shelter fence, J. Wind Eng. Ind. Aerod., 15, 145–156, https://doi.org/10.1016/0167-6105(83)90185-X, 1983.

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

Cakmur, R. V., Miller, R. L., and Torres, O.: Incorporating the ef- fect of small-scale circulations upon dust emission in an atmo- spheric general circulation model, J. Geophys. Res.-Atmos., 109, D7, https://doi.org/10.1029/2003JD004067, 2004.

Callot, Y., Marticorena, B., and Bergametti, G.: Geomor- phologic approach for modelling the surface features of arid environments in a model of dust emissions: appli- cation to the Sahara desert, Geodin. Acta, 13, 245–270, https://doi.org/10.1080/09853111.2000.11105373, 2000.

Chappell, A. and Webb, N. P.: Using albedo to reform wind erosion modelling, mapping and monitoring, Aeolian Res., 23, 63–78, https://doi.org/10.1016/j.aeolia.2016.09.006, 2016.

Chappell, A., Webb, N., Hennen, M., Zender, C., Ciais, P., Schep- anski, 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], https://doi.org/10.5194/gmd-2021-337, 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, https://doi.org/10.1111/j.1365-3091.1996.tb01509.x, 1996.

Chu, C. R., Parlange, M. B., Katul, G. G., and Albertson, J. D.: Probability density functions of turbulent velocity and temper- ature in the atmospheric surface layer, Water Resour. Res., 32, 1681–1688, https://doi.org/10.1029/96WR00287, 1996.

Claudin, P. and Andreotti, B.: A scaling law for aeo- lian dunes on Mars, Venus, Earth, and for subaque- ous ripples, Earth Planet. Sci. Lett., 252, 30–44, https://doi.org/10.1016/j.epsl.2006.09.004, 2006.

Comola, F., Gaume, J., Kok, J. F., and Lehning, M.: Cohesion- Induced Enhancement of Aeolian Saltation, Geophys. Res. Lett., 46, 5566–5574, https://doi.org/10.1029/2019GL082195, 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, https://doi.org/10.1029/2019GL085739, 2019b.

Cowie, S. M., Marsham, J. H., and Knippertz, P.: The im- portance of rare, high-wind events for dust uplift in northern Africa, Geophys. Res. Lett., 42, 8208–8215, https://doi.org/10.1002/2015GL065819, 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, https://doi.org/10.5194/soil-5-137-2019, 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 Forecast- ing (WRF) model: Assessment of dust emission parameter- izations and input parameters for source regions in Cen- tral and East Asia, J. Geophys. Res.-Atmos., 114, D14, https://doi.org/10.1029/2008JD011236, 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, https://doi.org/10.1029/2019GL086186, 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., Blan- ton, C., Dunne, K. A., Dupuis, C., Durachta, J., Dussin, R., Gau- thier, P. P. G., Griffies, S. M., Guo, H., Hallberg, R. W., Har- rison, M., He, J., Hurlin, W., McHugh, C., Menzel, R., Milly, P. C. D., Nikonov, S., Paynter, D. J., Ploshay, J., Radhakrish- nan, A., Rand, K., Reichl, B. G., Robinson, T., Schwarzkopf, D. M., Sentman, L. T., Underwood, S., Vahlenkamp, H., Win- ton, 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 Char- acteristics, J. Adv. Model. Earth Syst., 12, e2019MS002015, https://doi.org/10.1029/2019MS002015, 2020.

Dupont, S., Bergametti, G., Marticorena, B., and Simoëns, S.: Modeling saltation intermittency, J. Geophys. Res.-Atmos., 118, 7109–7128, https://doi.org/10.1002/jgrd.50528, 2013.

Durán, O., Claudin, P., and Andreotti, B.: On aeolian transport: Grain-scale interactions, dynamical mech- anisms and scaling laws, Aeolian Res., 3, 243–270, https://doi.org/10.1016/j.aeolia.2011.07.006, 2011.
Earthdata: Data Collections, Ges Dic [data set], https://disc.gsfc.nasa.gov/datasets?project=MERRA-2 (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, https://doi.org/10.1038/nature04058, 2005.

ESA: Land Cover CCI Product User Guide Version 2 Technical Report, European Space Agency, http://maps.elie.ucl.ac.be/CCI/viewer/download/ESACCI-LC-Ph2-PUGv2_2.0.pdf (last access: 21 June 2022), 2017.

ESA: Download CCI LC Products, ESA [data set], https://www.esa-landcover-cci.org/?q=node/164 (last access: 21 June 2022), 2023.

Evan, A. T.: Surface Winds and Dust Biases in Climate Models, Geophys. Res. Lett., 45, 1079–1085, https://doi.org/10.1002/2017GL076353, 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, https://doi.org/10.1002/2014GL060545, 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, https://doi.org/10.1016/j.aeolia.2015.01.001, 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, https://doi.org/10.1002/2016GL071016, 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, https://doi.org/10.1016/j.atmosenv.2006.09.048, 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, https://doi.org/10.1007/s00585-999-0149-7, 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, https://doi.org/10.1029/2021MS002909, 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, https://doi.org/10.1002/2016MS000823, 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, https://doi.org/10.1038/s41561-022-00901-w, 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, https://doi.org/10.1175/JCLI-D-16-0758.1, 2017.

Gillette, D. A. and Passi, R.: Modeling dust emission caused by wind erosion, J. Geophys. Res.-Atmos., 93, 14233–14242, https://doi.org/10.1029/JD093iD11p14233, 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, https://doi.org/10.1029/2000JD000053, 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, https://doi.org/10.1029/2012RG000388, 2012.

Greeley, R. and Iversen, J. D.: Wind as a Geological Process: On Earth, Mars, Venus and Titan, CUP Archive, 356 pp., https://doi.org/10.1017/CBO9780511573071, 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, https://doi.org/10.1029/97JE00518, 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, https://doi.org/10.1029/2004JD005037, 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, https://doi.org/10.1016/j.rse.2015.01.021, 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, https://doi.org/10.1080/014311698215333, 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, https://doi.org/10.1029/2019GB006448, 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, https://doi.org/10.1371/journal.pone.0105992, 2014.

Hengl, T., Heuvelink, G. B. M., Kempen, B., Leenaars, J. G. B., Walsh, M. G., Shepherd, K. D., Sila, A., MacMillan, R. A., Je-
sus, 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, https://doi.org/10.1371/journal.pone.0125814, 2015.

Hengl, T., Jesus, J. M. de, Heuvelink, G. B. M., González, M. R., Kilibarda, M., Blagotíæ, A., Shangguan, W., Wright, M. N., Geng, X., Bauer-Marschallinger, B., Guvara, 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, https://doi.org/10.1371/journal.pone.0169748, 2017.

Hillel, D.: Fundamentals of Soil Physics, Elsevier, https://doi.org/10.1016/C2009-0-03109-2, 1980.

Huang, X. and Foroutan, H.: Effects of Non-Photosynthetic Vegetation on Dust Emissions, J. Geophys. Res.-Atmos., 127, e2021JD035243, https://doi.org/10.1029/2021JD035243, 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, https://doi.org/10.5194/acp-11-7781-2011, 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, https://doi.org/10.1002/2016JD025939, 2017.

Iversen, J. D. and White, B. R.: Saltation threshold on Earth, Mars and Venus, Sedimentology, 29, 111-119, https://doi.org/10.1111/j.1365-3091.1982.tb01713.x, 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, https://doi.org/10.1126/science.1105959, 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, https://doi.org/10.1016/j.earscirev.2021.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, https://doi.org/10.2151/sola.2021-042, 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, https://doi.org/10.1029/2004JF000281, 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, https://doi.org/10.5194/acp-6-1815-2006, 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, https://doi.org/10.5194/essd-9-927-2017, 2017.

Klose, M. and Shao, Y.: Stochastic parameterization of dust emission and application to convective atmospheric conditions, Atmos. Chem. Phys., 12, 7309-7320, https://doi.org/10.5194/acp-12-7309-2012, 2012.

Klose, M. and Shao, Y.: Large-eddy simulation of turbulent dust emission, Aeolian Res., 8, 49-58, https://doi.org/10.1016/j.aeolia.2012.10.010, 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, https://doi.org/10.1002/2014JD021688, 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, https://doi.org/10.5194/gmd-14-6403-2021, 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, https://doi.org/10.5539/jgg.v9n3p1, 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, https://doi.org/10.1029/2010GL043646, 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, https://doi.org/10.1088/0034-4885/75/10/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, https://doi.org/10.5194/acp-14-13023-2014, 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, https://doi.org/10.5194/acp-14-13043-2014, 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, https://doi.org/10.1038/ngeo2912, 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, https://doi.org/10.1038/s41467-017-02620-y, 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, https://doi.org/10.5194/acp-21-8127-2021, 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, https://doi.org/10.5194/acp-21-8169-2021, 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, https://doi.org/10.1038/s43017-022-00379-5, 2023.

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

Laurent, B., Marticorena, B., Bergametti, G., Chazette, P., Maiganan, 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 POLD-ER/ADEOS 1 surface products, J. Geophys. Res.-Atmos., 110, D18, https://doi.org/10.1029/2004JD005013, 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, https://doi.org/10.5194/gmd-9-2973-2016, 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, https://doi.org/10.5194/gmd-12-131-2019, 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, https://doi.org/10.5194/gmd-16-1009-2023, 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 1 climate models - Part I: description and evaluation against inverse modeling emissions” (1.0), Zenodo [code], https://doi.org/10.5281/zenodo.7811426, 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], https://doi.org/10.5194/egusphere-2023-823, 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, https://doi.org/10.1002/jgrf.20040, 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, https://doi.org/10.5194/acp-21-3973-2021, 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, https://doi.org/10.5194/gmd-15-8181-2022, 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, https://doi.org/10.5194/essd-10-219-2018, 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, https://doi.org/10.2136/sssaj2014.06.0235, 2014.

Lin, S.-J.: A “Vertically Lagrangian” Finite-Volume Dynamical Core for Global Models, Mon. Weather Rev., 132, 2293–2307, https://doi.org/10.1175/1520-0493(2004)132<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, https://doi.org/10.5194/acp-6-3315-2006, 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, https://doi.org/10.1029/2000JD900839, 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, https://doi.org/10.1016/j geomorph.2004.03.009, 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, https://doi.org/10.1029/2004GB002402, 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 https://doi.org/10.1029/2005JD006653, 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, https://doi.org/10.5194/acp-10-10875-2010, 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, https://doi.org/10.5194/gmd-10-2397-2017, 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, https://doi.org/10.1029/95JD00690, 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, https://doi.org/10.1029/96JD02964, 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, https://doi.org/10.1080/0143116031000116976, 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, https://doi.org/10.1029/2006JF000462, 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, https://doi.org/10.1029/2017JF004416, 2018.

Martin, R. L. and Kok, J. F.: Size-Independent Susceptibility to Transport in Aeolian Saltation, J. Geophys. Res.-Earth Surf., 124, 1658–1674, https://doi.org/10.1029/2019JF005104, 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, https://doi.org/10.1016/j.rse.2006.01.016, 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, https://doi.org/10.1029/2019JD031807, 2020.

McKenna Neuman, C. and Sanderson, S.: Humidity control of particle emissions in aeolian systems, J. Geophys. Res.-Earth Surf., 113, F2, https://doi.org/10.1029/2007JF000780, 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 Com-

munity Land Model 5.1, Geosci. Model Dev., 15, 2365–2393, https://doi.org/10.5194/gmd-15-2365-2022, 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, https://doi.org/10.5194/acp-22-11889-2022, 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, https://doi.org/10.5194/gmd-14-4249-2021, 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, https://doi.org/10.1175/JTECH-D-17-0173.1, 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, https://doi.org/10.1002/jgrd.50313, 2013.

Miller, R. L. and Tegen, I.: Climate Response to Soil Dust Aerosols, J. Climate, 11, 3247–3267, https://doi.org/10.1175/1520-0442(1998)011<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, https://doi.org/10.1029/2003JD004085, 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, https://doi.org/10.1023/A:1022830119554, 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, https://doi.org/10.5194/gmd-5-581-2012, 2012.

Nikuradse, J.: Laws of flow in rough pipes, Tech. Rep. NACATM-1292, National Advisory Committee for Aeronautics, https://ntrs.nasa.gov/citations/19930093938 (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, https://doi.org/10.1029/2007JF000758, 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/TN-503+STR), https://doi.org/10.5065/D6RR1W7M, 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, https://doi.org/10.1029/2019RG000679, 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, https://doi.org/10.1007/BF02186086, 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, https://doi.org/10.1029/2019JD030248, 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, https://doi.org/10.5194/acp-22-8659-2022, 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, https://doi.org/10.5194/gmd-3-753-2010, 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, https://doi.org/10.1002/2014JF003177, 2014a.

Pierre, C., Bergametti, G., Marticorena, B., Abdourhamane-Touré, 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, https://doi.org/10.1016/j.catena.2014.06.006, 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, https://doi.org/10.1029/2004JD005370, 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, https://doi.org/10.5194/amt-5-2703-2012, 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, https://doi.org/10.1073/pnas.96.7.3396, 1999.

Pu, B. and Ginoux, P.: How reliable are CMIP5 models in simulating dust optical depth?, Atmos. Chem. Phys., 18, 12491–12510, https://doi.org/10.5194/acp-18-12491-2018, 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, https://doi.org/10.5194/acp-20-55-2020, 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, https://doi.org/10.1029/2018JD028925, 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, https://doi.org/10.1103/PhysRevE.105.054902, 2022.

Raupach, M. R.: Drag and drag partition on rough surfaces, Bound.-Lay. Meteorol., 60, 375–395, https://doi.org/10.1007/BF00155203, 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, https://doi.org/10.1029/92JD01922, 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, https://doi.org/10.1175/JCLI-D-16-0720.1, 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, https://doi.org/10.1029/2011JD016794, 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, https://doi.org/10.1002/grl.50409, 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, https://doi.org/10.5194/acp-14-5735-2014, 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, https://doi.org/10.1073/pnas.101122798, 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, https://doi.org/10.1029/97JD00341, 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, https://doi.org/10.1002/2013MS000293, 2014.

Shao, Y.: A model for mineral dust emission, J. Geophys. Res.-Atmos., 106, 20239–20254, https://doi.org/10.1029/2001JD900171, 2001.

Shao, Y.: Simplification of a dust emission scheme and comparison with data, J. Geophys. Res.-Atmos., 109, D4, https://doi.org/10.1029/2003JD004372, 2004.

Shao, Y.: Physics and Modelling of Wind Erosion, Springer Science & Business Media, 459 pp., https://doi.org/10.1007/978-1-4020-8895-7, 2008.

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