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

. 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 ﬂuctuations. 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 ﬂuxes 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 ﬁelds. 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 ﬁnd that our revised dust emission parameterization can substantially improve dust emission simulations in GCMs and LSMs.


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 dustradiation 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., , 2010Hamilton 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 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;. 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 . 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 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 landsurface roughness (e.g., Greeley et al., 1997;Roujean et al., 1997;Marticorena et al., 2004;Prigent et al., 2005Prigent et al., , 2012, especially over arid desert regions over which rocks, pebbles, and gravel dominate the roughness. However, with a few recent notable exceptions that attempted to represent the roughness effects of both rocks and vegetation (e.g., Darmenova et al., 2009;Foroutan et al., 2017;Klose et al., 2021), studies often omitted the drag partition effect either due to vegetation (e.g., Menut et al., 2013) or due to rocks (e.g., Wu et al., 2016;LeGrand et al., 2019;Tai et al., 2021). To resolve these issues, we propose a new approach that combines the drag partition effects of both elements, leveraging satellite scatterometer measurements to quantify the surface rock roughness and using observable vegetation and land-surface variables to quantify the surface vegetation roughness.
The third key piece of fundamental dust emission physics not accounted for by many models is the effect of turbulencedriven high-frequency wind variability on dust emissions. Most current GCMs assume a constant wind speed (and thus a constant emission flux) within the relatively large model time step, e.g., 30 min (e.g., Rahimi et al., 2019;Dunne et al., 2020). However, in reality, high-frequency turbulent fluctuations cause the wind speed to fluctuate within a time step (from seconds to minutes). Because dust emissions scale nonlinearly with wind speed, this causes highly uneven and fluctuating dust emission fluxes (Durán et al., 2011). Even more importantly, turbulent wind fluctuations can sweep across the dust emission threshold multiple times and shut off dust emissions intermittently within one model time step, resulting in strong dust emission intermittency (Comola et al., 2019b). Even regional climate models (RCMs), which typically use a smaller time step (e.g., < 1 min), do not resolve turbulence unless they are run in the computationally expensive large-eddy-simulation (LES) mode (e.g., WRF-LES). Omitting turbulence by GCMs and RCMs thus causes either an overestimate or an underestimate of dust emissions, especially over marginal source regions where winds fluctuate around the high emission threshold, as models do not account for the cessations or initiations of dust emissions due to turbulent fluctuations. To account for the instantaneous wind fluctuations, a dynamical approach is to derive a probability density function (PDF) for the instantaneous momentum flux using LES, which is then used for quantifying instantaneous dust emission fluctuations (Klose and Shao, 2012;Klose et al., 2014). A parameterization approach is to use the Monin-Obukhov similarity theory (MOST) to relate the standard deviation of the instantaneous wind to the boundary-layer dynamical variables (Comola et al., 2019b). In this study, we will account for turbulent dust emissions by following Comola et al. (2019b), which showed significant improvements in representing the small-magnitude saltation and dust fluxes that are particularly important over marginal source regions.
In addition to these issues of models missing some of the fundamental physics of dust emission, a central issue in modeling the global dust cycle is that dust emissions are gridresolution-dependent because of the nonlinear dependence of dust emissions to meteorological fields and land-surface variables. Since dust emission varies nonlinearly with wind speed and has an even more complex relation to the soil moisture (Gillette and Passi, 1988;Fécan et al., 1999;Shao, 2001;Kok et al., 2014a), the total regional and global emissions can vary significantly with grid resolution (Ridley et al., 2013;Meng et al., 2021). For instance, modeled emissions were found to increase by ∼ 29 % from a 1 • × 1 • to 0.25 • × 0.25 • resolution (Ridley et al., 2013;Feng et al., 2022). As a consequence, GCMs and ESMs often need to tune emissions separately for different grid resolutions to match observational dust budgets (Ginoux et al., 2001;Zender et al., 2003a;Albani et al., 2014;Kok et al., 2014b;Chappell et al., 2021). This issue occurs because current GCM grid sizes of ∼ 1 • or 100 km cannot resolve the spatial scales of ∼ 1 m to ∼ 1 km over which soil properties and wind speeds change (Ridley et al., 2013). When adopting coarse grid resolutions, coarser modeled meteorological fields (for GCMs) or spatially averaged input meteorological fields (for chemical transport models or CTMs) will smooth out the local wind extrema, possibly causing wind speeds to fall below the dust emission threshold. As a result, the coarse modeled winds usually result in strong GCM emissions underestimations (Ridley et al., 2013). The same smoothing problem also occurs for soil moisture for instance, with its maxima smoothed out leading to an overestimation of dust emissions. Thus, although some GCMs and ESMs recently implemented more physical schemes (Zhao et al., 2022), their inability to resolve the small scales still causes challenges for capturing the accurate spatial distributions of dust emissions (Meng et al., 2021). For the same reason, GCMs tend to neglect smallscale 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 scaledependence 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 turbulencedriven 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 6490 D. M. Leung et al.: A new desert dust emission scheme -Part 1: Description 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.

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.

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 ρ a , and soil moisture w (Greeley et al., 1997;Shao and Lu, 2000). Although there are multiple data sources of globally gridded products for ρ a and w, there are relatively few efforts on obtaining globally gridded 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 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): where u * ft0 is the "dry" fluid threshold friction velocity (in m s −1 ) on a smooth and bare surface as a function of air density ρ a (long, lat, t) (longitude, latitude, time) and D p (long, lat), which in this study will be the median diameter 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 ≥ 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: where g = 9.81 m s −2 is the gravitational acceleration, ρ p = 2650 kg m −3 is the typical soil particle density, ρ a is in kg m −3 , and A = 0.0123 and γ = 1.65 × 10 −4 kg s −2 are empirical constants accounting for the aerodynamic forces and interparticle forces, respectively. Assuming an air density ρ a = 1.225 kg m −3 , Eq.
(2) will yield the smallest u * ft0 of 0.21 m s −1 at D p = 80 µm. For larger sizes, the particles are heavier to lift; for smaller sizes, the particles are more strongly bound by interparticle forces. An alternative parameterization for u * ft0 is the Iversen and White (1982) scheme (hereafter I&W82). They derived a similar solution as S&L00 but further considered the effects of soil particle size to the airflows, characterized by the particle Reynolds number Re p . I&W82 calculates u * ft0 = u * ft0 ρ p , ρ 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 where ν 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 ρ 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 µm for the threshold calculation following Zender et al. (2003a) (Oleson et al., 2013; also see the latest version of CLM5.0 technical note on 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 ρ a in CESM (e.g., u * ft0 = 0.21 m s −1 for ρ a = 1.225 kg m −3 at D p = 75 µm).
Most models parameterize the effect of soil moisture f m on u * ft following Fécan et al. (1999): where w(long, lat, t) is the gravimetric soil moisture (kg kg −1 ) in the shallowest soil layer (see Sect. S1 and Fig. S1 for the relation between volumetric and gravimetric moisture); w t (long, lat) is the threshold gravimetric water content above which u * ft increases; f clay (long, lat) is the fraction of clay content in the topmost layer of soil between zero and one; %clay = 100f 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): 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 ≥ 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).

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 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 Marticorena and Bergametti, 1995;Marticorena et al., 1997Marticorena et al., , 2006Foroutan et al., 2017;Webb et al., 2020) u * ft = u * ft0 (D p , ρ a )f m (w)/F eff , where F 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 6492 D. M. Leung et al.: A new desert dust emission scheme -Part 1: Description 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 * : 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 mountains/orography, rocks, and plants, as well as down to soil particles, Menut et al. (2013) distinguished three distinct roughness lengths describing different sizes of roughness. First, the aerodynamic momentum roughness length z 0m mainly represents roughness due to large-scale orography, forests, and/or urbanization (with sizes of 10-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.
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 δ, they derived F eff as a function of z 0a and z 0s : 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  (hereafter R93), which primarily used the roughness density λ to quantify roughness. λ is defined as the total frontal area of roughness elements divided by the area of land A : λ = nhb/A, where h and b are the obstacle height and width, and n is the number of obstacles within the area. Knowing the geometric and aerodynamic properties of the roughness elements, R93 showed that the drag force of the exposed area τ S is related to the total drag force τ , given λ, the roughnesselement basal area-to-frontal area ratio σ , and the ratio of the roughness element-to-surface drag coefficient β: where m is a geometric parameter to account for the spatial variability of τ S on the erodible surface. Raupach then applied this ratio to the dust emission threshold (per Eq. 6a). Many later studies used the R93 parameterization for plants (specifically shrubs) with prescribed σ , m, and β (Darmenova et al., 2009;Xi and Sokolik, 2015). λ, however, is related to the abundance of obstacles and is thus spatially variable, and thus far there are no globally gridded datasets of λ available. Most studies thus related grid-scale λ to other grid-scale properties; for instance, Shao et al. (1996) linked λ to the vegetation cover fraction f v using in situ observations: where c λ is a proportionality constant. Gridded λ could thus be obtained from gridded satellite retrievals of vegetation cover (Gutman and Ignatov, 1998;Wu et al., 2016;Foroutan et al., 2017) or parameterized as a function of other gridded land-surface variables such as the leaf area index (LAI) (e.g., Klose et al., 2021). Later studies have attempted to improve Raupach's parameterization, and newer schemes relating F eff and λ have emerged (e.g., Okin, 2008). Many previous modeling studies have not accounted for the drag partition effects of both rocks and vegetation on dust emissions (e.g., Ginoux et al., 2001;Tegen et al., 2002;Zender et al., 2003a;Kok et al., 2014b). Many past studies either accounted for only the drag partitioning by rocks (e.g., Marticorena et al., 2006;Menut et al., 2013) or by vegetation (e.g., Shao et al., 2011b;Wu et al., 2016), mainly because it is very challenging to use proxies of both rocks and vegetation in either the M&B95 or R93 scheme. For instance, R93 was historically mostly used for modeling vegetation but not rock drag partitioning because there was no dataset of the λ due to rocks. Similarly, vegetation roughness is historically mostly represented by λ rather than 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 λ 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 λ 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 λ 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 satellitederived 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 λ from both schools of drag partition parameterizations, quantifying the drag partitions of rocks and plants into one hybrid drag partition factor F eff .

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: 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; ϕ = 10 13.4f 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 clay ; u * t is the dust emission threshold (in m s −1 ; Z03 used u * ft as u * t ); and f 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 bare , Mahowald et al. (2006) used a simple parameterization in which f bare is a pure function of LAI (neglecting the effects of other objects such as snow, rocks, and buildings): 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 ∝ f bare ) drops more rapidly with LAI such that the drag partition effect is also incorporated in the f bare term. However, since we are considering F eff in this study, we can set a more realistic LAI thr value such that f 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 √ ρ a /ρ a0 , which is u * ft scaled to the standard air density of ρ a0 = 1.22 kg m −3 : where C d (long, lat, t) is the time-varying dust emission coefficient or soil erodibility coefficient, C d0 = (4.4±0.5)×10 −5 , C e = 2.0±0.3, and u * st0 = 0.16 m s −1 . Furthermore, K14 derived a new dust emission equation for F d (kg m −2 s −1 ): where C κ = 2.7 ± 1.0, C tune = 0.05 is the proportionality constant, f bare is modeled by Eq. (11), u * s = u * in K14, and u * t is again the emission threshold (K14 assumed for simplicity that u * t = u * it = u * ft ). κ is the fragmentation exponent which quantifies the sensitivity of F d to u * s . Here we limit the value of κ to 3 in order to prevent excessive sensitivity of the model to wind speeds, which can be problematic around topography. From Eq. (12), 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).

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

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.

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 D p of the soil PSD. Thus, we ideally need a global gridded map of 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 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 D p , we use a compilation of literature measurements to explore and construct relationships between 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, 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, wetsieved soil texture f a , and other soil properties to investigate their statistical relations and infer a new global distribution of 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 D p by calculating the 50th percentile of the dry-soil PSD.

Deriving a global soil median diameter map
We classify the datasets into arid and nonarid groups, since we are primarily interested in D p over desert regions (although we also display the soil behaviors over nonarid regions). We follow past studies (Mahowald et al., 2006(Mahowald et al., , 2010Kok et al., 2014b) which defined arid (or dust emission) regions using the criterion of LAI smaller than a threshold LAI thr , 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 D p and the soil properties (see Figs. 1b and S2). Figure  1b shows a scatterplot of D p versus the sum of the soil component that produces substantial cohesion, namely the silt and clay fractions (f silt+clay = f silt + f clay ). The data exhibit distinctly different trends for nonarid versus arid soils: for nonarid soils, D p increases from 100 µm to greater than 1000 µm with f silt+clay (regression p value = 7.3 × 10 −6 ), likely due to increasing cohesion with increasing clay and silt content. In contrast, D p for arid soils shows a small and statistically insignificant increasing trend with f silt+clay (p value = 0.77) with a smaller D p variability (50-250 µm). This flat trend indicates that f silt+clay does not effectively explain the median diameter of aggregated soil particles in arid regions. We examined the relationships of D p with the individual fractions of sand, silt, and clay, as well as with other soil properties including SOC, pH, and 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 D p over the arid regions across different soil textures. We thus use a constant D p0 as an approximation for arid regions. From Fig. 1b, we summarize the relationship between D p and f silt+clay as where 0 = 7.8 ± 3.7 µm, 1 = 124 ± 36 µm, D p0 = 127 ± 47 µm, and LAI thr = 1 as specified in Eq. (11). This empirical formula suggests that some models' assumptions of the relationship between D p and soil texture were inaccurate (e.g., Table 2 of Laurent et al. (2008) assumed D p decreases with f 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 * ft0 from 0.204 to 0.216 m s −1 (given ρ a = 1.225 kg m −3 ), which slightly decreases global dust emissions by 18 % (see Sect. 4.1). The uncertainty in D p0 = 127±47 µm translates to an uncertainty of u * ft0 between 0.204 to 0.234 m s −1 . We then project our derived relation between D p and f silt+clay on the available soil texture and properties database. We employ global soil properties data from the SoilGrids database (Hengl et al., 2014(Hengl et al., , 2015(Hengl et al., , 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).  regions, D p increases with f silt+clay . Our derived D p values are largely consistent with the site 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 D p -f silt+clay relation over nonarid regions (blue points). Figure 1d shows that Eq. (14) gives satisfactory agreement in predicting global D p , but dust emission modeling will depend exclusively on the predicted 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 D p and various soil properties over arid regions.

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

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 • z 0a 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 z 0a yielded an empirical formula for extrapolating a global dataset of backscattering signal to global z 0a . We use here the global aeolian z 0a dataset from Prigent et al. (2005) (hereafter Pr05), which contains the climatological monthly mean z 0a (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 z 0a 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 z 0a 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 z 0a (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 z 0a and z 0s . They assumed that this equation is valid for roughness elements that are not too closely spaced (small wake), i.e., z 0a < 1 cm (Darmenova et al., 2009). Here we use their semiempirical equation to quantify the drag partition due to rocks, f eff,r (also see Eq. 8): 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) 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 z 0s from our de-rived global dataset of D p in Sect. 3.1 using Eq. (7). When nonerodible roughness elements are abundant over a surface, z 0a z 0s and f eff,r 1, causing the sheltering of the bare soil from the wind; when there are few roughness elements, z 0a is small and close to z 0s , and thus f eff,r approaches 1. In Fig. 2a, the red areas with very small z 0a are the most susceptible regions for dust emission. Figure 2a also shows that most arid and semiarid regions have z 0a < 0.2 cm, such that Eq. (15) can follow the criterion (z 0a < 1 cm) in Darmenova et al. (2009) well. Figure 2b shows the global f eff,r over arid regions, which is dominated by the spatial pattern of z 0a in Fig. 2a given f eff,r is governed purely by z 0a .

Drag partition due to roughness of vegetation
Unlike the very static and slowly evolving rock roughness, vegetation changes temporally. To include the effect of these dynamic vegetation changes on the drag partition, we follow the approach of Okin (2008) (hereafter O08), which uses unvegetated gap size (the distance between neighboring plants) to characterize the variability of the reduced wind stress. O08 argued that his scheme represents an advancement over the classical R93 scheme, since R93 uses the roughness density (or lateral cover) λ, which only quantifies how much roughness is on a surface but not how that roughness is spatially distributed. O08 pointed out that, given the same λ, roughness elements divided into small blocks spread over the soil surface would be more effective than elements stacked up like a telephone pole in partitioning wind stress (see Fig. 3 in Okin, 2008). Okin argued that since Raupach's model neglects the spatial variability of λ, the resulting simulated emission flux using the R93 scheme in Okin's paper decreased rapidly with increasing λ and unrealistically reached zero at relatively low λ. To partially compensate for this error, R93 introduced a tuning parameter m (Eq. 9a), serving to reduce the effective λ and thereby reducing the rapid decrease in dust flux. However, m is a tuning parameter not derived from first principles, and it is not clear how m changes over different surface conditions. Therefore, we use the O08 model here to better characterize the spatial variability of wind stress and the resulting dust emissions.
Here we describe the O08 scheme and adapt it for use in LSMs and GCMs. O08 assumes u * drops significantly when encountering a roughness element (plant) and gradually recovers at the lee (downwind region) of the plant as a function of distance x, following where x/ h is the dimensionless downwind distance from an obstacle normalized by vegetation height h (m), f 0 = 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  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.
of distance x is Note that in the limit of x/ h → ∞, u * s → 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: 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 ≡ 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 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 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 where we assume LAI thr = 1 (in Eq. 11).
The assumption of f v ∼ 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 −2 s −1 ) by locally integrating the emission F d following the spatially varying u * s over the normalized distance x/ h: where F d (kg m −2 s −1 ) is the local emission as a function of u * s which is itself a function of x/ h. In the integration, 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 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 f eff,v 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: 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: This elegant formula conveys a clear physical intuition: if the obstacle does not effectively dissipate momentum (f 0 → 1), f eff,v → 1; if land is densely covered by vegetation (gap length K → 0), f eff,v → 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 . 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.

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 ( 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 : .
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 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é Depression, 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.

Parameterizing the dust emission intermittency
The above improvements enable a more accurate calculation of emission when wind speeds are sufficient to initiate dust emission. Next, we will improve the calculation of the resulting dust emission flux by accounting for the effects of boundary-layer turbulence on dust emission intermittency. Dust emission intermittency exists because saltation is driven by turbulent surface winds, which exhibit strong spatiotemporal fluctuations in speed and direction. Instantaneous winds can thus pass within short timescales across the emission thresholds for initiating or ceasing saltation (Martin and Kok, 2018). Consequently, saltation can be highly intermittent (Comola et al., 2019b), with pronounced variability on timescales of seconds to hours (Dupont et al., 2013). In contrast, existing dust emission parameterizations describe saltation as uniform in time and space and driven by a constant downward momentum flux within a model time step. The disconnect between the reality of intermittent dust emissions and uniform emissions in current theories is likely contributing to the poor performance of dust emission simulations (Barchyn et al., 2014;Todd et al., 2008). Comola et al. (2019b) argued that the intermittency effect is more prevalent for regions with low-intensity dust emissions when u * s is regularly fluctuating around the threshold to turn on or shut off dust emissions. Neglecting intermittent dust emissions in current models thus likely degrades the accuracy of dust emission simulations for arid regions during low-wind periods and for marginal dust source regions such as semiarid areas (since soil cohesion increases u * ft but does not affect u * it ), which dominate in much of the Southern Hemisphere (Ginoux et al., 2012;Ito and Kok, 2017).
Accounting for the intermittency effect on dust emission fluxes is complicated by the hysteresis of dust emission due to the existence of double thresholds for dust emission physics. The instantaneous wind at the saltation level u s (we use the tilde to denote instantaneous quantities and take away the asterisk to denote winds at the saltation level of z sal ∼ 0.1 m using Eq. (S4a) instead of a velocity scale) needs to exceed the fluid threshold u ft (also defined at the saltation level) to initiate saltation, but it only needs to exceed a smaller impact threshold u it to sustain it (Kok et al., 2012;Martin and Kok, 2018;Comola et al., 2019b). Wheñ u s at a moment lies between both thresholds (u it <ũ s < u ft ), saltation is active if transport was more recently initiated (ũ s > u ft ) and inactive if transport was more recently terminated (ũ s < u 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ũ s within a model time step) is between u it and u 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 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 it and the Monin-Obukhov length L. C19 showed that when a dust emission equation employs u it and accounts for the intermittency effect, it can successfully capture the magnitudes of small dust fluxes otherwise missed by models using u 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 where C d is still a function of u * st , and u * st = u * ft √ ρ a /ρ a0 is the same standardized fluid threshold as in the default K14, and u * it is computed using Eq. (5). Because u * it <u * ft , this modified equation allows more small dust fluxes over the marginal source regions that are otherwise missed by employing u * ft as the threshold (see Figs. 7g-h).
Next, we account for the intermittency effect by introducing the intermittency factor η, which is the fraction of time that saltation is active in a model time step (e.g., ∼ 30 min). η corrects the horizontal sand saltation flux, which scales with dust emission flux (Shao et al., 1993), thereby also representing the fraction of time that dust emission is active in a model time step. C19 accounts for the effect of intermittency by multiplying dust emission by η as follows: where η ∈ [0, 1]. Note that C19 parameterizes η using wind information at the typical saltation height of z sal = 0.1 m instead of the velocity scales. η 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 σũ s of the instantaneousũ s (Eq. 23): σũ s is defined given thatũ s can be described by a normal distribution (Chu et al., 1996), with its mean being the model time step mean (at 0.1 m) u s and its standard deviation σũ s . From Eq. (22c), η → 1 when u s u ft (active emission for the whole time step) and η → 0 when u s 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 u s sweeping across the thresholds and the more dichotomous η behaves (either zero or 1). If u s is very close to u ft or u it , or is indeed between them (u it <ũ s < u ft ), the frequency of crossing the threshold is determined by the magnitude of the turbulent fluctuation σũ s . σũ s is parameterized using the similarity theory (Panofsky et al., 1977): 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 * , ρ a , sensible heat flux H , and temperature T for our simulations (see Sect. S3). In boundary-layer dynamics, turbulence is generated by mechanical shear and buoyancy (Stull, 1988). From Eq. (23), high-frequency wind fluctuations σũ s increase with shear (u * s > 0) and buoyancy (L < 0). A larger σũ s makes it easier forũ s to sweep across u it and shut off dust emission, leading to η < 1. In a time step, if u s u ft + σũ s , u s will be unlikely to sweep across u it , and η will approach 1. If u s is slightly larger than u ft , the instantaneousũ s will be likely to sweep across u it , leading to η < 1. In the hysteresis regime (u it < u s < u ft ), η will be around 0.3-0.7, sincẽ u s will sweep across both thresholds given σũ s , leading to a reduced emission flux (meanwhile, other parameterizations predict a zero emission flux since they use u ft only). When u s < u it , η could also be greater than zero when σũ s is large enough so that the instantaneousũ 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 η to yield the dust emission with intermittency effect F d,η 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 D p0 = 127 µm. Its spatial variability is purely a function of ρ a (Eq. 2). u * it is around ∼ 0.16-0.24 m s −1 , and higher u * it implies higher altitude. A lower u * it leads to a smaller aerodynamic drag force from the airflows given the same wind speed; conversely, there will be a large u * it for soils over low ρ a regions. Figure 5b shows u * ft (Eq. 1). It varies between 0.2 and 0.9 m s −1 . Its spatial variability is dictated by the spatial variability of soil moisture w (see Fig. S1). Regions with the lowest u * ft are the driest places in the world, which are all deserts. Regions with the highest u * ft are wet soils covered by rainforests, boreal forests, tundras/permafrosts, and snow. Figure 5c shows the ratio of 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 ∼ 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. Semiarid 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 shows hourly mean η as a function of the hourly mean u s . It demonstrates the properties of η discussed above: e.g., when u s > u ft , η → 1 and when u s < u it , η → 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 , η 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 η values of ∼ 0.5 only, and for the highest η 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 σũ s ≈ 12 1/3 u * s . Figure 6b shows the global spatial distribution of the annual mean η for the year 2006, averaged across time steps during which saltation and dust emissions are occurring over the grid (and thus η during time steps when F d = 0 is not counted). Most marginal sources have small η < 0.3 (red color), indicating dust emissions are fluctuating and intermittent. These emissions may not be existent in other LSMs employing u * ft in the dust emission equation (but also dependent on their threshold tuning). Over these regions, because of the intermittent shut-offs of emissions, the emissions need to be scaled down by the fraction of time η which is also missed by other LSMs. Intermittency is thus critical in accounting for emissions over semiarid regions. Regions with high η (more continuous emission; light blue color) include El Djouf; the Bodélé Depression; the Libyan-Nubian Desert over Libya, Egypt, and Sudan; the Rub' al Khali Desert within the greater Arabian Desert; the Lut Desert in Iran; the Taklamakan Desert; and the Strzelecki Desert in Australia. During saltation, these regions tend to have wind episodes further away from the thresholds, leading to a high η. However, regions with high η are not necessarily regions with the highest dust emissions (Fig. 7a-b) because their saltation frequency could be low, or the strength of their emission fluxes is limited by other factors such as the fragmentation exponent and soil moisture. In Fig. S5, we further show the factor η averaged over all time steps of 2006, including periods of no emissions over the grid. Since η is close to 0 when there is no emission, Fig. S5 shows much smaller η compared to Fig. 6b.

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  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 lowresolution dust emissions to match the spatial variability of high-resolution emissions.

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 D p0 ∼ 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;, 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 20 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 η. 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 D p0 = 127 µm (experiment II). As described in Sect. 3.1, a globally larger 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 D p0 is relatively mild (of the largest order of 0.001 kg m −2 yr −1 ) compared to the drag partition effect and intermittency ( Fig. 7e and g; of the largest order of ∼ 0.1 kg m −2 yr −1 ). Figure 7d shows the ratio map of the normalized emissions of experiment II to those of experiment I. Employing a new globally constantD 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 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 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 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 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 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 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 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 η (experiment V). Figure 7gh 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 κ (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. Stud-ies reported that many models greatly underestimate highlatitude dust emissions (Bullard et al., 2016;Meinander et al., 2022), and the use of u * it will mitigate this issue. Figure 7g shows that although the changes over nonarid regions are high in ratios (e.g., > 1000 times over Canada in Fig. 7h), the emission redistributions in magnitude are still small (e.g., ∼ 10 −5 kg m −2 yr −1 in Canada) compared to the major sources (e.g., > 10 −3 kg m −2 yr −1 in the Sahara) because the magnitudes in experiment III are too small over nonarid regions. The unnormalized total emissions (Fig. S6d) vastly increase to ∼ 13 200 Tg yr −1 when employing u * it instead of u * ft , more than 4 times that of experiment III.
On the other hand, the effect of multiplying the emission flux with the intermittency factor η is less dramatic than the effect of using u * it . η tends to scale down small emissions (Fig. 5b), so fluxes from major sources are only moderately reduced. In contrast, fluxes from marginal source regions (e.g., high-latitude boreal forests) are typically reduced by ∼ 1-3 orders of magnitude (blue areas in Fig. S6e). However, experiment V has a global total of ∼ 11 700 Tg yr −1 , which is only 11 % smaller than experiment IV, because those remote regions as such already have small emissions. All in all, accounting for both the impact threshold and the intermittency factor will increase the global total emission from ∼ 2880 Tg yr −1 in experiment III to ∼ 11 700 Tg yr −1 in experiment V, which is about a 4-fold increase. Figure 7h shows that C19 mainly increases marginal emissions; the overall effect of C19 is thus to move emissions from the hyper-arid regions to semiarid regions. Summing up the absolute magnitude of changes in Fig. 7g, out of 5000 Tg yr −1 there is 3163 Tg yr −1 of dust redistributed within the source regions, indicating that the intermittency scheme induces a similar magnitude of changes compared to employing F 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 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 (experiment I) to our scheme (experiment V). Figure 7i-j show that compared to K14, our scheme's emission fluxes over densely vegetated regions (e.g., equatorial Africa and northern Australia) are reduced due to the drag partition effect, while there are increases in marginal sources like the Arctic and mid-latitude boreal forests due to the intermittency effect. The figures show essentially a combination of drag partition ( Fig. 7e-f) and intermittency (Fig. 7g-h) effects. Major sources are more affected by the drag partition effect (e.g., the Bodélé Depression and El Djouf), while marginal sources are more dominated by turbulence and intermittency (e.g., the Arctic). For regions where both effects take place, more vegetated semiarid regions are more affected by the drag partition (e.g., western United States), while less vegetated semiarid regions are more affected by the intermittency (e.g., Patagonia, the Great Plains of the United States, and southern Australia). For unnormalized emissions, our scheme's global total of 11 700 Tg yr −1 is ∼ 60 % smaller than the K14 emission. A notable feature is that the new mechanisms favor the emissions over the Horn of Africa (HOA) the most, with an emission increase of ∼ 2 kg m −2 yr −1 (as seen in Figs. 7g and S7f). This is because in C19, the HOA has low u * it , high summertime u * , low roughness element cover (and thus high F eff ), and moderate soil moisture (and thus high dust emission 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.

The grid-scale dependence of our new dust emission scheme
The spatial resolution of a GCM strongly affects the budget and spatiotemporal variability of dust emissions because modeled emissions scale nonlinearly with input meteorological fields (Ridley et al., 2013;Meng et al., 2021). Many current GCMs have a horizontal resolution of ∼ 1-2 • (Zhao et al., 2022). Most of the datasets employed in this paper, such as the 0.5 • MERRA-2 fields and the even finer aeolian roughness length 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 emission has nonlinear dependences on multiple variables such as u * and w, using a simple, area-weighted spatial average of u * to calculate dust emission would be inaccurate, as it is different from an area-weighted average of high-resolution dust emissions per se, i.e., for any n > 1 u n * s < u n * s leads to F d u n * s < F d (u n * s ) (Ridley et al., 2013) (here we use the bar as a spatial average from a finer to a coarser grid). This inequality applies to our model and datasets as well: simply taking the area-weighted mean of high-resolution u * s or w in a model grid box will omit the locally high u * s (or low w) values that can produce locally extremely high emissions, resulting in an underestimation of emissions relative to direct area-weighted averaging of the emissions. Additionally, the presence of thresholds in dust emission parameterizations further intensifies the scale dependence, because spatially averaging u * s might cause u * s < u * it which leads to zero emission over a coarse grid box, whereas fine emissions could be large than zero when u * s > u * it in any fine grid. Dust emission flux will thus be strongly dependent on model resolution more than other linear processes (Zender et al., 2003a;Ridley et al., 2013;Meng et al., 2021), which is undesirable. There is a need to better upscale low-resolution dust emissions to match the variability of high-resolution emissions such that dust emissions tend to be less resolution dependent. To address the problem of missing emissions due to the smoothing of the subgrid wind maxima, a common approach is to employ a grid-by-grid Weibull distribution to the GCM winds to represent the subgrid wind maxima and thus obtain the sub- 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 highresolution 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 areaweighted mean of all input meteorological and land-surface variables to various coarser resolutions. Starting from the native 0.5 • × 0.625 • resolution of MERRA-2, we regrid fields to 0.9 • × 1.25 • , 1.9 • × 2.5 • , and 4 • × 5 • . Then we use these input fields with our new scheme to compute hourly dust emissions for 2006 across these four resolutions. We then compare the emission outputs across these resolutions.
We examine the global, regional, and grid-level scale dependence of our dust emission scheme in Fig. 8. At the regional level, Fig. 8a shows the unnormalized annual total emissions (in Tg yr −1 ) of nine major dust source regions. The source regions are defined following Kok et al. (2021a) (see Figs. S9 and 8c-d). The regional dust abundance in Fig. 8a is mainly consistent with the regional dust emissions in Kok et al. (2021b;see their Fig. 2). The highest emissions occur over the Middle East/Central Asia, followed by the Sahel and northern Africa. Smaller emissions are over East Asia, North and South America, Australia, and southern Africa. The dust emission scheme also shows scale dependence across different resolutions. Some regions may have a sharper decrease in emissions from 0.9 • × 1.25 • to 1.9 • × 2.5 • (e.g., northeastern Africa), and some regions may have a sharper decrease from 1.9 • × 2.5 • to 4 • × 5 • (e.g., Sahel). This difference is contingent upon the degree of smoothing of the input fields such as u * s at a particular resolution. The emissions will drastically drop when the local extrema of u * s are smoothed out and no longer can be represented at a particular resolution, and over different regions this cutoff may occur at different resolutions depending on the spatial heterogeneity of the local-scale meteorological fields. Nevertheless, in general, the coarser the resolution, the worse the model can represent the local variability of u * s and other input fields and subsequently the emissions, and thus the magnitudes of the emissions decrease with resolution. Figure 8a shows that the relative differences in regional emissions can be different in different resolutions (e.g., Sahel can have larger emissions than the Middle East/Central Asia in the 4 • × 5 • simulation), which will subsequently affect the spatial variability of other major dust cycle variables (such as DAOD or deposition). Therefore, it is always preferential for GCMs to simulate the dust cycle in high resolutions. Figure 8b shows the scale de-pendence of the global emissions, which also shows the same decrease in emissions from fine to coarse grid resolutions.
We also examine the scale dependence of dust emissions at the grid level. Figure 8c-d show the spatial distributions of 0.5 • × 0.625 • and 1.9 • × 2.5 • unnormalized emissions (with a global total of 11 700 and 5450 Tg yr −1 , respectively). The 0.5 • × 0.625 • simulation shows a more detailed local spatial variability compared to the coarser 1.9 • × 2.5 • simulation. The 1.9 • × 2.5 • simulation fails to capture some of the high emission regimes in the 0.5 • × 0.625 • simulation such as over the Taklamakan, and it fails to simulate the local emission peaks over marginal sources such as the Chihuahuan Desert and Patagonia. Coarse-gridded simulations lose emissions in both major and marginal dust sources. We calculate in Fig. 8e the grid-by-grid ratio of unnormalized 0.5 • × 0.625 • emissions to 1.9 • × 2.5 • emissions to show the emissions missed by the 1.9 • × 2.5 • simulation in each grid. Both the coarse (y axis) and fine (x axis) emissions are plotted in the log 10 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 lowresolution simulation can capture the largest emissions well (top righthand corner), as predicted by the high-resolution simulation. However, the smaller the emission, the more significant the difference in emissions, and for minimal emissions (< 10 −5 kg m −2 yr −1 ), the difference can be up to 3 orders of magnitude. There are very few grids with lowresolution emissions higher than high-resolution emissions. Those are exceptional cases due to the spatial variability of moisture w more heterogeneous than that of u * s , leading to the smoothing of local maxima of u * ft instead of u * s and thus the overestimation of 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 highresolution 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 factorsK c to account for the differences in the spatial variability of dust emissions between high-and low-resolution simulations: K c (long, lat) = F d,f (long, lat) /F d,c (long, lat) .
We obtain a global map of correction factorsK c and apply this map to all simulated coarse emissions from both historical and future simulations, correcting their spatial variability to match that of the high-resolution emissions. We show the seasonal variability ofK 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 ofK c . However,K c varies modestly across season because the subgrid variability of the meteorological and land-surface fields is partly determined by spatial structures such as orography or land use/land cover, which change across a relatively long enough timescale (decadal to multidecadal or longer). We note that theoretically this approach will fail to work if, over the remote regions, the low-resolution emission is zero throughout the entire simulation period while the high-resolution emission is a small positive definite, sincẽ 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 highresolution 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 • (K c,2 ) generally have bigger magnitudes than those in Fig. 9a for 0.9 • × 1.25 • (K c,1 ). These maps indicate that in coarse-gridded simulations, dust emissions are overall underrepresented over marginal sources and overrepresented over major sources. Therefore, we propose implementing these maps into GCMs of ∼ 1 • or coarser resolution to correct the dust emission spatial variability accordingly. Scaling all simulations across dif-6510 D. M. Leung et al.: A new desert dust emission scheme -Part 1: Description ferent spatial resolutions to the finest spatial resolution will move our dust emission scheme toward a scale-aware and grid-independent formulation.

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 (Dust-COMM) 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 . 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 3 Tg yr −1 for dust PM 20 . The global total of ∼ 4700 Tg yr −1 is close to the 5000 Tg yr −1 we adopted in this study for normalization, and we again normalized the DustCOMM global emissions to 5000 Tg yr −1 . Note that Kok et al. (2021a, b) only constrained the emissions and other dust variables for each broad region, but its subregional spatial distribution of dust is a multimodel mean and thus unconstrained. Therefore, in Fig. 10b, we sum up the emissions of DustCOMM to regional total emissions, which are constrained by the regional DAOD from Ridley et al. (2016). We also sum up the emissions of all other schemes to regional levels to evaluate each scheme's regional spatial variability against DustCOMM. Tables S3 and S4 summarize the global total emissions and regional emissions of Dust-COMM and all other schemes.
We first compare the global emission maps between Dust-COMM 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 highlatitude 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 Dust-COMM 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 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, Emissions of all resolutions are constrained and normalized to have a global total of 5000 Tg yr −1 before calculating the ratios. use remotely sensed aerosol optical depth, DustCOMM and MERRA-2 emissions show relatively similar spatial distributions. Figure S12a shows MERRA-2 annual mean emissions, which show very similar features of grid-by-grid variability especially over Australia, East Asia, and the Middle East. For the same reason, MERRA-2 emissions have a high correlation of r = 0.86 with DustCOMM.
Next, we aggregate the emission maps to regional total emissions to examine their regional variability. In Fig. 10b, we compare the simulated dust emissions summed over each of nine source regions against the regional DustCOMM emissions. We also summed the emissions outside all rectangular boxes as the high-latitude emissions, to yield the 10 data points in Fig. 10b. High-latitude emissions from our scheme (Fig. 7b) mainly include emissions from Alaska, Canada, Greenland, and Iceland, and there are no emissions from Antarctica because of the lack of necessary input data there (e.g., soil texture and roughness due to rocks). However, since K21 does not provide emission estimates outside of the nine defined source regions, we compare against the estimate of high-latitude dust emissions from Bullard et al. (2016) (hereafter B16) obtained from GCM results. Their definition of high-latitude emission not only includes emissions from the abovementioned regions but also from Patagonia, so we add Patagonia (39 • -56 • S of S. America) emissions as part of the high-latitude emissions in Fig. 10b. B16 estimated that high-latitude emissions (without error estimates) accounted for 4 %-5 % of their assumed 2000 Tg yr −1 of global total emissions. We normalized their estimate to match our global total of 5000 Tg yr −1 , yielding a highlatitude emission range of 200-250 Tg yr −1 . We take the average, which is 225 Tg yr −1 . For all schemes and datasets discussed here, Table S4 provides a list of regional emission contributions to the global total emission. Figure 10b shows that our scheme's emissions are in overall better agreement with DustCOMM emissions than the K14 scheme. Some of the most significant differences in emissions between our scheme and K14 are over regions including Australia, North America, and southern Africa, where the vegetation drag partitioning causes strong reductions in winds and emission fluxes (Fig. 7f) from K14 to our scheme. Our scheme's East Asian emission is significantly higher than K14's (also shown in Fig. 7j), primarily due to the switch from using u * ft to u * it in the dust emission equation (Figs. 7h and S8d). Emissions over South America, the Middle East, and the three regions of northern Africa have relatively small and negligible differences between K14 and our scheme. This occurs because both the drag partitioning and intermittency effects create only minimal changes to the emissions over these regions (Fig. 7j). The results with our scheme (blue color) better match the DustCOMM regional emissions than results with the K14 scheme, lying substantially closer to the 1 : 1 (black) line over most regions including Africa, Asia, and Australia. There are two notable exceptions where our scheme has less agreement than K14 with DustCOMM, namely North America and the high-latitude emissions. Our scheme generates fewer dust emissions over the Mojave-Sonoran-Chihuahuan deserts over the United States-Mexico border compared to the K14 emissions (Fig. 7a), because of the high LAI (annual mean > 0.4) over the western United States that leads to the strong wind drag partitioning. Meanwhile, our scheme generates significant high-latitude emissions over the Arctic, which were not captured by K14 emissions due to the very high u * ft . Using the emission intermittency parameterization, our scheme represents one of the earliest attempts to successfully capture significant levels of high-latitude emissions. Our high-latitude emissions account for 262 Tg yr −1 (without Patagonia) and 308 Tg yr −1 (with Patagonia), in total accounting for 5 %-6 % of a global sum of 5000 Tg yr −1 , which is very close to the percentage B16 suggested. Our scheme has an 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).

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 modifi-cation made in this paper. We provide discussions on further limitations for each modification in Sect. S6. For the new dry-soil median diameter D p representation, we used a globally constant D p0 = 127 ± 47 µm for calculating dust emission thresholds. In theory, 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 D p over arid regions because (1) there were no statistically significant correlations between 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 D p measurements over arid regions found 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 ρ a = 1.225 kg m −3 ). Our analysis of a compilation of D p measurements also suggests that the spatial variability of D p over deserts is relatively small compared to D p over nonarid regions. Our results thus surprisingly suggest that the D p parameterization can reasonably be much simplified by using a constant value over arid regions. We also compared our approach in deriving D p against some other derived 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 aeolian roughness. We also note that ESA CCI vegetation cover fractions are annual values and do not exhibit seasonality. The seasonal variability of the simulated 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 , σũ s , u * it , and u * ft (see Sect. S2) and is thus very sensitive to the accuracy of the GCM simulations of these four variables. For instance, if the thresholds are overestimated by the employed threshold schemes, not only will emissions be underestimated but η from C19 will also be close to zero and further worsen the low bias of the simulated dust emissions. Therefore, a prerequisite of employing the C19 scheme is that the wind friction velocity u * s and the thresholds u * it and u * ft should be adequately simulated and have reasonable ranges of variability throughout the year. We also note that C19 was designed to be used in climate models that run in the Reynolds-averaged ) based on emissions from six different models that were adjusted using inverse modeling to match constraints on particle size distribution, extinction efficiency, and regional dust aerosol optical depth. The rectangular regions specify the nine source regions defined by Kok et al. (2021a) as also shown in Fig. 8c-d. (b) DustCOMM regional emissions (based on fractional emissions reported in the fifth column of Table 1 in K21b scaled to a global total of 5000 Tg yr −1 ) versus the regional emissions computed by the K14 scheme and our new scheme. (c, d) Comparison of the DustCOMM regional emissions with the regional emissions computed by the (c) Z03-Z and Z03-G schemes and with (d) MERRA-2 dust emissions. The regional emissions of all simulations are obtained from DustCOMM following the nine source regions in panel (a), with one extra data point representing the high-latitude emissions estimated by Bullard et al. (2016). The error bars show 1 standard error, except that the B16 high-latitude emission does not include an error estimate. The black line shows the 1 : 1 line. Navier-Stokes (RANS) mode. If climate models are resolving turbulence (i.e., in the LES mode), there is minimal need to use C19 to account for intermittent emissions.
For the dust emission scaling method (Sect. 4.3), we note that the scaling factors neglect seasonal variability, which Fig. S10 indicates is moderate. However, employing an annual scaling map like Fig. 9 will already address a large part of the scale-dependence problem. Although we suggested the use of an annual scaling map, ESMs and CTMs that focus on present-day simulations can also perform multiyear simulations in both high and native (coarser) grid resolutions to obtain their own monthly climatological maps of scaling factors. Afterwards, ESMs only need to read in the climatological monthly scaling maps to rescale the native grid dust emissions every month before passing the dust emissions to the atmospheric model component. We further discussed our scaling approach versus some other existent approaches (e.g., using a Weibull distribution for the winds) in Sect. S6.4.
Finally, all emission maps produced in this paper depend on the accuracy of the representations of the input 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 final simulation of our scheme outperforms other dust emission schemes employed in this paper (Fig. 10), some features in the dust emissions map are unrealistic. First, for instance, the Australian dust emissions are of a comparable order of magnitude to the East Asian emissions (and even larger than East Asian dust emissions in coarser resolutions, per Fig. 8a), which might be because the soil moisture over Australia is slightly underestimated by MERRA-2. Second, northeastern China has larger emissions than northwestern Chinese deserts (Fig. 7b), and similar spatial variability is also seen in other past studies (e.g, Kok et al., 2014b), which might be due to the stronger friction velocity over northeastern China (annual mean > 0.3 m s −1 in MERRA-2) than over northwestern China (annual mean ∼ 0.2 m s −1 in MERRA-2) in the GCMs or input MERRA-2 fields (see Fig. S15 for the spatial distribution of MERRA-2 u * over 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 * (∼ 1 m s −1 in the summer), low soil moisture, and aeolian roughness (see Fig. 7i). However, this emission peak is not consistent with dust aerosol optical depth (DAOD) observations. There might be several reasons for the HOA emissions peak, including that (i) the input fields are biased over the HOA, (ii) some unknown mechanisms are responsible for suppressing the HOA dust, and (iii) normalized emissions outside of the HOA are overly suppressed relative to the HOA by the drag partitioning and the intermittency effect.

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 gridresolution 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 D p as a critical parameter that determines the dust emission thresholds u * it (impact) and u * ft (fluid). We found that over the arid desert regions (LAI < 1), D p can be approximated as a global constant of 127 ± 47 µm, and over nonarid regions D p increases linearly with silt and clay content. This finding indicates that past dust modeling approaches which parameterized 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 boundarylayer 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 highresolution 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 landsurface 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 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.
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. Review statement. This paper was edited by Susannah Burrows and reviewed by two anonymous referees.