Improving annual fine mineral dust representation from the surface to the column in GEOS-Chem 14.4.1

Accurate representation of mineral dust remains a challenge for global air quality or climate models due to inadequate parametrization of the emission scheme, removal mechanisms, and size distribution. While various studies have constrained aspects of dust emission fluxes and/or dust optical depth, annual mean surface dust concentrations still vary by factors of 5–10 among models. In this study, we focus on improving the annual simulation of fine dust in the GEOS-Chem chemical transport model, leveraging recent mechanistic understanding of dust source and removal, and reconciling the size differences between models and ground-based measurements. Specifically, we conduct sensitivity simulations using GEOS-Chem in its high performance configuration (GCHP) version 14.4.1 to investigate the effects of mechanism or parameter updates on annual mean concentrations. The results are evaluated by comparisons versus Deep Blue satellite-based aerosol optical depth (AOD) and AErosol RObotic NETwork (AERONET) ground-based AOD for total column abundance, and versus the Surface Particulate Matter Network (SPARTAN) for novel measurements of surface PM_2.5 dust concentrations. Reconciling modelled geometric diameter versus measured aerodynamic diameter is important for consistent comparison. The two-fold overestimation of surface fine dust in the standard model is alleviated by 39 % without degradation of total column abundance by implementing a new physics-based dust emission scheme with better spatial distribution. Further reduction by 20 % of the overestimation of surface PM_2.5 dust is achieved through reducing the mass fraction of emitted fine dust based on the brittle fragmentation theory, and explicit tracking of three additional fine mineral dust size bins with updated parametrization for below-cloud scavenging. Overall, these developments reduce the normalized mean difference against surface fine dust measurements from SPARTAN from 94 % to 35 %, while retaining comparable skill of total column abundance against satellite and ground-based AOD.

1 Introduction

Mineral dust exerts significant impacts on air quality as the most abundant aerosol component by mass globally (Kok et al., 2021b), on ecosystem health through nutrient transport and deposition such as phosphorous (Bayon et al., 2024; Swap et al., 1992) and iron (Jickells et al., 2005), and on climate through its direct scattering and absorbing of radiation and indirect modifications of cloud properties (Kok et al., 2017; Liao and Seinfeld, 1998; Mahowald et al., 2014). Despite its importance, accurate representation of long-term concentrations of mineral dust remains a challenge for global air quality or climate models due to inadequate parametrization of the emission scheme (Darmenova et al., 2009; Kok, 2011; Leung et al., 2023), removal mechanisms (Jones et al., 2022; Petroff and Zhang, 2010; Ryu and Min, 2022; Wang et al., 2014b; Zhang and Shao, 2014; Zhang et al., 2001), and size distribution (Kok et al., 2017; Mahowald et al., 2014). Observational constraints from satellite have been applied to reduce the large uncertainty of simulated mineral dust and its emissions (Mytilinaios et al., 2023; Ridley et al., 2016). However, intercomparison projects with various models still suggest large variability within a factor of 2 for the annual mean total column abundance of mineral dust, with even larger variability in surface concentrations and deposition by factors of 5–10 (Huneeus et al., 2011; Uno et al., 2006; Wu et al., 2020).

In addition to total column observations, ground-level measurements of mineral dust offer another promising opportunity to understand mechanisms affecting the accuracy of the surface concentration simulation and the variable performance from the surface to the total column in intercomparison projects. The Surface PARTiculate mAtter Network (SPARTAN, https://www.spartan-network.org/, last access: 4 February 2025) is a globally distributed monitoring network that measures the chemical components of fine particulate matter (PM_2.5), including in arid environments (Liu et al., 2024; Snider et al., 2015). These ground-based measurements of mineral dust in PM_2.5 offer new data to evaluate, understand, and improve fine dust simulation in global models.

Dust emissions play a central role in controlling the surface and total column abundance of mineral dust (Kok et al., 2014; Leung et al., 2023; Tian et al., 2021). The predicted spatial distribution of dust emissions particularly affects the downwind dust concentrations through long-range transport and deposition (Prospero, 1999). A new physics-based dust emission scheme (Leung et al., 2023) includes recent developments in the parametrization of the threshold of friction velocity for dust mobilization (Martin and Kok, 2018), combined drag partitioning effects due to rocks (Marticorena and Bergametti, 1995) and vegetation (Pierre et al., 2014a) for a better representation of exerted surface friction velocity (Leung et al., 2023), and intermittent dust mobilization due to high-frequency turbulence (Comola et al., 2019). This dust emission scheme has achieved better spatial correlations of dust column abundance against ground-based and satellite-derived dust optical depth in the Community Earth System Model version 2 (CESM2) (Leung et al., 2023, 2024). However, the effects of these new developments of dust emission scheme on the bias against ground-based measurements of surface fine dust concentrations are less well known and require further investigation.

The source and removal of dust in the size bins used in dust parametrizations can vary by orders of magnitude across the broad size range of mineral dust (Kok, 2011; Wang et al., 2014b; Zhang et al., 2001). Accounting for this size heterogeneity among dust bins could enable better representation of the global dust cycle. Prior studies have found an underestimation of coarse dust emissions and an overestimation of fine dust (Cakmur et al., 2006; Kok, 2011; Kok et al., 2017). While various studies have focused on developing the representation of coarse or super coarse dust (Kok et al., 2017; Meng et al., 2022), investigation of the effects of different emission size distributions on ambient fine dust are needed through comparison with in situ fine dust measurements. In addition, the developments and improvements of parallel computing in air quality or climate models (Eastham et al., 2018; Harris et al., 2020; Hu et al., 2018; Martin et al., 2022) offer computational capabilities to extend dust size bins with explicit treatments that could enable better representation of dust, especially with rapid variation in processes across different sizes. While the parametrization of dry deposition has been revisited and evaluated against observations (Emerson et al., 2020), below-cloud or washout scavenging has been generally limited to lumped treatments for fine and coarse aerosols in the bulk models (Jones et al., 2022; Wang et al., 2011, 2014a). Developments of the size-resolved parametrization for below-cloud (washout) scavenging (Wang et al., 2014b) are promising to improve the wet deposition of fine dust, which is especially important in distant downwind regions due to long-range transport.

Many studies have examined daily dust variability for the purpose of short-term prediction (Amato et al., 2013; Tindan et al., 2023; Yu et al., 2021). Our study focuses on a different objective of accuracy of annual mean concentrations.

In this study, we implement recent developments of a new dust emission scheme with further refinements including the clay content and wetness in the topsoil layer; reducing the dust emissions over wet, snow and vegetation covered land surfaces; while constraining the global and regional source with satellite aerosol optical depth (AOD). We revisit the size distribution of emitted dust, explicitly track mineral dust with geometric diameter less than 2 µm in four size bins, and update the parametrization for size-resolved washout scavenging. We conduct sensitivity simulations using the GEOS-Chem chemical transport model in its high performance configuration (GCHP) to investigate the effects of these developments. We focus on improving the annual fine dust representation in GCHP from the surface to the column, by comparisons against ground-level fine dust measurements, and against the ground-based and satellite-retrieved AOD over dusty regions of the Sahara, the Middle East and Asia.

2 Data sources and model description

2.1 Data sources

Ground-based AOD measurements are obtained from the Aerosol Robotic Network (AERONET) Version 3 Level 2 database with improved cloud screening (Giles et al., 2019). The median number of days with AERONET measurements is 168 d for each site. We average daily AERONET AOD to an annual mean in the year of 2018. We use satellite retrievals of AOD from the Deep Blue algorithm (Hsu et al., 2019) based on Collection 6.1 of the Moderate Resolution Imaging Spectroradiometer (MODIS) instrument aboard the satellite platforms of Terra with local overpass around 10:30 and of Aqua around 13:30, and the Version 2.0 Deep Blue aerosol global product of the Visible Infrared Imaging Radiometer Suite (VIIRS) instruments aboard the joint NASA/NOAA Suomi National Polar-orbiting Partnership (Suomi NPP) and NOAA-20 satellites with local overpass around 13:30 (Cao et al., 2014). We choose the Deep Blue aerosol product due to its optimization for the retrieval of aerosol properties over bright surfaces, which is typical over arid regions. We average daily Deep Blue aerosol data for the year 2018. Simulated AOD is coincidently sampled with available daily Deep Blue AOD. We compare simulated AOD over mainly dusty regions (defined as ${\mathrm{AOD}}_{\mathrm{Dust}}/\mathrm{AOD}>\mathrm{0.5}$ from simulations) against satellite and AERONET AOD to reduce the effects of errors in other AOD components and focus on the performance of mineral dust.

We use the Version 4.2 Level 3 gridded cloud-free tropospheric aerosol extinction profile product during daytime and nighttime of the last 15 years (2007–2021) retrieved from the Cloud–Aerosol Lidar with Orthogonal Polarization (CALIOP) on board the Cloud–Aerosol Lidar Infrared Pathfinder Satellite Observations (CALIPSO) satellite for climatological aerosol profiles (Young et al., 2018).

We use global ground-based data from the Surface Particulate Matter Network (SPARTAN; https://www.spartan-network.org/, last access: 4 February 2025) with filter-based PM_2.5 chemical composition data (Liu et al., 2024; Snider et al., 2015). Particles with aerodynamic diameter less than 2.5 µm are collected on Teflon filters using AirPhoton SS5 sampling stations with a sharp-cut cyclone (SCC) 1.829 that operates at a target flow rate of 5 L min⁻¹. The sampling station follows either a standard sampling protocol or the National Aeronautics and Space Administration (NASA) – Italian Space Agency (ASI) Multi-Angle Imager for Aerosols (MAIA) sampling protocol. Under the standard sampling protocol, PM_2.5 is collected at staggered 3-hour intervals over a 9 d period, generating a 24 h PM_2.5 sample covering a full diel cycle. Under the MAIA sampling protocol, PM_2.5 is collected continuously for 24 h from 09:00 a.m. to 09:00 a.m. at a mission-defined frequency, which has been typically every 3 d during the sampling periods used here. The starting dates for MAIA sites are listed in Table A1 in the Appendix. SPARTAN samples are analyzed for fine mineral dust concentrations using X-ray Fluorescence (XRF) and a global mineral dust equation (Eq. A1 in the Appendix; Liu et al., 2022) including correction of attenuation effects due to mass loading. The 5-year averaged surface fine dust concentrations from SPARTAN sites are listed in Table A1. We use data from sites with at least 10 samples for the 5-year (2019–2023) period after the network began using XRF with samples. A sensitivity analysis requiring at least 50 samples per site is also conducted. This study used 2296 filters from 25 SPARTAN sites for a total of 10 072 observational days.

Ground-based observations of PM_2.5 dust over North America are constructed with a global dust (Eq. A1; Liu et al., 2022) and the elemental measurements from the Air Quality System (AQS) database for speciated PM_2.5 observations in the United States (https://aqs.epa.gov/aqsweb/airdata/download_files.html#Daily, last access: 8 April 2025) and from the National Air Pollution Surveillance Program in Canada (https://data-donnees.az.ec.gc.ca/data/air/monitor/national-air-pollution-surveillance-naps-program/Data-Donnees/2018/?lang=en, last access: 8 April 2025). The AQS database includes measurements from both the Interagency Monitoring of Protected Visual Environments (IMPROVE) and Chemical Speciation Network (CSN) networks.

2.2 GEOS-Chem chemical transport model

We use the GEOS-Chem chemical transport model (http://www.geos-chem.org, last access: 4 February 2025) in its high-performance configuration (Eastham et al., 2018) version 14.4.1 (The International GEOS-Chem User Community, 2024) with improved performance and usability (Martin et al., 2022). The model is driven by meteorological inputs from GEOS Forward Processing (GEOS-FP; https://gmao.gsfc.nasa.gov/, last access: 4 February 2025) with a fine resolution 0.25° × 0.3125° (∼ 25 km) and 72 hybrid sigma-pressure vertical levels up to 0.01 hPa. GEOS-FP uses dynamic near-real-time assimilation algorithms compared to consistent static assimilation algorithms used in Modern-Era Retrospective analysis for Research and Applications Version 2 (MERRA-2; https://gmao.gsfc.nasa.gov/GMAO_products/, last access: 19 April 2025). We choose GEOS-FP over MERRA-2 for this study since GEOS-FP offers finer resolution for dust emission calculations.

GEOS-Chem simulates detailed oxidant-aerosol chemistry in the troposphere and stratosphere, with gas-phase mechanism of HO_x-NO_x-BrO_x-VOC-O₃ chemistry (Bey et al., 2001; Wang et al., 2021), coupled to aerosol chemistry for sulfate-nitrate-ammonium (SNA) aerosol (Park et al., 2004), black carbon (BC) (Wang et al., 2014a), and primary and secondary organic aerosol (Pai et al., 2020), sea salt (Jaeglé et al., 2011), and natural and anthropogenic dust (Fairlie et al., 2007; Meng et al., 2021; Philip et al., 2017; Zhang et al., 2013). The gas-aerosol partitioning for SNA is computed by the HETP v1.0 thermodynamic module (Miller et al., 2024). We use the simple, irreversible, direct yield scheme for secondary organic aerosol production (Pai et al., 2020). The effects of aerosol on photolysis rates are computed with relative humidity dependent aerosol size distributions and optical properties for hydrophilic aerosols with improved parametrization for the effective radii of inorganic and organic aerosols (Latimer and Martin, 2019; Ridley et al., 2012; Zhu et al., 2023) and updated optical properties for aspherical hydrophobic mineral dust (http://geoschemdata.wustl.edu/ExtData/CHEM_INPUTS/CLOUD_J/v2025-01/FJX_scat-aer.dat, last access: 7 April 2025) for different dust size bins as calculated by Singh et al. (2024) using the T-matrix method for an equiprobable mixture of prolate and oblate spheroids with varying aspect ratios using complex refractive indices from Sinyuk et al. (2003).

The standard dry deposition scheme in GEOS-Chem accounts for gravitational settling, aerodynamic resistance with respect to turbulent transport within the surface layer, and surface resistance to particle-surface contact due to Brownian diffusion, impaction, and interception with an observation constrained parametrization (Emerson et al., 2020; Zhang et al., 2001). Wet deposition includes separate algorithms for scavenging in convective updrafts, and in-cloud and below-cloud scavenging from precipitation (Liu et al., 2001; Wang et al., 2011, 2014a).

Emissions for GEOS-Chem are configured using the Harmonized Emissions Component (HEMCO) module v3.9.1 (Lin et al., 2021). Global anthropogenic emissions are from the Community Emissions Data System (CEDS) v2 at 0.5° × 0.5° resolution (Feng et al., 2020). Offline emissions of lightning NO_x (Murray et al., 2012), biogenic VOCs, soil NO_x, sea salt (Weng et al., 2020) and mineral dust (Sects. 2.3 and 4.2) at 0.25° × 0.3125° resolution are included to represent emission processes at the finest available resolution and to enable consistent emission fluxes across model resolutions. Open fire emissions are from the daily Global Fire Emissions Database (GFED) v4.1s (Giglio et al., 2013) at 0.25° × 0.25° resolution. Other default emission inventories in GCHP v14.4.1 include volcanic SO₂ emissions (Fisher et al., 2011), marine emissions of dimethylsulfide (DMS) (Breider et al., 2017) at 1° × 1° resolution, and ammonia at 0.25° × 0.25° resolution (Bouwman et al., 1997; Croft et al., 2016). We conduct GCHP simulations at C48 (∼ 200 km) resolution for the full year of 2018 following a one-month spin-up.

2.3 Default dust emission scheme

The default dust emission scheme in GEOS-Chem (hereafter GC Dust) originally implemented by Fairlie et al. (2007) is based on the semi-empirical Mineral Dust Entrainment and Deposition (DEAD) emission scheme (Zender et al., 2003) and the GOCART topographical source function (Ginoux et al., 2001) updated to a fine resolution of 0.25° × 0.25° (Meng et al., 2021). The total dust emission flux in kg m⁻² s⁻¹ is calculated based on Zender et al. (2003) and Fairlie et al. (2007):

$$\begin{array}{}\text{(1)}& F_{\mathrm{d}}=C_{\mathrm{g}}{C}_{\mathrm{NA}}{f}_{\mathrm{bare}}S\mathit{\phi }{Q}_{\mathrm{s}}\end{array}$$

where C_g is a global scaling factor and C_NA is a regional scaling factor in North America for total annual emissions of ∼ 2000 Tg yr⁻¹ as optimized by Meng et al. (2021); f_bare is the bare ground fraction as specified in Zender et al. (2003) to reduce dust emissions over wet, snow and vegetation covered surfaces:

$$\begin{array}{}\text{(2)}& f_{\mathrm{bare}}=(\mathrm{1}-A_{\mathrm{l}}-A_{\mathrm{wl}})\left(\mathrm{1}-A_{\mathrm{snow}}\right)\left(\mathrm{1}-{\displaystyle \frac{\mathrm{LAI}}{{\mathrm{LAI}}_{\mathrm{thr}}}}\right)\end{array}$$

where A_l, A_wl, and A_snow is the fraction of land covered by lakes, wetlands, and snow, respectively. LAI is the leaf area index, and LAI_thr is the threshold LAI to reduce the bare soil fraction due to vegetation cover, which is set to 0.3 m² m⁻² by default.

S is a topographical source function (Ginoux et al., 2001) updated at fine resolution of 0.25° × 0.25° and multiplied by the fraction of bare surface within each grid cell (Meng et al., 2021); φ is the sandblasting efficiency to convert horizontal saltation flux to vertical dust flux (Marticorena and Bergametti, 1995):

$$\begin{array}{}\text{(3)}& \mathit{\phi }={\mathrm{10}}^{\mathrm{13.4}{f}_{\mathrm{clay}}-\mathrm{4}}\end{array}$$

where f_clay is the clay content in the topsoil layer and is set to a global constant value of 0.2 to reduce excessive sensitivity of dust emission fluxes to f_clay (Zender et al., 2003). Q_s is the horizontal saltation flux as described in Sect. A2 in the Appendix.

2.4 Size distribution of emitted dust

The default size distribution of emitted dust in GEOS-Chem implemented by Zhang et al. (2013) is based on the Brittle Fragmentation Theory (Kok, 2011) with fitted parameter values for better agreement of dust observations from the Interagency Monitoring of Protected Visual Environments (IMPROVE) ground-based monitoring network in the United States:

$$\begin{array}{}\text{(4)}& {\displaystyle \frac{\mathrm{d}{V}_{\mathrm{d}}}{\mathrm{d}\ln D_{\mathrm{d}}}}={\displaystyle \frac{D_{\mathrm{d}}}{c_V}}\left[\mathrm{1}+\mathrm{erf}\left({\displaystyle \frac{\ln \left(D_{\mathrm{d}}/\stackrel{\mathrm{‾}}{D_{\mathrm{s}}}\right)}{\sqrt{\mathrm{2}}\ln {\mathit{\sigma }}_{\mathrm{s}}}}\right)\right]\exp \left[-{\left({\displaystyle \frac{D_{\mathrm{d}}}{\mathit{\lambda }}}\right)}^{\mathrm{3}}\right]\end{array}$$

where V_d is the normalized volume for emitted dust aerosols in diameter of D_d in µm; c_V is the normalization constant to make the integration total of V_d of 1; $\stackrel{\mathrm{‾}}{D_{\mathrm{s}}}=\mathrm{3.4}$ µm is the median diameter of soil particles; σ_s=3.0 is the geometric standard deviation of soil particles; λ is the side crack propagation length, whose value is 8 µm in the default particle size distribution (PSD) used in the GEOS-Chem (GC PSD), and is 12 µm in the Kok PSD (Kok, 2011).

Table 1 The binning of mineral dust in 4-bin and 7-bin simulations using GEOS-Chem. The geometric diameter range is listed in the bracket adjacent to each size bin in unit of µm.

Dust aerosols are conventionally separated into several dust bins to compromise between accuracy and computational expense (Ginoux et al., 2001; Zender et al., 2003). Table 1 summarizes the binning of mineral dust in 4-bin and 7-bin simulations. In the GEOS-Chem standard bulk configuration used here, 4 dust size bins are used including DST1 to DST4 covering geometric diameter of 0.2–12.0 µm (Fairlie et al., 2007). For DST1, 4 sub-bins of 0.2–0.36, 0.36–0.6, 0.6–1.2, and 1.2–2.0 µm are further separated for heterogeneous chemistry and AOD calculations, with shared emission, transport and deposition altogether as DST1 (Fairlie et al., 2007). To improve submicron dust representation, we implement full separation of the 7 dust bins for coupled physical and chemical processes in GEOS-Chem, as discussed in Sect. 4.3.2.

2.5 Reconciling geometric and aerodynamic diameter

A recent study has emphasized the importance of reconciling the geometric diameter used in models and the aerodynamic diameter used in ground-based measurements, especially for mineral dust with higher particle density of ∼ 2500 kg m⁻³ than the standard density of 1000 kg m⁻³ and with aspherical shapes observed in the atmosphere (Huang et al., 2021). We harmonize the differences between geometric diameter and aerodynamic diameter based on Reid et al. (2003):

$$\begin{array}{}\text{(5)}& D_{\mathrm{aer}}=D_{\mathrm{geo}}\sqrt{{\displaystyle \frac{{\mathit{\rho }}_{\mathrm{d}}}{\mathit{\chi }{\mathit{\rho }}_{\mathrm{0}}}}}\end{array}$$

where D_aer is the aerodynamic diameter; D_geo is the geometric diameter; ρ_d=2500 kg m⁻³ is the dust density; ρ₀=1000 kg m⁻³ is the standard spherical particle density; χ is the dynamic shape factor calculated by $\mathit{\chi }=\frac{\mathrm{1}}{\mathrm{2}}\left(F_{\mathrm{s}}^{\mathrm{1}/\mathrm{3}}+\frac{\mathrm{1}}{F_{\mathrm{s}}^{\mathrm{1}/\mathrm{3}}}\right)$ and F_s is Stokes form factor (Bagheri and Bonadonna, 2016; Huang et al., 2020) which can be calculated by $\mathrm{HWR}(\frac{\mathrm{1}}{\mathrm{AR}}{)}^{\mathrm{1.3}}$ where $\mathrm{AR}=\mathrm{1.70}\pm \mathrm{0.03}$ is the particle length to width ratio, and $\mathrm{HWR}=\mathrm{0.40}\pm \mathrm{0.07}$ is the particle height to width ratio (Huang et al., 2021). With this conversion, the aerodynamic diameter of 2.5 µm corresponds to the geometric diameter of 1.7 µm. The mass fraction of each simulated dust size bin to the total fine dust mass concentrations can be calculated by the integration of the dust size distribution of Eq. (4) with the λ value of 8 µm of the default PSD used in the GEOS-Chem (GC PSD), which is 68 % of DST1 with diameter of 0.2–2.0 µm.

In addition to harmonizing different size types used in models and measurements, prior studies also suggested that the sharpness of size cut-off of different inlets used to collect PM_2.5 samples can affect the measured concentrations (Kenny et al., 2000; Peters et al., 2001). To evaluate the effects, we obtain the dust size distributions of different inlets by multiplying their penetration efficiencies (Peters et al., 2001) and GC PSD (Eq. 4).

Figure 1 shows the effects of the sharpness of size cut on the size distribution of collected dust PM_2.5 samples. All four inlets have a penetration efficiency of near unity for dust with geometric diameter less than 1.0 µm, which diminishes to 0.5 at a geometric diameter of 1.7 µm and further diminishes with increasing diameter. The Well Impactor Ninety-Six (WINS) referenced by the Federal Reference Method (FRM) exhibits the sharpest size cut. The corresponding dust PSD is sharply attenuated for geometric diameters greater than 1.7 µm. The resultant effects on the mass fractions of the dust size bin to be included in dust PM_2.5 are small, with the mass fraction of DST1 ranging from 65 %–70 %. The mass fraction based on SCC 1.829 as used by SPARTAN differs by only −0.4 % from that based on the original GC PSD without inlet penetration correction. In our Base simulation using the standard version of GEOS-Chem, we calculate surface PM_2.5 dust as 67.6 % of DST1 to account for both aerodynamic diameter and inlet collection efficiency. Neglect of these effects would have increased simulated PM_2.5 dust concentrations by a factor of 2.

Figure 1 Normalized particle size distribution (PSD) used by default in GEOS-Chem (GC PSD) in solid black with left axis; penetration efficiencies for different types of PM_2.5 inlets shown in dashed colored lines with right axis, including the Well Impactor Ninety-Six (WINS), and three types of Sharp-Cut Cyclone (SCC) inlets; Solid colored lines show the adjusted GC PSD collected by different inlets. Grey dash-dotted line indicates the corresponding geometric diameter of 1.7 µm for the aerodynamic diameter of 2.5 µm. Filled rectangles indicate size ranges of 4 dust size bins. Percentages adjacent to GC PSD and different inlets are mass fractions of DST1 for the calculation of PM_2.5 dust concentrations.

3 Strong overestimation of surface fine dust

Figure 2 shows the spatial distributions of the annual total column AOD and surface PM_2.5 dust from AERONET, SPARTAN, and the Base simulation using the standard version of GEOS-Chem in the year of 2018. Mineral dust largely determines the AOD in AERONET and GEOS-Chem over and downwind of the main dust source regions including the Sahara, Middle East, and the Taklamakan and Gobi deserts in Asia. The simulated AOD over dusty regions (simulated ${\mathrm{AOD}}_{\mathrm{Dust}}/\mathrm{AOD}>\mathrm{0.5}$) exhibits a high degree of consistency versus the ground-based observations of AERONET AOD with the regression slope near unity and R² of 0.7. However, the simulated surface PM_2.5 dust exhibits a pronounced overestimation by a factor of 2.4 compared to the ground-based measurements of SPARTAN. Simulated PM_2.5 dust is overestimated at the dusty sites of Abu Dhabi in the United Arab Emirates by 163 %, Ilorin in Nigeria by 108 %, and Kanpur in India by 96 %.

Figure 2 Annual simulated aerosol optical depth (AOD) and comparison against ground-based observations from AERONET over dusty regions (simulated ${\mathrm{AOD}}_{\mathrm{Dust}}/\mathrm{AOD}>\mathrm{0.5}$) (top) in the year of 2018; Annual simulated surface PM_2.5 dust and comparison against ground-based measurements from SPARTAN (bottom) from the Base simulation in the year of 2018. Filled circles on the maps represent ground-based observations from SPARTAN and AERONET. Inset values at the bottom right of the maps are arithmetical mean with 5th and 95th percentiles in the square brackets. Regression statistics including reduced-major-axis linear regression equation, coefficient of variation (R²), total number of points (N), normalized mean difference (NMD), and normalized root-mean-square difference (NRMSD) are listed at the bottom right of the scatter plots. Major source regions over land are outlined in red including: (1) the Sahara – SA, (2) Middle East – ME, and (3) Asia – AS. Major dust outflow regions over ocean are outlined in green including: (4) the Caribbean Sea – CRB, (5) the tropical Atlantic Ocean – TAT, (6) the Mediterranean Sea – MED, (7) the Arabian Sea – ARB, (8) the tropical Indian Ocean and the Bay of Bengal – IND, and (9) the northwestern Pacific Ocean – NWP.

Figure 3 shows the vertical profile of the aerosol extinction normalized by AOD over the main dust source regions and associated downwind regions, to understand the significant performance difference between the surface and the column, with the absolute extinction profile shown in Fig. A1 in the Appendix. The simulated vertical profile exhibits overall agreement against the 15-year (2007 to 2021) climatological mean extinction vertical profile from the CALIOP, with no evidence of a model overestimate in the lower mixed layer versus aloft, indicating the vertical distribution of mineral dust is not the main driver of the performance discrepancy between the surface and the column. However, further evaluations of the vertical profile near the surface are needed as CALIOP retrievals are challenging at lower altitudes especially below 100 m.

Figure 3 Comparisons of the annual extinction vertical profile normalized by total column aerosol optical depth from the Base simulation in the year of 2018 against the 15-year (2007 to 2021) climatological mean extinction vertical profile from the CALIOP over different regions including the major dust source regions over land of the Sahara – SA, Middle East – ME, and Asia – AS, and the major dust outflow regions over ocean of the Caribbean Sea – CRB, the tropical Atlantic Ocean – TAT, the Mediterranean Sea – MED, the Arabian Sea – ARB, the tropical Indian Ocean and the Bay of Bengal – IND, and the northwestern Pacific Ocean – NWP.

4 Model revisions to reduce the overestimation of surface fine mineral dust

To reduce the overestimation of surface PM_2.5 dust, we (1) implement a new dust emission scheme with further refinements for soil properties including the clay content and soil wetness in the top soil layer and the threshold of leaf area index, (2) revisit the size distribution of emitted dust, (3) explicitly track dust with geometric diameter less than 2 µm in four size bins, and (4) update the parametrization for size-resolved below-cloud scavenging.

4.1 Sensitivity simulation setup

Figure 4 summarizes the setup of sensitivity simulations to evaluate the effects of algorithmic modifications and their performance versus satellite-retrieved AOD and surface dust measurements. The default dust simulation (Base) in GEOS-Chem as implemented by Fairlie et al. (2007) uses the DEAD emission scheme (Zender et al., 2003) with a topographical source function (Ginoux et al., 2001; Meng et al., 2021) for natural dust (GC Dust) with 4 dust sizebins for emission, transport and removal with 7 dust size bins for dust optical depth calculation and heterogeneous chemistry. To improve the spatial distributions of dust total column abundance, we implement a new dust emission scheme developed by Leung et al. (2023) (DustL23; Emis). Additional modifications on top of the original DustL23 emission scheme include (1) reducing the sensitivity of soil clay content by eliminating the multiplication of the factor of the capped soil clay content $f_{\mathrm{clay}}^{\prime }$ (EmisClay); (2) halving the topmost soil wetness in the layer of 0–5 cm to approximate the soil wetness in the top 1–2 cm layer which is most pertinent to dust emissions (Darmenova et al., 2009; Wu et al., 2022) (EmisClayWet); and (3) reducing the threshold of LAI_thr from 1.0 to 0.5 m² m⁻² (EmisClayWetLAI_thr or Emis^*). To further improve the surface fine dust simulation, we update the GEOS-Chem particle size distribution (PSD) with the PSD developed by Kok (2011) (Emis^*PSD) with a larger value for the side crack propagation length of λ (12 µm versus 8 µm) which reduced the mass fraction of emitted fine dust. The Kok PSD was shown to have excellent agreement versus various soil size measurements (Kok, 2011), especially for fine dust distributions (González-Flórez et al., 2023). Lastly, we allow for the four dust bins with geometric diameter less than 2 µm to have separate emission, transport, and dry and wet deposition while halving anthropogenic dust emissions from AFCID (Emis^*PSD7Bins0.5AD), and with updated below-cloud or washout scavenging parametrization (Emis^*PSD7Bins0.5ADWetDep). Each of these changes is examined below.

Figure 4 Sensitivity simulation setup. The grey box indicates default settings with the default dust emission scheme used in GEOS-Chem (GC Dust) with 4 dust size bins (Base). The orange box indicates the implementation of a modified dust scheme based on DustL23 (Emis^*). Modifications based on the original DustL23 scheme with the soil texture dataset from the Global Soil Dataset for use in Earth System Models (GSDE) (Emis) include the soil clay content (EmisClay), soil wetness (EmisClayWet), and threshold leaf area index (EmisClayWetLAI_thr). The simulation setup for EmisClayWetLAI_thr is the same as that for Emis^*. The blue box indicates the modification of size distribution of emitted dust (Emis^*PSD). The green boxes indicate the improvements for fine dust including explicit tracking of dust with diameter less than 2 µm with a total of 7 dust size bins with halved anthropogenic fugitive, combustion, and industrial dust (AFCID) emissions (Emis^*PSD7Bins0.5AD), and updating below-cloud (washout) scavenging coefficients (Emis^*PSD7Bins0.5ADWetDep).

The total global annual source strength for each sensitivity simulation is scaled to achieve unity slope versus Deep Blue AOD (Fig. A2) over major dust source regions. The surface PM_2.5 dust concentrations are calculated by accounting for aerodynamic diameter and inlet penetration efficiency (Section 2.5) as 0.676 DST1 for 4-bin simulations, and DSTbin1 + DSTbin2 + DSTbin3 + 0.546 DSTbin4 for 7-bin simulations. We focus our evaluation on the skill in representing in situ PM_2.5 dust concentrations measured by SPARTAN, and in representing the spatial variation in annual mean AOD. Regression equations are calculated using reduced-major-axis linear regression (Smith, 2009) to account for uncertainties in both simulations and measurements.

4.2 Improving the spatial distribution of mineral dust with updated emission scheme

We implement into GEOS-Chem a new physics-based dust emission scheme developed by Leung et al. (2023) (DustL23) to replace the default dust emission scheme (Sect. 2.3) used in GEOS-Chem (GC Dust). The spatial distributions of DustL23 in the Community Earth System Model version 2 (CESM2) exhibited better correlation against dust optical depth datasets and AERONET AOD than the DEAD scheme (Leung et al., 2024). We modify DustL23 for implementation into GEOS-Chem by (1) reducing dust emissions over wet, snow, and vegetation covered surfaces of semi-arid regions using Eq. (7) below, (2) eliminating the multiplication of the capped clay content of the topsoil in Eq. (8) below, (3) halving the soil wetness in the layer of 0–5 cm to represent the soil wetness in the top 1–2 cm layer which is most pertinent to dust emissions (Darmenova et al., 2009; Wu et al., 2022), (4) applying a regional scaling factor of 0.6 over the Sahara to reduce its emissions (Eq. 8), and (5) scaling the global total emission flux to achieve unity regression slope of simulated AOD versus Deep Blue AOD over dusty regions.

We begin with the formulation for total dust emission flux F_d in kg m⁻² s⁻¹ following Leung et al. (2024):

$$\begin{array}{}\text{(6)}& \begin{array}{rl}& F_{\mathrm{d}}=\mathit{\eta }{C}_{\mathrm{tune}}{C}_d{f}_{\mathrm{bare}}{f}_{\mathrm{clay}}^{\prime }{\displaystyle \frac{{\mathit{\rho }}_{\mathrm{a}}\left(u_{\ast \mathrm{s}}^{\mathrm{2}}-u_{\ast \mathrm{it}}^{\mathrm{2}}\right)}{u_{\ast \mathrm{st}}}}{\left({\displaystyle \frac{u_{\ast \mathrm{s}}}{u_{\ast \mathrm{it}}}}\right)}^{\mathit{\kappa }}\\ & \mathrm{for}{u}_{\ast \mathrm{s}}>u_{\ast \mathrm{it}}\end{array}\end{array}$$

where η is an intermittency factor, C_tune is a global tuning factor for the emission strength, C_d is the time-varying soil erodibility coefficient, f_bare is the bare ground fraction, $f_{\mathrm{clay}}^{\prime }$ is the clay content in the topmost soil layer of f_clay capped at 0.2, ρ_a is the surface air density in kg m⁻³, u_∗s is the soil surface friction velocity in m s⁻¹ corrected from the surface friction velocity of u_∗ by the drag partitioning effects of F_eff, u_∗it is the dynamic or impact threshold friction velocity in m s⁻¹, u_∗st is the standardized wet fluid threshold friction velocity in m s⁻¹, and κ is the fragmentation exponent. We use u_∗st in the denominator of Eq. (6) following Kok et al. (2014) instead of u_∗it following Leung et al. (2023) for tuning purpose. The parametrization details for these factors following Leung et al. (2023) can be found in the Appendix, Sect. A3.

We modify the DustL23 scheme (Leung et al., 2023) by adopting the equation for the bare ground fraction in Zender et al. (2003) to reduce dust emissions over wet, snow and vegetation covered surfaces with the dry erodible land fraction taken from satellite-based land cover:

$$\begin{array}{}\text{(7)}& f_{\mathrm{bare}}=A_{\mathrm{erod}}\left(\mathrm{1}-A_{\mathrm{snow}}\right)\left(\mathrm{1}-{\displaystyle \frac{\mathrm{LAI}}{{\mathrm{LAI}}_{\mathrm{thr}}}}\right)\end{array}$$

where A_erod is the area fraction of erodible surfaces including barren and sparsely vegetated land cover taken from the MODIS Land Cover Climate Modeling Grid (CMG) (MCD12C1) Version 6.1 data product; A_snow is the area fraction of snow cover, LAI is the leaf area index (Yuan et al., 2011), and LAI_thr is the threshold LAI to reduce the bare soil fraction due to vegetation cover. We set an intermediate value of LAI_thr=0.5 m² m⁻² instead of 1.0 m² m⁻² in Leung et al. (2023) to represent the reduction in dust emissions from sparse vegetation over semi-arid regions, which is more similar to the value of 0.3 used in prior work (Mahowald et al., 1999; Zender et al., 2003).

The enhancement factor f_m≥1 for the wet fluid threshold friction velocity due to soil wetness is calculated using Eqs. (A8) and (A9), but with spatially varying clay content f_clay in the topsoil layer. The gridded f_clay dataset is taken from the Global Soil Dataset for use in Earth System Models (GSDE) with various inputs from global and regional soil database (Shangguan et al., 2014), rather than the machine-learning trained Soil Grids v2.0 dataset with very few observations over arid regions (Poggio et al., 2021) used in Leung et al. (2023). In addition, we reduce the effects of clay content on dust emissions by eliminating the multiplication of the capped clay content $f_{\mathrm{clay}}^{\prime }$. Soil wetness is taken from the parent meteorological inputs of GEOS-FP (Koster et al., 2020) which targets the top 5 cm layer that desiccates more slowly following precipitation than the soil wetness in the top 1–2 cm layer (Swenson and Lawrence, 2014) that is most pertinent to dust emissions; we halve the soil wetness in an attempt to represent this process (Darmenova et al., 2009; Wu et al., 2022).

The global scaling factor C_tune is determined by the reduced-major-axis linear regression slope of simulated AOD versus satellite-retrieved AOD over dusty regions ($\frac{{\mathrm{AOD}}_{\mathrm{Dust}}}{\mathrm{AOD}}>\mathrm{0.5}$) in this study to constrain the intensity of dust emissions, whose values corresponding to different emission schemes are listed in Table A2. Additionally, a regional scaling factor of 0.6 over the Sahara (C_sah) and unity elsewhere is applied to reduce regionally excessive dust emissions that may be influenced by the tendency for global models to overrepresent emissions from large source regions compared with smaller sources (Kok et al., 2021a; Zhao et al., 2022).

The final formulation for dust emission flux is:

$$\begin{array}{}\text{(8)}& \begin{array}{rl}& F_{\mathrm{d}}=\mathit{\eta }{C}_{\mathrm{sah}}{C}_{\mathrm{tune}}{C}_d{f}_{\mathrm{bare}}{\displaystyle \frac{{\mathit{\rho }}_{\mathrm{a}}\left(u_{\ast \mathrm{s}}^{\mathrm{2}}-u_{\ast \mathrm{it}}^{\mathrm{2}}\right)}{u_{\ast \mathrm{st}}}}{\left({\displaystyle \frac{u_{\ast \mathrm{s}}}{u_{\ast \mathrm{it}}}}\right)}^{\mathit{\kappa }}\\ & \mathrm{for}{u}_{\ast \mathrm{s}}>u_{\ast \mathrm{it}}\end{array}\end{array}$$

The calculated offline hourly dust emissions at 0.25° × 0.3125° resolution using Eq. (8) are then used to drive GCHP simulations at C48 resolution. The spatial distributions predicted from different emission schemes are evaluated against satellite-based Deep Blue AOD, ground-based AERONET AOD, and SPARTAN surface PM_2.5 dust measurements.

Figure 5 shows the spatial distributions of annual dust emission fluxes and dust optical depth predicted from different emission schemes, with Fig. 6 showing the comparisons against Deep Blue satellite AOD globally and over major dust source regions. Comparison of the Base and Emis schemes reveals that the latter captures more secondary dust emission spots, especially over the Sahara, and inland dust sources in Saudi Arabia. However, the comparison against Deep Blue AOD over the Sahara is degraded versus the default scheme (Fig. 6). As suggested by prior studies, soil clay content is an important factor affecting the threshold friction velocity (Fécan et al., 1999; Tian et al., 2021; Zender et al., 2003) and sandblasting efficiency (Zender et al., 2003), and is often tuned for the optimization of dust emissions (Leung et al., 2024; Tian et al., 2021). Eliminating the multiplication of the capped clay content of $f_{\mathrm{clay}}^{\prime }$ reduces the effects of the clay content, increasing emissions from the Bodélé Depression in Chad and El Djouf across the border of Mauritania and Mali over the Sahara, from the Rub' al Khali desert in the inland Saudi Arabi, and Taklamakan desert in the northwest China (Fig. 5, EmisClay). Correspondingly, the R² from the linear regression against Deep Blue AOD is improved from 0.60 to 0.70 over the Sahara, from 0.68 to 0.77 over the Middle East, and from 0.35 to 0.56 over Asia (Fig. 6). The other two modifications of halving soil wetness (EmisClayWet) and setting LAI_thr to 0.5 m² m⁻² (EmisClayWetLAI_thr) slightly improve the spatial distribution of dust emissions by reducing the underestimation in Asia while retaining the agreements in the Sahara and Middle East (Fig. 6). Using the same dusty region of the Base (Fig. A3) or EmisClayWetLAI_thr (Fig. A4) scheme for the comparisons of all dust emission schemes versus Deep Blue AOD confirms similarly slight improvements of regional dust emissions. Together these refinements exhibit comparable global performance as the Base simulation versus Deep Blue AOD with improvements to the relative regional magnitude of dust across the Sahara, Middle East and Asia as indicated by more comparable regression slopes (Fig. 6).

Figure 5 Annual dust emission flux (left) and simulated dust optical depth (AOD_Dust; right) in the year of 2018 zoomed in over dusty regions of the Sahara, Middle East, and Asia from different emission schemes as described in Fig. 4. Inset values are the regional arithmetical mean with 5th and 95th percentiles in the square brackets.

Figure 6 Comparisons of annual simulated aerosol optical depth (AOD) versus the Deep Blue satellite AOD globally (GL) and over main dust source regions of the Sahara – SA, Middle East – ME, and Asia (AS) with different emission schemes. Regression statistics including reduced-major-axis linear regression equation, coefficient of variation (R²), total number of points (N), normalized mean difference (NMD), and normalized root-mean-square difference (NRMSD) are in the top left. Note the total number of points varies across different schemes.

Figure 7 shows the evaluation of the Emis^* (or EmisClayWetLAI_thr) simulation with ground-based observations from AERONET and SPARTAN. The overestimation of surface PM_2.5 dust against the ground-based measurements of SPARTAN is reduced from 94 % (Fig. 2) to 55 % (Fig. 7), reflecting regional improvements of the spatial distributions especially over the Middle East (Fig. 6). The skill in representing AOD in the Emis^* simulation remains comparable to that in the Base simulation shown in Fig. 2.

Figure 7 Annual simulated aerosol optical depth (AOD) and comparison against ground-based observations from the AERONET over dusty regions (${\mathrm{AOD}}_{\mathrm{Dust}}/\mathrm{AOD}>\mathrm{0.5}$) (top); Annual simulated surface PM_2.5 dust and comparison against ground-based measurements from the SPARTAN from the Emis^* simulation in the year of 2018 (bottom). Filled circles on the maps represent ground-based observations from SPARTAN and AERONET. Inset values at the bottom right of the maps are arithmetical mean with 5th and 95th percentiles in the square brackets. Regression statistics including the reduced-major-axis linear regression equation, coefficient of variation (R²), total number of points (N), normalized mean difference (NMD), and normalized root-mean-square difference (NRMSD) are listed at the bottom right of the scatter plots.

4.3 Improving the representation of fine mineral dust

As the size distribution of mineral dust is particularly important for the performance discrepancy between simulated AOD over dusty regions and surface PM_2.5 dust, we focus on improving its size-resolved source and sink.

4.3.1 Revisiting the size distribution of emitted mineral dust

Figure 8a shows different PSDs including the default PSD used in the GEOS-Chem (GC PSD) based on the brittle fragmentation theory with the side crack propagation length λ of 8 µm (Zhang et al., 2013), the Kok PSD with λ of 12 µm (Kok, 2011), and the Meng PSD focusing on the optimization for coarse to super coarse dust (Meng et al., 2022), in comparison with the observed PSD from the 2011 Fennec campaign (Ryder et al., 2013). While all modelled PSDs are within the wide range of PSD from the Fennec campaign, the fraction of emitted DST1 from the Kok PSD exhibits greater consistency with the Fennec observations than the other two PSDs. The larger discrepancy for the size distribution with diameter less than ∼  0.4 µm between the observed PSD from Fennec and parametrized PSDs is possibly due to anthropogenic aerosol influence (González-Flórez et al., 2023). In addition, a recent field study in the Moroccan Sahara (González-Flórez et al., 2023) indicated overall agreement of emitted dust size distributions against the Kok PSD especially at the fine diameter range. Therefore, we adopt the Kok PSD with λ of 12 µm for the size distribution of emitted mineral dust in GEOS-Chem. Figure 8b shows the spatial distribution from the Emis^*PSD simulation which remains similar to that from the Emis^* simulation in Fig. 7. Reduced emissions from DST1 by using the Kok PSD reduces the overestimation of surface PM_2.5 dust from 55 % to 33 % compared to the ground-based measurements from SPARTAN (Fig. 8c).

Figure 8 (a) Normalized particle size distribution (PSD) of emitted dust based on default PSD used in GEOS-Chem (GC PSD) (Zhang et al., 2013), the Kok PSD (Kok, 2011), the Meng PSD (Meng et al., 2022), and the Fennec PSD (Ryder et al., 2013). All PSDs are normalized for a total volumetric integration of 1 within the diameter range of 0.2 to 12 µm used in GEOS-Chem. The grey shades show the minimum and maximum PSD curves from the Fennec 2011 campaign. The grey dash-dotted line indicates the corresponding geometric diameter of 1.7 µm for the aerodynamic diameter of 2.5 µm. Filled rectangles indicate size ranges of 4 dust size bins. Percentages adjacent to each PSD are mass fractions of emitted DST1 over total dust emission flux within diameter range of 0.2 to 12 µm. (b) Simulated annual surface PM_2.5 dust from the Emis^*PSD simulation in the year of 2018. Filled circles on the map represent ground-based observations from SPARTAN and AERONET. Inset values at the bottom right of the maps are arithmetical mean with 5th and 95th percentiles in the square brackets. (c) Comparison of simulated PM_2.5 dust versus observed fine dust from SPARTAN. Regression statistics including the reduced-major-axis linear regression equation, coefficient of variation (R²), total number of points (N), normalized mean difference (NMD), and normalized root-mean-square difference (NRMSD) are listed at the bottom right.

4.3.2 Improving the size-resolved dry and wet deposition of mineral dust

The default below-cloud (washout) scavenging of dust by rain and snow in GEOS-Chem is separated for fine (DST1) and coarse dust (DST2 to DST4) (Wang et al., 2011). However, washout scavenging coefficients strongly depend on aerosol size (Wang et al., 2014b). To improve the size-dependent washout treatment of dust, we update washout rates by rain and snow for 7 dust size bins by (Wang et al., 2014b):

$$\begin{array}{}\text{(9)}& \mathrm{\Lambda }=A\left(D_{\mathrm{d}}\right){\left({\displaystyle \frac{P_d}{f_r}}\right)}^{B\left(D_{\mathrm{d}}\right)}\end{array}$$

where Λ is the washout scavenging coefficient in s⁻¹ by either rain or snow; P_d is the precipitation rate in mm h⁻¹ falling form upper layers; f_r is the area fraction of precipitation within each grid box; A and B are empirical constants dependent on dust size D_d. Using the same semi-empirical equations for A and B as Wang et al. (2014b), the updated values for different dust size bins with different effective spherical radii are summarized in Table 2.

Table 2 Values of A and B for washout parametrizations by rain and snow for different dust size bins.

Figure 9 shows the size-dependent variations of mineral dust dry and wet deposition. We adopt the effective radii for 7 dust size bins for dry and wet deposition. The dry deposition velocity can vary by a factor of 4.9 among Bin1 to Bin4 with the minimum near the geometric diameter of 0.5 µm. The washout scavenging coefficient can vary by a factor of 2.6 among Bin1 to Bin4 with the minimum near the geometric diameter of 0.4 µm. Given the steep increasing strength of emitted dust from Bin1 to Bin4 (Fig. 8), there is need to explicitly track dust within DST1. We evaluate these developments by examining their effects on the fractional contributions of fine dust to total dust.

Figure 9 Size-resolved dry deposition velocity over desert (left y-axis) and washout scavenging coefficient by rain (right y-axis). Dry deposition velocity is calculated with the friction velocity of 0.4 m s⁻¹ and the particle density of 2500 kg m⁻³ with the default dry deposition scheme used in the GEOS-Chem. Washout scavenging coefficient is calculated with the precipitation rate of 0.1 mm h⁻¹ with the updated washout parametrization. Solid lines indicate theoretical parametrization. Dashed lines indicate the default discrete treatment. Dotted lines indicate the updated discrete treatment. Grey dash-dotted line indicates the corresponding geometric diameter of 1.7 µm for the aerodynamic diameter of 2.5 µm.

Figure 10 shows the fractional contributions of fine dust with geometric diameter less than 2 µm to total dust (${\mathrm{AOD}}_{\mathrm{FineDust}}/{\mathrm{AOD}}_{\mathrm{Dust}}$) from the simulations with a total of 7 dust bins for dry deposition with updated washout scavenging parametrization and their differences. Due to the dominance of dry deposition over arid dusty regions, the explicit tracking of fine dust dry deposition slightly reduces ${\mathrm{AOD}}_{\mathrm{FineDust}}/{\mathrm{AOD}}_{\mathrm{Dust}}$ over major dust source regions. However, the anthropogenic contributions to fine dust are correspondingly enhanced over urban and industrial regions, leading to degraded comparison against SPARTAN measurements (Fig. A5). We thus halve the AFCID emissions to reduce the excessive contributions from this uncertain source (Emis^*PSD7Bins0.5AD). In addition, accounting for the steep washout scavenging efficiency across DSTbin5 to DSTbin7 (Fig. 9) with updated washout parametrization would induce enhanced fractional contributions especially for DSTbin5 (Fig. A6) and thus relatively reduce fractional contributions from fine dust with geometric diameter less than 2 µm to total dust (${\mathrm{AOD}}_{\mathrm{FineDust}}/{\mathrm{AOD}}_{\mathrm{Dust}}$). Figure 11 shows the overall performance with all revisions from the simulation of Emis^*PSD7Bins0.5ADWetDep. The reduced-major-axis linear regression slope is further reduced from 1.68 (Fig. 8) to 1.59 with comparable values of NMD against SPARTAN measurements.

Figure 10 Fractional contributions of fine dust with geometric diameter less than 2 µm to total dust column abundance (${\mathrm{AOD}}_{\mathrm{FineDust}}/{\mathrm{AOD}}_{\mathrm{Dust}})$ from the (a) Emis^*PSD7Bins0.5ADWetDep, (b) Emis^*PSD7Bins0.5AD, (c) Emis^*PSD and their absolute differences. Inset values at the bottom right are arithmetical mean with 5th and 95th percentiles in the square brackets.

Figure 11 Annual simulated aerosol optical depth (AOD) and comparison against ground-based observations from AERONET over dusty regions (${\mathrm{AOD}}_{\mathrm{Dust}}/\mathrm{AOD}>\mathrm{0.5}$) (top); Annual simulated surface PM_2.5 dust and comparison against ground-based measurements from SPARTAN from the Emis^*PSD7Bins0.5ADWetDep simulation in the year of 2018 (bottom). Filled circles on the maps represent ground-based observations from SPARTAN and AERONET. Inset values at the bottom right of the maps are arithmetical mean with 5th and 95th percentiles in the square brackets. Regression statistics including the reduced-major-axis linear regression equation, R², total number of points (N), normalized mean difference (NMD), and normalized root-mean-square difference (NRMSD) are listed at the bottom right of the scatter plots.

Comparisons against other surface dust datasets also show improved or comparable performance compared to the Base simulation. Figure A7 shows the comparison against ground observations over North America. Using the refined new dust emission scheme with the replacement of the size distribution from the Kok PSD, explicitly tracking submicron bins for dry deposition, and updating the washout scavenging parametrization contribute to a comparable extent to reduce the overestimation over North America from 43 % of the Base simulation to 15 % of the Emis^*PSD7Bins0.5ADWetDep simulation. Comparisons against surface concentrations and total deposition of PM₁₀ dust (Li et al., 2022b) for the Emis^*PSD7Bins0.5ADWetDep simulation are also comparable with the Base simulation (Figs. A8 and A9). Consistent with prior studies about the spatial sensitivity of dust emissions (Leung et al., 2023; Meng et al., 2021), fine-resolution meteorological fields are needed to capture dust emission hotspots. If the dust emissions were calculated with C48 meteorological fields, the global dust distribution would become more concentrated in the major global source regions with the elimination of marginal dust sources, and the R² versus SPARTAN surface PM_2.5 dust would diminish to 0.83 (Fig. A10). Overall comparisons for the seasonal mean between the Base and the Emis^*PSD7Bins0.5ADWetDep simulations confirm largely reduced overestimation for the surface fine dust against SPARTAN, while retaining comparable skill for the total column AOD against AERONET (Figs. A11 to A14). The reduction of surface overestimation is especially prominent over dusty seasons in Spring (from 73 % to 48 %) and Summer (from 138 % to 50 %), while further improvements are needed for surface overestimation in Fall (from 140 % to 95 %).

Table 3 summarizes the effects of different modifications on the model performance of total column AOD and surface fine mineral dust in this study. Strong overestimation of surface PM_2.5 dust concentrations exist in the Base simulation by a factor of 2.4 versus SPARTAN measured dust. Updating the dust emission scheme with further refinements in the soil properties reduces the overestimation of surface PM_2.5 dust by 39 %. The surface overestimation by 55 % is reduced to 35 % by updating the size distribution of emitted dust, explicitly tracking dust with diameter less than 2 µm in 4 bins, and updating the parametrization of below-cloud scavenging. The comparisons of simulated AOD versus AERONET and Deep Blue AOD are comparable for all simulations with the correlation coefficient of 0.8–0.9, and NMDs from −9 % to 31 %. The emissions between the Base and Emis^* simulations are comparable with the global annual dust emission of ∼ 2000 Tg yr⁻¹, which is within the range of 1000–5000 Tg yr⁻¹ from intercomparison projects (Huneeus et al., 2011; Wu et al., 2020). As the Kok PSD reduces the mass fraction of fine dust, the total emitted mass is enhanced to ∼ 3000 Tg yr⁻¹ with larger contributions from coarse dust. The reduction of surface PM_2.5 dust overestimation with these revisions is confirmed if SPARTAN sites are restricted to those with at least 50 samples as well (Table A3).

Table 3 Effects of different modifications on the model performance of simulated annual surface PM_2.5 dust versus SPARTAN, and simulated annual aerosol optical depth (AOD) versus AERONET AOD and Deep Blue satellite AOD in terms of the correlation coefficient (r), the reduced-major-axis linear regression slope, and the normalized mean difference (NMD), with associated annual dust optical depth (AOD_Dust) and total dust emissions in the year of 2018.

5 Conclusions

In summary, we evaluate and improve the annual mineral dust simulation in the GEOS-Chem model by building upon recent ground-based measurements from SPARTAN of mineral dust in PM_2.5 over land, together with total column AOD from AERONET measurements and from the MODIS and VIIRS Deep Blue satellite products. We devote attention to the representation of aerodynamic diameter when comparing with ground-based PM_2.5 measurements, since representation as geometric diameter in models would introduce a two-fold bias. We nonetheless find that the standard GEOS-Chem chemical transport model much better represents columnar AOD with a slope near unity than surface PM_2.5 dust concentrations which are overestimated by a factor of two. Comparison of simulated extinction profiles versus the 15-year climatological CALIOP extinction profiles yields overall consistency in the vertical shape (Fig. 3), indicating the importance of other dominant factors.

We develop the mineral dust representation in GEOS-Chem with attention to its sources, size distribution, and sinks. We implement a new dust emission scheme based on Leung et al. (2023) with further refinements to the clay content and wetness in the topsoil layer, threshold leaf area index, and reducing dust emissions over snow and vegetation covered land surfaces. The NMD versus surface measurements is reduced by 39 % while the simulated AOD better represents the spatial distribution of Deep Blue AOD over dusty regions. To further improve the fine dust representation in GEOS-Chem, we revisit the size distribution of emitted dust and find that the Kok particle size distribution (PSD; Kok, 2011) better represents the mass fraction of fine dust measured during the Fennec field campaign over Northern Africa than does the default PSD despite the uncertainties from the Fennec observations. The implementation of the Kok PSD into GEOS-Chem reduces the surface overestimation of PM_2.5 dust by 22 %. We also enable explicit tracking of mineral dust with geometric diameter less than 2 µm in 4 size bins for emission, transport, and deposition with updated parametrization for below-cloud scavenging, which further reduces the overestimation of surface PM_2.5 dust concentrations to within 35 %.

Despite these advances, challenges remain in mineral dust development and evaluation. The performance of AOD against satellite and AERONET observations over dusty regions may still be affected by other aerosol components which may benefit from further evaluations and developments. Although the simulations are only for a single year, we average the multi-year observational data from the CALIOP extinction profile and SPARTAN measured surface dust concentrations for long-term representativeness. This approach benefits from the weak interannual variability of annual mean mineral dust concentrations (Li et al., 2017; Song et al., 2021). Nonetheless, additional observational data will enable further evaluation of the performance of mineral dust simulations. In addition, knowledge gaps remain for mechanistic representation of mineral dust emissions. We call for further developments on the parametrization of dust emissions, particularly for the uncertainties in global and regional dust emission strength and further constraints on the effects of soil wetness on the threshold friction velocity. Future examination of daily variability would also be valuable for short-term predictability.

These investigations indicate the importance of size type reconciliation in models versus measurements, the spatial distribution of dust emissions, the size distribution of emitted dust, and the explicit tracking of fine dust bins for more accurate simulation of fine dust abundance from the surface to the column.

Appendix A: Additional details about dust emission parametrizations, SPARTAN dust, and complementary figures

A1 A global dust equation

We follow a global dust equation for the calculation of surface PM_2.5 dust concentrations from SPARTAN (Liu et al., 2022):

$$\begin{array}{}\text{(A1)}& \begin{array}{rl}\mathrm{Dust}& =\left[\mathrm{1.89}\mathrm{Al}\times \left(\mathrm{1}+\mathrm{MAL}\right)+\mathrm{2.14}\mathrm{Si}+\mathrm{1.40}\mathrm{Ca}\right.\\ & \left.+\mathrm{1.36}\mathrm{Fe}+\mathrm{1.67}\mathrm{Ti}\right]\times \mathrm{CF}\end{array}\end{array}$$

where 1.89, 2.14, 1.40, 1.36, and 1.67 are the mass conversion ratios for corresponding mineral oxides; MAL is the mineral-to-aluminum mass ratio of $({\mathrm{K}}_{\mathrm{2}}\mathrm{O}+\mathrm{MgO}+{\mathrm{Na}}_{\mathrm{2}}\mathrm{O})/{\mathrm{Al}}_{\mathrm{2}}{\mathrm{O}}_{\mathrm{3}}$; CF is a correction factor (CF) to account for other missing compounds.

A2 Horizontal saltation flux in standard version of GEOS-Chem

The default horizontal saltation flux Q_s in GEOS-Chem is based on the parametrization of White (1979):

$$\begin{array}{}\text{(A2)}& Q_{\mathrm{s}}=C_z{\displaystyle \frac{{\mathit{\rho }}_{\mathrm{a}}}{g}}{u}_{\ast \mathrm{s}}^{\mathrm{3}}\left(\mathrm{1}-{\displaystyle \frac{u_{\ast \mathrm{ft}}}{u_{\ast \mathrm{s}}}}\right){\left(\mathrm{1}+{\displaystyle \frac{u_{\ast \mathrm{ft}}}{u_{\ast \mathrm{s}}}}\right)}^{\mathrm{2}}\mathrm{for}{u}_{\ast \mathrm{s}}>u_{\ast \mathrm{ft}}\end{array}$$

where C_z=2.61 is the saltation constant; ρ_a is the air density in kg m⁻³; g=9.81 m s⁻² is the gravitational acceleration; the drag partitioning effects are ignored by default and thus $u_{\ast \mathrm{s}}=u_{\ast }$, where u_∗ is calculated from the wind speed at 10 m u_10 m based on the logarithmic wind profile within the boundary layer under adiabatic conditions (Marticorena and Bergametti, 1995):

$$\begin{array}{}\text{(A3)}& u_{\ast }={\displaystyle \frac{k{u}_{\mathrm{10}\mathrm{m}}}{\ln \left(z_{\mathrm{0}}/z_{\mathrm{0}\mathrm{a}}\right)}}\end{array}$$

where k=0.4 is the von Kármán constant; u_10 m is the wind speed at 10 m; z₀=10 m is the reference height; $z_{\mathrm{0}\mathrm{a}}={\mathrm{10}}^{-\mathrm{4}}$ m is the surface roughness height. The wet fluid threshold friction velocity of u_∗ft is the minimum surface friction velocity required to initiate the saltation from the bare soil (Fécan et al., 1999):

$$\begin{array}{}\text{(A4)}& u_{\ast \mathrm{ft}}=u_{\ast \mathrm{ft}\mathrm{0}}\cdot f_m\end{array}$$

where u_∗ft0 is the dry fluid threshold friction velocity following Iversen and White (1982):

$$\begin{array}{}\text{(A5)}& u_{\ast \mathrm{ft}\mathrm{0}}=\left\{\begin{array}{l}\frac{\mathrm{0.129}\mathrm{K}}{\sqrt{\mathrm{1.928}{\mathit{\text{Re}}}^{\mathrm{0.092}}-\mathrm{1}}},\\ \mathrm{0.03}<\mathit{\text{Re}}<\mathrm{10}\\ \\ \mathrm{0.12}\mathrm{K}[\mathrm{1}-\mathrm{0.0858}{e}^{-\mathrm{0.0617}\left(\mathit{\text{Re}}-\mathrm{10}\right)}],\\ \mathit{\text{Re}}\ge \mathrm{10}\end{array}\right.\end{array}$$

where:

$$\begin{array}{}\text{(A6)}& {\displaystyle }K=\sqrt{{\displaystyle \frac{{\mathit{\rho }}_{\mathrm{p}}g{D}_{\mathrm{p}}}{{\mathit{\rho }}_{\mathrm{a}}}}\left(\mathrm{1}+{\displaystyle \frac{\mathrm{0.006}}{{\mathit{\rho }}_{\mathrm{p}}g{D}_{\mathrm{p}}^{\mathrm{2.5}}}}\right)}\text{(A7)}& {\displaystyle };;;LETEX-DUMMY;;;\mathit{\text{Re}}=\mathrm{1331}{D}_{\mathrm{p}}^{\mathrm{1.56}}+\mathrm{0.38}\end{array}$$

Where D_p=75 µm is the diameter of soil particle which corresponds to the minimum dry fluid threshold velocity of u_∗ft0 (Iversen and White, 1982).

The enhancement factor f_m≥1 is a function of soil wetness (Fécan et al., 1999):

$$\begin{array}{}\text{(A8)}& f_m=\left\{\begin{array}{ll}\mathrm{1},& w\le w_{\mathrm{t}}\\ \sqrt{\mathrm{1}+\mathrm{1.21}{\left[\mathrm{100}\left(w-w_{\mathrm{t}}\right)\right]}^{\mathrm{0.68}}},& w>w_{\mathrm{t}}\end{array}\right.\end{array}$$

where w is the gravimetric soil moisture (kg kg⁻¹) in the shallowest soil layer; w_t is the threshold gravimetric water content above which u_∗ft increases with soil wetness (Fécan et al., 1999):

$$\begin{array}{}\text{(A9)}& w_{\mathrm{t}}=\mathrm{0.01}a\left(\mathrm{17}{f}_{\mathrm{clay}}+\mathrm{14}{f}_{\mathrm{clay}}^{\mathrm{2}}\right)\end{array}$$

where a is a tuning factor which is taken as $\mathrm{1}/f_{\mathrm{clay}}=\mathrm{5}$ by default.

A3 Additional details about the new dust emission scheme

The variables used in the calculation for the total dust emission flux F_d (Eq. 6) can be categorized into meteorological fields including η, ρ_a, and u_∗, land surface properties including f_bare, $f_{\mathrm{clay}}^{\prime }$, F_eff, and u_∗it, intrinsic soil erodibility properties including u_∗st, C_d, and κ, and a global tuning factor of C_tune.

Intermittency effects due to the fluctuation of instantaneous soil friction velocity ${\stackrel{\mathrm{̃}}{u}}_{\mathrm{s}}$ are reflected in the intermittency factor of η, which is denoted by the temporal fraction of active dust emission ranging from 0 to 1 within a transport time step. The parametrization of η is based on Comola et al. (2019):

$$\begin{array}{}\text{(A10)}& \mathit{\eta }=\mathrm{1}-P_{\mathrm{ft}}+\mathit{\alpha }\left(P_{\mathrm{ft}}-P_{\mathrm{it}}\right)\end{array}$$

where P_ft and P_it are the cumulative probability of instantaneous friction velocity larger than a wet fluid threshold, and an impact threshold, respectively; α is the fraction of ${\stackrel{\mathrm{̃}}{u}}_{\mathrm{s}}$ crossing a wet fluid threshold over the total fraction crossing a wet fluid threshold and an impact threshold.

The calculation of η is based on velocity at the saltation height of z_sal=0.1 m. Thus the surface friction velocity of u_∗s, and threshold velocities of u_∗ft and u_∗it are first calculated at the saltation height based on (Marticorena and Bergametti, 1995):

$$\begin{array}{}\text{(A11)}& u_X\left(\mathrm{sal}\right)={\displaystyle \frac{u_{\ast X}}{k}}\ln \left({\displaystyle \frac{z_{\mathrm{sal}}}{z_{\mathrm{0}\mathrm{a}}}}\right)\end{array}$$

where the subscript X can be ft, it or s, $z_{\mathrm{0}\mathrm{a}}={\mathrm{10}}^{-\mathrm{4}}$ m, and k=0.386 is the von Kármán constant.

Assuming a normal distribution of instantaneous soil friction velocity ${\stackrel{\mathrm{̃}}{u}}_{\mathrm{s}}\sim N(u_{\mathrm{s}},{\mathit{\sigma }}_{{\stackrel{\mathrm{̃}}{u}}_{\mathrm{s}}}^{\mathrm{2}})$, a standard deviation of instantaneous friction velocity ${\mathit{\sigma }}_{{\stackrel{\mathrm{̃}}{u}}_{\mathrm{s}}}$ is a central parameter to calculate the fraction of active dust emissions within a time step for transportation. ${\mathit{\sigma }}_{{\stackrel{\mathrm{̃}}{u}}_{\mathrm{s}}}$ is calculated based on the similarity theory (Panofsky et al., 1977):

$$\begin{array}{}\text{(A12)}& {\mathit{\sigma }}_{{\stackrel{\mathrm{̃}}{u}}_{\mathrm{s}}}=u_{\ast \mathrm{s}}{\left(\mathrm{12}-\mathrm{0.5}{\displaystyle \frac{z_i}{L}}\right)}^{\mathrm{1}/\mathrm{3}}\end{array}$$

where z_i is the planetary boundary layer height, and L is the Monin-bukhov length calculated by (Panofsky et al., 1977):

$$\begin{array}{}\text{(A13)}& L=-{\displaystyle \frac{{\mathit{\rho }}_{\mathrm{a}}{c}_{p}T{u}_{\ast }^{\mathrm{3}}}{kgH}}\end{array}$$

where c_p=1005 J kg⁻¹ K⁻¹ is the specific hear capacity of air under constant pressure; T is surface air temperature; u_∗ in m s⁻¹ is the original surface friction velocity without the drag partitioning correction; g=9.81 m s⁻² is the gravitational acceleration; H is the sensible heat flux from turbulence in W m⁻².

Given that a normal distribution is assumed, cumulative probabilities of P_ft and P_it can be calculated by $P_{\mathrm{ft}}=\mathrm{0.5}[\mathrm{1}+\mathrm{erf}(\frac{u_{\mathrm{ft}}-u_{\mathrm{s}}}{\sqrt{\mathrm{2}}{\mathit{\sigma }}_{{\stackrel{\mathrm{̃}}{u}}_{\mathrm{s}}}}\left)\right]$, and $P_{\mathrm{it}}=\mathrm{0.5}[\mathrm{1}+\mathrm{erf}(\frac{u_{\mathrm{it}}-u_{\mathrm{s}}}{\sqrt{\mathrm{2}}{\mathit{\sigma }}_{{\stackrel{\mathrm{̃}}{u}}_{\mathrm{s}}}}\left)\right]$. α is the number of crossing rate of ${\stackrel{\mathrm{̃}}{u}}_{\mathrm{s}}$ across the wet fluid threshold C_ft over the total number of crossing rate of ${\stackrel{\mathrm{̃}}{u}}_{\mathrm{s}}$ across the wet fluid threshold C_ft and the impact threshold C_it (Comola et al., 2019):

$$\begin{array}{}\text{(A14)}& \mathit{\alpha }={\displaystyle \frac{C_{\mathrm{ft}}}{C_{\mathrm{ft}}+C_{\mathrm{it}}}}\end{array}$$

The crossing fraction of α is approximated by $\mathit{\alpha }\approx {\left[\exp \left(\frac{u_{\mathrm{ft}}^{\mathrm{2}}-u_{\mathrm{it}}^{\mathrm{2}}-\mathrm{2}{u}_{\mathrm{s}}\left(u_{\mathrm{ft}}-u_{\mathrm{it}}\right)}{\mathrm{2}{\mathit{\sigma }}_{{\stackrel{\mathrm{̃}}{u}}_{\mathrm{s}}}^{\mathrm{2}}}\right)+\mathrm{1}\right]}^{-\mathrm{1}}$ as suggested by Comola et al. (2019).

The soil surface friction velocity of u_∗s is calculated by (Leung et al., 2023; Marticorena and Bergametti, 1995; Webb et al., 2020):

$$\begin{array}{}\text{(A15)}& u_{\ast \mathrm{s}}=u_{\ast }{F}_{\mathrm{eff}}\end{array}$$

where u_∗ is the surface friction velocity taken directly from the parent meteorological fields; F_eff is the drag partitioning effects due to the presence of non-erodible elements including rocks and vegetation.

Drag partitioning effects are calculated following Leung et al. (2023):

$$\begin{array}{}\text{(A16)}& F_{\mathrm{eff}}={\left(A_{\mathrm{r}}{f}_{\mathrm{eff},\mathrm{r}}^{\mathrm{3}}+A_{\mathrm{v}}{f}_{\mathrm{eff},\mathrm{v}}^{\mathrm{3}}\right)}^{\mathrm{1}/\mathrm{3}}\end{array}$$

where A_r is the fraction of barren and sparsely vegetated land cover approximated by A_erod; A_v is the fraction of short vegetation land cover taken from the MCD12C1 Version 6.1 land cover product; f_eff,r is the drag partitioning effects due to rocks (Marticorena and Bergametti, 1995):

$$\begin{array}{}\text{(A17)}& f_{\mathrm{eff},\mathrm{r}}=\mathrm{1}-{\displaystyle \frac{\ln \left(\frac{z_{\mathrm{0}\mathrm{a}}}{z_{\mathrm{0}\mathrm{s}}}\right)}{\ln \left[b_{\mathrm{1}}{\left(\frac{X}{z_{\mathrm{0}\mathrm{s}}}\right)}^{b_{\mathrm{2}}}\right]}}\end{array}$$

where z_0a is the aeolian roughness length which the surface roughness of overlaying nonerodable elements and was taken as the minimum of monthly mean gridded aeolian roughness length (Prigent et al., 2005); $z_{\mathrm{0}\mathrm{s}}=\frac{D_{\mathrm{p}}}{\mathrm{15}}$ is the smooth roughness length which quantifies the roughness of a bed of fine soil particles in the absence of roughness elements (Pierre et al., 2014b); b₁=0.7, b₂=0.8, and X=10 m are empirical constants (Leung et al., 2023). f_eff,v is the drag partitioning effects due to vegetation (Pierre et al., 2014a):

$$\begin{array}{}\text{(A18)}& f_{\mathrm{eff},\mathrm{v}}={\displaystyle \frac{K+f_{\mathrm{0}}c}{K+c}}\end{array}$$

where f₀=0.32 and c=4.8 are empirical constants (Okin, 2008); K is calculated by $\frac{\mathit{\pi }}{\mathrm{2}}\left(\frac{\mathrm{1}}{\mathrm{LAI}/{\mathrm{LAI}}_{\mathrm{thr}}}-\mathrm{1}\right)$ (Leung et al., 2023; Okin, 2008).

The wet fluid threshold velocity u_∗ft is calculated using Eq. (A4), except the dry fluid threshold velocity u_∗ft0 is calculated by (Shao and Lu, 2000):

$$\begin{array}{}\text{(A19)}& u_{\ast \mathrm{ft}\mathrm{0}}=\sqrt{A\left({\mathit{\rho }}_{\mathrm{p}}g{D}_{\mathrm{p}}+\mathit{\gamma }/D_{\mathrm{p}}\right)}/{\mathit{\rho }}_{\mathrm{a}}\end{array}$$

where A=0.0123 and $\mathit{\gamma }=\mathrm{1.65}\times {\mathrm{10}}^{-\mathrm{4}}$ kg s⁻² are empirical constants (Darmenova et al., 2009; Leung et al., 2023); $D_{\mathrm{p}}=\mathrm{127}\pm \mathrm{47}$ µm is the median diameter of soil particle as evaluated from various field measurements in Leung et al. (2023).

Once the saltation is initialized, the threshold velocity required to maintain the saltation diminishes, which is defined as the dynamic or impact threshold friction velocity u_∗it in m s⁻¹ (Martin and Kok, 2018):

$$\begin{array}{}\text{(A20)}& u_{\ast \mathrm{it}}=B_{\mathrm{it}}{u}_{\ast \mathrm{ft}\mathrm{0}}\end{array}$$

where B_it=0.82. A prior study suggested that the impact threshold primarily governed the saltation flux (Martin and Kok, 2018) and thus u_∗it is adopted as the governing threshold in Eq. (14).

The standardized wet fluid threshold friction velocity u_∗st was proposed and argued as a central factor to characterize soil aridity by a prior study (Kok et al., 2014):

$$\begin{array}{}\text{(A21)}& u_{\ast \mathrm{st}}=u_{\ast \mathrm{ft}}\sqrt{{\mathit{\rho }}_{\mathrm{a}}/{\mathit{\rho }}_{\mathrm{a}\mathrm{0}}}\end{array}$$

where ρ_a0=1.225 kg m⁻³ is the standard surface air density.

The fragmentation exponent of κ quantifies the sensitivity of F_d to u_∗s and is capped at 3 to prevent excessive sensitivity of the model to wind speeds according to (Kok et al., 2014; Leung et al., 2024):

$$\begin{array}{}\text{(A22)}& \mathit{\kappa }=C_{\mathit{\kappa }}{\displaystyle \frac{\left(u_{\ast \mathrm{st}}-u_{\ast \mathrm{st}\mathrm{0}}\right)}{u_{\ast \mathrm{st}\mathrm{0}}}}\end{array}$$

where $C_{\mathit{\kappa }}=\mathrm{2.7}\pm \mathrm{1.0}$ and $u_{\ast \mathrm{st}\mathrm{0}}=\mathrm{0.16}$ m s⁻¹ are constants.

The time-varying soil erodibility coefficient is a function of u_∗st only (Kok et al., 2014):

$$\begin{array}{}\text{(A23)}& C_d=C_{d\mathrm{0}}\exp \left(-C_e{\displaystyle \frac{u_{\ast \mathrm{st}}-u_{\ast \mathrm{st}\mathrm{0}}}{u_{\ast \mathrm{st}\mathrm{0}}}}\right)\end{array}$$

where $C_{d\mathrm{0}}=(\mathrm{4.4}\pm \mathrm{0.5})\times {\mathrm{10}}^{-\mathrm{5}}$ and $C_e=\mathrm{2.0}\pm \mathrm{0.3}$ are empirical constants.

Table A1 The mean and standard deviation (σ) of surface PM_2.5 dust measured from 25 SPARTAN sites with at least 10 samples in 5 years from 2019 to 2023 globally. Sites are sorted by the mean surface PM_2.5 dust concentrations.

^* The number of days when SPARTAN sampling occurred for a part of the day according to either the standard 9 d protocol or the MAIA protocol.

Table A2 The values of a global tuning factor C_tune used for different simulations.

Table A3 Effects of different modifications on the model performance of simulated annual surface PM_2.5 dust versus SPARTAN over sites with >50 samples in terms of the correlation coefficient (r), the reduced-major-axis linear regression slope, and the normalized mean difference (NMD).

Table A4 Computational demand of 4-bin and 7-bin dust simulations.

^* Calculated on Intel^® Xeon^® Gold 6154 with the clock speed of 3.00 GHz.

Figure A1 Same as Fig. 3 but for the absolute extinction vertical profile.

Figure A2 Annual aerosol optical depth (AOD) from the Deep Blue satellite retrieval and comparison against ground-based observations from AERONET in the year of 2018. Filled circles on the map represent ground-based observations from AERONET. Inset values at the bottom right of the map are arithmetical mean with 5th and 95th percentiles in the square brackets. Regression statistics including the reduced-major-axis linear regression equation, coefficient of variation (R²), total number of points (N), normalized mean difference (NMD), and normalized root-mean-square difference (NRMSD) are listed at the top left of the scatter plot.

Figure A3 Same as Fig. 6 but over the same dust source regions for the Base scheme for all dust emission scheme comparisons versus Deep Blue AOD.

Figure A4 Same as Fig. 6 but over the same dust source regions for the EmisClayWetLAI_thr scheme for all dust emission scheme comparisons versus Deep Blue AOD.

Figure A5 Differences of the fractional contributions of fine dust with geometric diameter less than 2 µm to total dust column abundance (${\mathrm{AOD}}_{\mathrm{FineDust}}/{\mathrm{AOD}}_{\mathrm{Dust}})$ between the Emis^*PSD7Bins and Emis^*PSD simulations (left); Comparison between simulated PM_2.5 dust against SPARTAN measurements from the Emis^*PSD7Bins simulation with color coded by the differences of ${\mathrm{AOD}}_{\mathrm{FineDust}}/{\mathrm{AOD}}_{\mathrm{Dust}}$ between the Emis^*PSD7Bins and Emis^*PSD simulations over SPARTAN sites. Open circles in the map indicate SPARTAN sites. Inset values at the bottom right of the map are arithmetical mean with 5th and 95th percentiles in the square brackets. Regression statistics including the reduced-major-axis linear regression equation, coefficient of variation (R²), total number of points (N), normalized mean difference (NMD), and normalized root-mean-square difference (NRMSD) are listed at the bottom right of the scatter plot.

Figure A6 Fractional contributions of DSTbin5 to total dust column abundance (${\mathrm{AOD}}_{\mathrm{DSTbin}\mathrm{5}}/{\mathrm{AOD}}_{\mathrm{Dust}})$ from the (a) Emis^*PSD7Bins0.5ADWetDep, (b) Emis^*PSD7Bins0.5AD, (c) Emis^*PSD and their absolute differences. Inset values at the bottom right are arithmetical mean with 5th and 95th percentiles in the square brackets.

Figure A7 Comparisons of simulated annual surface PM_2.5 dust against ground-based observations in the year of 2018 over North America from the Base (top), Emis^*PSD (second), Emis^*PSD7Bins0.5AD (third), and Emis^*PSD7Bins0.5ADWetDep (bottom) simulations. Filled circles represent ground-based observations of surface PM_2.5 dust concentrations. Inset values at the bottom left are arithmetical mean with 5th and 95th percentiles in the square brackets. Regression statistics including the reduced-major axis linear regression equation, coefficient of variation (R²), total number of points (N), normalized mean difference (NMD), and normalized root-mean-square difference (NRMSD) are listed at the top left of right panels.

Figure A8 Annual simulated surface PM₁₀ dust concentrations in the year of 2018 from the simulations of (a) Base, (b) Emis^*, (c) Emis^*PSD, and (d) Emis^*PSD7Bins0.5ADWetDep. Filled circles represent ground-based observations of surface PM₁₀ dust concentrations. Inset values at the bottom right are arithmetical mean with 5th and 95th percentiles in the square brackets. Dash lines in the scatter plots indicate variations within a factor of 5. Regression statistics including the reduced-major-axis linear regression equation, coefficient of variation (R²), total number of points (N), normalized mean difference (NMD), and normalized root-mean-square difference (NRMSD) are listed at the top left of right panels.

Figure A9 Annual simulated total deposition of PM₁₀ dust within the troposphere in the year of 2018 from the simulations of (a) Base, (b) Emis^*, (c) Emis^*PSD, and (d) Emis^*PSD7Bins0.5ADWetDep. Filled circles represent ground-based observations of surface PM₁₀ dust deposition. Inset values at the bottom right are arithmetical mean with 5th and 95th percentiles in the square brackets. Dash lines in the scatter plots indicate variations within a factor of 5. Regression statistics including the reduced-major-axis linear regression equation, coefficient of variation (R²), total number of points (N), normalized mean difference (NMD), and normalized root-mean-square difference (NRMSD) are listed at the top left of right panels.

Figure A10 Annual simulated aerosol optical depth (AOD) and comparison against ground-based observations from AERONET over dusty regions (${\mathrm{AOD}}_{\mathrm{Dust}}/\mathrm{AOD}>\mathrm{0.5}$) (top); Annual simulated surface PM_2.5 dust and comparison against ground-based measurements from SPARTAN from the Emis^*PSD7Bins0.5ADWetDep simulation with the dust emissions calculated at C48 resolution in the year of 2018 (bottom). Filled circles on the maps represent ground-based observations from SPARTAN and AERONET. Inset values at the bottom right of the maps are arithmetical mean with 5th and 95th percentiles in the square brackets. Regression statistics including the reduced-major-axis linear regression equation, coefficient of variation (R²), total number of points (N), normalized mean difference (NMD), and normalized root-mean-square difference (NRMSD) are listed at the top left of the scatter plots.

Figure A11 Simulated seasonal mean (March, April, and May or MAM) aerosol optical depth (AOD; a, c) and surface PM_2.5 dust (b, d) from the Base and Emis^*PSD7Bins0.5ADWetDep simulations. Filled circles on the maps represent ground-based observations from SPARTAN and AERONET. Inset values at the bottom right of the maps are arithmetical mean with 5th and 95th percentiles in the square brackets. Comparisons of simulated AOD versus AERONET AOD over dusty sites (${\mathrm{AOD}}_{\mathrm{Dust}}/\mathrm{AOD}>\mathrm{0.5}$), and simulated surface PM_2.5 dust versus SPARTAN observations are shown in the right panels. Regression statistics including the reduced-major-axis linear regression equation, coefficient of variation (R²), total number of points (N), normalized mean difference (NMD), and normalized root-mean-square difference (NRMSD) are listed at the bottom right of the scatter plots.

Figure A12 Same as Fig. A11 but for the seasonal mean of June, July, and August (JJA).

Figure A13 Same as Fig. A11 but for the seasonal mean of September, October, and November (SON).

Figure A14 Same as Fig. A11 but for the seasonal mean of December, January, and February (DJF).