Importance of different parameterization changes for the updated dust cycle modeling in the Community Atmosphere Model (version 6.1)

The Community Atmosphere Model (CAM6.1), the atmospheric component of the Community Earth System Model (CESM; version 2.1), simulates the life cycle (emission, transport, and deposition) of mineral dust and its interactions with physio-chemical components to quantify the impacts of dust on climate and the Earth system. The accuracy of such quantifications relies on how well dust-related processes are represented in the model. Here we update the parameterizations for the dust module, including those on the dust emission scheme, the aerosol dry deposition scheme, the size distribution of transported dust, and the treatment of dust particle shape. Multiple simulations were undertaken to evaluate the model performance against diverse observations, and to understand how each update alters the modeled dust cycle and the simulated dust direct radiative effect. The model–observation comparisons suggest that substantially improved model representations of the dust cycle are achieved primarily through the new more physically-based dust emission scheme. In comparison, the other modifications induced small changes to the modeled dust cycle and model–observation comparisons, except the size distribution of dust in the coarse mode, which can be even more influential than that of replacing the dust emission scheme. We highlight which changes introduced here are important for which regions, shedding light on further dust model developments required for more accurately estimating interactions between dust and climate.

1 Introduction

Mineral dust accounts for most aerosol mass in the Earth's atmosphere and plays an important role in different aspects of the coupled Earth–human–climate system. For example, dust modifies the radiative budget and atmospheric dynamics via direct, semi-direct, and indirect interactions with radiation (Sokolik and Toon, 1996; Miller and Tegen, 1999; Pérez et al., 2006; Li and Sokolik, 2018a) and clouds (DeMott et al., 2003; Rosenfeld et al., 2001; Shi and Liu, 2019). In addition, the deposition of mineral dust perturbs the energy budget by darkening snow and glacial ice sheets directly due to the relatively darker color of dust particles (Skiles et al., 2018; Sarangi et al., 2020) and indirectly by providing nutrients (e.g., phosphorus) to snow algae (McCutcheon et al., 2021). Dust deposited onto land and ocean can also affect the biogeochemistry by adding nutrients (iron and phosphorus) and/or pollutants to ecosystems (Martin et al., 1990; Swap et al., 1992; Shinn et al., 2000; Tie and Cao, 2009; Mahowald, 2011; Mahowald et al., 2017, 2010; Hamilton et al., 2020).

To quantify the climate and biogeochemical impacts of dust, accurately reproducing the dust cycle (e.g., emission, transport, and deposition) with models is required. However, previous studies have shown substantial differences between the modeled dust cycle and observations (e.g., surface dust concentration and dust deposition) (Albani et al., 2014; C. Wu et al., 2020). These uncertainties in the dust cycle modeling, as well as uncertainties in optical properties due primarily to dust size and mineral composition, suggest a large uncertainty in estimating the dust direct radiative effect (DRE) (Kok et al., 2017; Li et al., 2021).

The difficulty in modeling dust results primarily from a limited understanding of the processes that control the emission, aging, and removal of dust during transport (Sokolik et al., 2001). Past studies have documented a nonlinear response of dust emission to the soil surface state and meteorological fields (Kok et al., 2012), strong regional variation of the erodible soil composition (Claquin et al., 1999; Journet et al., 2014), complex chemical and physical aging of dust during transport (Cwiertny et al., 2008; Usher et al., 2003) at varied time and spatial scales, a wide range of dust particle size (Mahowald et al., 2014), and irregular shape of dust aerosol particles (E. Reid et al., 2003; Wang et al., 2015). These complexities impose a great challenge to parameterizing dust-related processes (e.g., dust emissions and dust deposition) and thus to accurately simulating the dust cycle in climate models. In addition, in situ or station-based measurements of dust aerosols are highly limited at both temporal and spatial scales, which makes representation of those measurements challenging, especially considering the episodic character of dust events (Mahowald et al., 2009). As such, the modeling community is still moving toward better parametrizing the different phases of the dust cycle.

To account for regional variations in dust composition and the resultant dust optical properties in estimating the dust DRE, several common and radiatively important minerals found in dust from major dust sources were introduced to the Community Atmosphere Model versions 4 (CAM4) and 5 (CAM5) (Scanza et al., 2015) and migrated to CAM6.1 (Li et al., 2021), which are the atmospheric components of the Community Earth System Model (CESM: version 1 and 2, respectively). Including the ability to resolve dust speciation along with the addition of an atmospheric iron cycle module (Scanza et al., 2018; Hamilton et al., 2019) facilitates the study of dust impacts on biogeochemical cycles (Hamilton et al., 2020).

As one of the widely used climate models, the CAM contains several weaknesses in modeling the dust cycle, as follows for example:

1. The default scheme in CAM6.1 (Zender et al., 2003a; Dust Entrainment And Deposition DEAD model, referred to as DEAD) relies on an empirical geomorphic dust source function, created based on satellite retrievals of dust source regions, to model dust emissions.

2. The current default CESM2.1 is using the dry deposition scheme Zhang et al. (2001; Z01 hereafter) developed for particle deposition over smooth and non-vegetated surfaces. This scheme, however, underemphasizes the interception loss, the mechanism of which is less influential over the other surfaces such as grassland. The use of the Z01 in the current default CESM2.1 is thus very likely overestimating the dry deposition velocity of fine-sized aerosols (diameter < 1.0 µm; referring to the geometric diameter herein unless stated otherwise) and slightly underestimating that of aerosols with diameter > 5.0 µm (Wu et al., 2018), especially over non-vegetated surfaces (Petroff and Zhang, 2010).

3. One of the changes from CAM5 to CAM6.1 was replacing the size distribution of aerosols in the coarse mode in CAM5 with the one that has a much narrower width in CAM6.1 (Table 1). This change was to accommodate stratospheric aerosols in the coarse mode (e.g., volcanic sulfate) compared to an early officially released version of this model (Mills et al., 2016). A recent model evaluation against satellite retrievals (M. Wu et al., 2020) suggests that CESM2.1-CAM6.1 worsened the dust cycle representation and stands out in simulating the relative importance of wet to dry deposition, compared with the other global climate models or model versions, such as CESM1-CAM5, due partially to the narrow coarse-mode geometric standard deviation.

4. Dust aerosols are typically aspherical particles in shape. The dust asphericity could lengthen the dust lifetime by ∼ 20 %, compared to modeling dust as spherical particles (Huang et al., 2020). Still, CAM6.1 simulates dust as spherical particles, though the impact of dust asphericity on optical depth and resulting DRE of dust (Kok et al., 2017) has been previously introduced to CAM6.1 (Li et al., 2021).

Correspondingly, this paper describes several updates to the dust representation in CAM6.1 on the four aspects and evaluates whether and for what conditions they improve the dust model comparison to observations in the present climate. Specifically, we do the following:

1. We replace DEAD with a new more physically based dust emission scheme, Kok et al. (2014a; referred to as BRIFT) previously developed for the climate models within the framework of DEAD. This scheme performs well against observations in CESM-CAM4 (Kok et al., 2014b) without the aid of the empirical geomorphic dust source function.

2. We replace Z01 with the dry deposition scheme developed by Petroff and Zhang (2010) (PZ10 hereafter) to mediate the overestimation of the dry deposition velocity of fine-sized aerosols.

3. We revert size distribution of dust aerosol particles in the coarse mode to the one previously employed in CAM5.

4. We account for the lifetime effect of dust asphericity by decreasing the modeled gravitational settling velocity.

These updates are based on up-to-date knowledge of the dust cycle and are thus more physically realistic than the default dust parameterizations in CAM6.1/Community Land Model (version 5; CLM5).

Table 1 Mode parameters for the Modal Aerosol Module version 4 (MAM4) used in CAM5 (CAM5 size) and CAM6.1 (CAM6 size) by default: geometric standard deviations (σ) and initialized geometric mean diameter (GMD) and its ranges. Values in parentheses if present are for CAM6.1; cells without parentheses are kept the same between CAM5 and CAM6.1.

We organize the paper as follows: Sect. 2 describes the model (Sect. 2.1–2.3), the modifications we made to the model (Sect. 2.5), and the experiment we conducted (Sect. 2.6) under present climate conditions to achieve our purpose. Section 3 presents the observation and semi-observation for model evaluation in current climate. Section 4 describes metrics used to assess the model performance. Section 5 evaluates the performance of the updated model by comparing simulated dust properties (e.g., surface dust concentrations, deposition fluxes, vertical distribution, and size distribution of transported dust) against measurements, retrievals, and model–observation integration (Sect. 5.1), quantifies the influence of each modification on those simulated dust properties (Sect. 5.2), documents the influence of those modifications on the estimate of the dust DRE (Sect. 5.3), and compares these changes in order to recommend which are the most important for other models to consider (Sect. 5.4). Section 6 shows the difference between the bulk- (consisting of single “bulk” composition, so no composition distinguished between particles) and speciated-dust (consisting of mineral components) models on the dust cycle modeling and the resultant dust climatic effects. Furthermore, we discuss limitations in the model–observation comparison in Sect. 7, and discussions and conclusions in Sect. 8.

2 Model descriptions

We used CAM6.1 (Sect. 2.1), embedded within the National Center for Atmospheric Research CESM2.1, to simulate the dust cycle in all the numerical experiments. This section describes bulk- (Sect. 2.2) and speciated-dust model (Sect. 2.3), dust optical properties and radiation flux diagnostics in CAM6.1 (Sect. 2.4), and our modifications to the base code (Sect. 2.5): the new dust emission scheme and change to the aerosol dry deposition and gravitational settling velocity to include dust asphericity. Two sets of simulations with offline dynamics (meteorology field nudged toward reanalysis data) were conducted (Table 2; Sect. 2.6) using bulk and speciated dust. A total of nine experiments were conducted to evaluate the performance of each development that a future version of official model release will likely include on reproducing the dust cycle against that of the current schemes and observations. Three out of the nine experiments quantify how the size treatment for transported dust affects the dust cycle modeling. We do not evaluate the model performance on simulating the dust cycle in the preindustrial considering the scarcity of measurements relative to the current climate (Mahowald et al., 2010).

Table 2 Simulations performed in this study for years 2006–2011. Treatment of dust tracers: speciated dust with separate tracers (mineralogy), or no dust speciation (bulk); the dust emission scheme: Zender et al. (2003a; DEAD) or Kok et al. (2014a; BRIFT); with or without accounting for the lifetime effect of dust asphericity (asp versus sph); dry deposition scheme: Zhang et al. (2001; Z01) or Petroff and Zhang (2010; PZ10); parameters for size distribution taken from the released version of CAM5 and CAM6.1 (see Table 1 for CAM5 and CAM6 sizes, respectively); additional tests on dust size distribution using the coarse-mode σ=1.2 from the released version of CAM6.1 and the other parameters (e.g., boundaries of the geometric mean diameter) from the released version of CAM5; dust tuning parameter includes the CAM namelist variable (dust_emis_fact) and “b” used in the calculation of the threshold gravimetric water content (see Sect. 2.5.1). The variable f_clay denotes the clay fraction in CLM5. CAM6.1 and CAM6.α in bold refer to the default model and proposed new model versions, respectively, with bulk dust. Note negligible influence on the dust cycle modeling and corresponding direct radiative effect by changing the size parameters of the accumulation mode between CAM5 and CAM6 size.

2.1 Aerosol representation

We used the Modal Aerosol Model version 4 (MAM4) in the CESM2.1-CAM6.1 (Liu et al., 2016). We consider both the default DEAD dust emission scheme (Zender et al., 2003a) in the current officially released version of CAM6.1 model as well as that of Kok et al. (2014a) (Sect. 2.5.1).

CAM6.1 simulates the advection, deposition, and aerosol microphysics (e.g., coagulation and nucleation) during transport via MAM4 using four log-normal size modes (Liu et al., 2016): Aitken (containing dust, sulfate, sea salt, and secondary organic matter), accumulation (containing sulfate, secondary organic matter, primary organic matter, black carbon, sea salt, and soil dust), coarse (containing dust, sea salt, and sulfate), and a primary carbonaceous mode (primary organic matter and black carbon). Within each mode, aerosol tracers are transported as an internal mixture of the species present, while aerosol species from different modes are externally mixed. Also advected in each of the four modes is the number concentration of aerosol particles (Liu et al., 2016), allowing an effective radius to be calculated and the effect of aerosol–cloud interactions to be diagnosed. The removal of dust aerosols is mainly through dry deposition and wet deposition (including in- and below-cloud processes, as detailed in Neale et al., 2010). In the formation of precipitating clouds, dust particles can serve as cloud condensation nuclei (CCN) and/or ice nucleating particles (INPs) and thus can be removed via nucleation scavenging (Zender et al., 2003a). In addition, the model accounts for the in-cloud scavenging of dust in the Aitken mode by Brownian diffusion, but neglects the other scavenging processes (Easter et al., 2004), which are relatively slow (Pruppacher et al., 1998), such as thermophoresis. Below the cloud, dust particles can be removed by the sub-cloud scavenging. This sub-cloud scavenging of dust aerosols follows a first order loss as the product of the precipitation flux, dust mass mixing ratio, and the scavenging efficiency (Dana and Hales, 1976), for example. The wet deposition rate thus depends on the hygroscopicity of dust (=  0.068; Scanza et al., 2015) as CCN/INPs and the prescribed scavenging coefficient (Neale et al., 2010), both of which are currently constant with respect to the dust size (and composition for speciated dust) in CAM6.1. This size independency of the scavenging coefficient may be an oversimplification, since measurements suggest that it can vary intensively on an order or two even within a size mode (Wang et al., 1978).

The geometric standard deviation of each mode is prescribed, and default values for CAM5 and CAM6.1 are given in Table 1, along with the initialized geometric mean diameter (GMD), based on which the model predicates the GMD online, and its ranges. Note that the current default CAM6.1 employs a narrow coarse-mode size distribution but a broad boundary width (high bound minus low bound), likely resulting in the GMD bounds less in effect, compared to that in CAM5. The narrower set of the coarse-mode size distribution was designed to accommodate for stratospheric aerosols (e.g., volcanic sulfate) (Mills et al., 2016), but was not previously compared to dust aerosol observations in detail.

2.2 Bulk-dust modeling

Parameterization of the default dust emissions in DEAD generally follows the dust mobilization mechanism developed by Marticorena and Bergametti (1995) (referred to as DEAD hereafter as well). As a component of CESM2, the CLM initiates dust entrainment once the near-surface friction velocity exceeds the soil threshold friction velocity, which primarily depends on the physical characteristics of the soil (e.g., soil moisture content and grain size distribution) and land cover (Kok et al., 2012; Shao, 2008). The downwind transfer of wind momentum to the surface soil to produce dust emissions is assumed to be completely prevented by vegetation when the leaf area index (LAI) exceeds a threshold value, 0.3 m² m⁻² (Mahowald et al., 2006a). Below the threshold value, the fraction of a grid cell capable of releasing dust aerosols is parameterized as an inverse and linear function of LAI (Mahowald et al., 2006a). The inhibition of soil moisture on dust deflation, and thus dust emission, activates when the near-surface soil gravimetric water content exceeds a threshold value, determined by the static mass fraction of the clay soil, and is parameterized in the land model according to a semi-empirical relation (Fécan et al., 1999) (see Sect. 2.5.1).

The size distribution of the emitted dust was derived using the brittle fragmentation theory developed by Kok (2011b) distributing 0.1 %, 1.0 %, and 98.9 % percentage of dust mass into Aitken, accumulation, and coarse modes, respectively, independent of the friction velocity upon dust emissions (Kok, 2011a).

2.3 Speciated-dust aerosol modeling

The bulk-dust model (Sect. 2.2) has previously been modified to speciate the bulk dust into eight mineral tracers, which allows more detailed optical properties as a function of minerals (Scanza et al., 2015; Li et al., 2021). Using the approach of Claquin et al. (1999), Li et al. (2021) estimated a mean mineralogical composition in the soil at each model grid cell for the minerals illite, kaolinite, montmorillonite, hematite, quartz, calcite, feldspar, and gypsum (Fig. S2 of Li et al., 2021). These minerals represent the most common classes for clay- (soil grain diameter < 2 µm including the first five minerals) and silt-sized (diameter between 2–63 µm including the last five minerals) soil categories (Claquin et al., 1999). As detailed in Scanza et al. (2015) and Li et al. (2021), additional modifications include (1) the mineral components in soil types of Gypsic Xerosols and Yermosols, Gleyic and Orthic Solonchaks and salt flats were normalized to unity; (2) the same amount of hematite in the clay- and silt-sized categories was prescribed with equal and opposite change to the illite percentage; (3) the nearest neighborhood algorithm was applied to fill in the grid cells for dust emission; and (4) the soil mineralogy was converted to that of the dust aerosol following the brittle fragmentation theory (Kok, 2011b), as detailed in Scanza et al. (2015).

The distribution of the mass flux for each mineral into the three emission modes follows that of the bulk-dust modeling (Sect. 2.2). The sum of the masses of the eight considered minerals equals the total bulk-dust mass without dust speciation. Each of the mineral aerosols are treated as a separate tracer in the same manner as bulk dust, experiencing advection, deposition, and aerosol microphysics (e.g., coagulation).

2.4 Dust optical properties and radiation flux diagnostic

We show results of the DRE calculations from two code versions: one with the bulk dust and the other with speciated dust. Aerosol optical properties (e.g., single scattering albedo and asymmetry factor) of the internal mixture in an aerosol mode are parameterized based upon the complex refractive index (CRI) of the mixture, which is calculated as the volume-weighted CRI of each component, including water (Ghan and Zaveri, 2007) in that mode. The wet size due to growth of aerosol particles by adsorbing water vapor follows the κ-Köhler theory with a time-invariant hygroscopicity for each aerosol species (Petters and Kreidenweis, 2007). CAM6.1 computes the net radiative flux based on the radiation fluxes diagnosed for each model layer at 14 shortwave and 16 longwave spectral bands per model hour. The DRE by dust aerosols under all-sky conditions was then determined by calculating the difference of the net radiative flux with and without dust at the top of the atmosphere. We augmented the longwave DRE from the model by 51 % to account for the dust scattering (Dufresne et al., 2002). The DRE efficiency, which we used to evaluate the model performance on simulating the dust optical properties, is defined as the ratio of dust DRE to dust optical depth (DOD) under clear conditions. This study does not consider the indirect radiative effect which is subject to substantially larger uncertainty due to the complexity involved in cloud microphysics (IPCC, 2021).

2.5 Changes to the dust parameterizations in CAM6.1/CLM5

The model developments introduced in this section are closely related to the three major components of the dust cycle (emission, transport, and removal mechanisms) and the radiative effects. Specifically, we incorporate into CAM6.1 the relatively new dust emission scheme originally developed by Kok et al. (2014a, b), the dry deposition scheme developed by Petroff and Zhang (2010) and incorporated in CAM5 by Wu et al. (2018), and the influence of dust non-spherical shape on the removal rate of dust aerosol particles (Huang et al., 2020).

2.5.1 Dust emission schemes

The vertical flux of dust emitted by wind erosion in a model grid cell is represented by

$$\begin{array}{}\text{(1)}& {\mathit{\phi }}_{\mathrm{d}}=\mathit{\lambda }{S}_{\mathrm{f}}{F}_{\mathrm{d},\mathrm{CLM}\mathrm{5}},\end{array}$$

where λ is a global tuning factor, S_f is the source function that shifts the dust emission to the most erodible sources, such as the Bodélé depression in North Africa (Zender et al., 2003a), and F_d,CLM5 is the vertical emission flux predicted by the dust emission scheme in CLM5 (Kok et al., 2014a).

As part of the DEAD scheme (Zender et al., 2003a), dust sources are strongly associated with the erodible soils (Ginoux et al., 2001). These source regions are parameterized using information contained in the time invariant geomorphology map (Zender et al., 2003b) which was optimized (Albani et al., 2014) to match the observed DOD. The vertical dust emission in CLM5 occurs when the friction velocity ($u_{\ast }^{\prime }$) exceeds the threshold friction velocity ($u_{\ast \mathrm{t}}^{\prime }$) on the entire surfaces including nonerodible roughness elements, which is parameterized in DEAD as

$$\begin{array}{}\text{(2)}& \begin{array}{rl}& F_{\mathrm{d},\mathrm{CLM}\mathrm{5}}=\\ & \left\{\begin{array}{l}{C}_{\mathrm{MB}}\mathit{\chi }{f}_{\mathrm{bare}}\frac{{\mathit{\rho }}_{\mathrm{a}}}{g}{u}_{\ast }^{\prime \mathrm{3}}\left(\mathrm{1}-\frac{u_{\ast \mathrm{t}}^{\prime \mathrm{2}}}{u_{\ast }^{\prime \mathrm{2}}}\right)\left(\mathrm{1}+\frac{u_{\ast \mathrm{t}}^{\prime }}{u_{\ast }^{\prime }}\right),u_{\ast }^{\prime }>u_{\ast \mathrm{t}}^{\prime }\\ \mathrm{0},u_{\ast }^{\prime }\le u_{\ast \mathrm{t}}^{\prime }\end{array}\right.,\end{array}\end{array}$$

where C_MB is a dimensionless proportionality constant; f_bare is the bare soil fraction; ρ_a is the atmospheric density; and the sandblasting efficiency χ is written as a function of the clay fraction (f_clay),

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

Kok et al. (2014a) developed a new dust emission scheme for climate models based on the brittle fragmentation theory (Kok, 2011b), which avoided the use of a static soil erodibility map (the source function, S_f in Eq. 1) while improving the accuracy of the dust cycle modeling (Kok et al., 2014b); although even dust modeling with DEAD can be improved if optimized against observations (Kok et al., 2021b). Improvements are likely achieved because, compared to that in DEAD, the dust emission in BRIFT tends to be more sensitive to the soil's threshold friction velocity and thus to the surface physical conditions when soil becomes more erodible, owing to the introduced dust emission coefficient (Kok et al., 2014a) and the new method of calculating the threshold gravimetric water content in the top soil layer (see Eq. 4 of Kok et al., 2014b). Below we briefly introduce the new dust emission scheme.

In this new scheme the vertical dust emission in CLM5 is expressed as

$$\begin{array}{}\text{(4)}& \begin{array}{rl}& F_{\mathrm{d},\mathrm{CLM}\mathrm{5}}=\\ & \left\{\begin{array}{l}{C}_{\mathrm{d}}{f}_{\mathrm{bare}}{f}_{\mathrm{clay}}\frac{{\mathit{\rho }}_{\mathrm{a}}\left(u_{\ast }^{\mathrm{2}}-u_{\ast \mathrm{t}}^{\mathrm{2}}\right)}{u_{st}}{\left(\frac{u_{\ast }}{u_{\ast \mathrm{t}}}\right)}^{C_{\mathit{\alpha }}\frac{u_{\ast \mathrm{st}}-u_{\ast \mathrm{st}\mathrm{0}}}{u_{\ast \mathrm{st}\mathrm{0}}}},u_{\ast }>u_{\ast \mathrm{t}}\\ \mathrm{0},u_{\ast }\le u_{\ast \mathrm{t}}\end{array}\right.\end{array},\end{array}$$

where u_∗ and u_∗t are friction velocity and threshold friction velocity on the bare erodible soil surfaces; u_∗st is the threshold friction velocity standardized according to atmospheric density and the standard atmospheric density of ρ_a0=1.225 kg m⁻³,

$$\begin{array}{}\text{(5)}& u_{\ast \mathrm{st}}=u_{\ast \mathrm{t}}\sqrt{{\displaystyle \frac{{\mathit{\rho }}_{\mathrm{a}}}{{\mathit{\rho }}_{\mathrm{a}\mathrm{0}}}}};\end{array}$$

and u_∗st0 is the minimal value of u_∗st equaling 0.16 m s⁻¹.

Since the presence of non-erodible roughness elements, such as rocks, is not considered in CLM5,

$$\begin{array}{}\text{(6)}& {\displaystyle }{u}_{\ast }=u_{\ast }^{\prime },\text{(7)}& {\displaystyle }{u}_{\ast \mathrm{t}}=u_{\ast \mathrm{t}}^{\prime }.\end{array}$$

Because of the neglection of the non-erodible elements, u_∗t is mostly determined by soil moisture content, which means that the augmentation factor of u_∗t is

$$\begin{array}{}\text{(8)}& f_{\ast \mathrm{t}}=\left\{\begin{array}{l}\sqrt{\mathrm{1}+\mathrm{1.21}{\left(w-w^{\prime }\right)}^{\mathrm{0.68}}};w>w^{\prime }\\ \mathrm{1};w\le w^{\prime }\end{array}\right.,\end{array}$$

where w and w^′ are soil moisture content and the threshold gravimetric water content of the top soil layer in percentage.

Fécan et al. (1999) parameterized the threshold gravimetric water content (w) of the top soil layer by

$$\begin{array}{}\text{(9)}& w^{\prime }=b\left(\mathrm{17}{f}_{\mathrm{clay}}+\mathrm{14}{f}_{\mathrm{clay}}^{\mathrm{2}}\right),\end{array}$$

where b is a tuning factor.

Equations (8) and (9) are also used in DEAD with an equivalent tuning factor b set to be $f_{\mathrm{clay}}^{-\mathrm{1}}$ which in BRIFT is set as unity. The clay fraction is taken from the FAO et al. (2012) soil database (see Fig. S1 of Kok et al., 2014b).

The dust emission coefficient, C_d, in Eq. (4) is expressed as

$$\begin{array}{}\text{(10)}& C_{\mathrm{d}}=C_{\mathrm{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_d0 equals $\mathrm{4.4}\times {\mathrm{10}}^{-\mathrm{5}}$ and C_e equals 2.0.

2.5.2 Dry deposition schemes

The default dry deposition scheme, Z01

As is typical among aerosols dry deposition resistance models, CAM6.1 includes parameterizations of gravitational settling (V_g), and aerodynamic (R_a) and surface resistance (R_s).

The gravitational settling is parameterized following

$$\begin{array}{}\text{(11)}& V_{\mathrm{g}}={\displaystyle \frac{{\mathit{\rho }}_{\mathrm{p}}g{C}_{\mathrm{c}}{d}_{\mathrm{p}}^{\mathrm{2}}}{\mathrm{18}\mathit{\mu }}},\end{array}$$

where ρ_p is the particle density (unit: kg m⁻³), d_p is the particle diameter (unit: m), g the acceleration of gravity, C_c is the Cunningham correction factor as a function of d_p and the mean free path of air molecules, and μ is the viscosity coefficient of air.

Aerodynamic resistance is parameterized following

$$\begin{array}{}\text{(12)}& R_{\mathrm{a}}={\displaystyle \frac{\ln \left(Z_{\mathrm{R}}/Z_{\mathrm{0}}\right)-{\mathit{\phi }}_H}{\mathit{\kappa }{u}_{\ast }}},\end{array}$$

where Z_R is the reference height; Z₀ is the roughness length; φ_H is the stability function, κ is the von Karman constant set as 0.4; and u_∗ is the friction velocity.

The surface resistance dominates over aerodynamic resistance under turbulent conditions, and is written as

$$\begin{array}{}\text{(13)}& R_{\mathrm{s}}={\displaystyle \frac{\mathrm{1}}{{\mathit{\epsilon }}_{\mathrm{0}}{u}_{\ast }\left(E_{\mathrm{B}}+E_{\mathrm{IM}}+E_{\mathrm{IN}}\right)R_{\mathrm{1}}}},\end{array}$$

where ε₀ is an empirical constant set as 3.0; R₁ is the factor to represent particle rebound; E_B,E_IM, and E_IN are collection efficiencies due to Brownian diffusion, impaction, and interception, parameterized respectively as

$$\begin{array}{}\text{(14)}& {\displaystyle }{E}_{\mathrm{B}}=S{c}^{-\mathit{\gamma }},\text{(15)}& {\displaystyle }{E}_{\mathrm{IM}}={\left({\displaystyle \frac{St}{St+\mathit{\alpha }}}\right)}^{\mathit{\beta }},\text{(16)}& {\displaystyle }{E}_{\mathrm{IN}}=\mathrm{0.5}{\left({\displaystyle \frac{d_{\mathrm{p}}}{A}}\right)}^{\mathrm{2}},\end{array}$$

where Sc is the Schmidt number, defined as the ratio of the kinematic viscosity of air ($\mathrm{1.57}\times {\mathrm{10}}^{-\mathrm{5}}$ m² s⁻¹) to the particle Brownian diffusivity (D); γ depends on land use categories, typically ranging between [0.50, 0.67]; α depends on the land use categories; β=2; A is the characteristic radius of collectors depending on land use categories; and St is the Stokes number parameterized following

$$\begin{array}{}\text{(17)}& St=\left\{\begin{array}{ll}\frac{V_{\mathrm{g}}{u}_{\ast }}{gA};& \text{Vegetated surfaces}\\ \frac{V_{\mathrm{g}}{u}_{\ast }^{\mathrm{2}}}{v}.& \text{Smooth surfaces or surfaces with}\\ & \text{bluff roughness elements}\end{array}\right.\end{array}$$

According to Eq. (13), the surface resistance consists of three processes, two applicable to all land types (Brownian diffusion and impaction), and one only to non-smooth surfaces (interception). All the three processes are a function of aerosol size through empirical coefficients constrained by matching the modeled dry deposition velocity with field and laboratory measurements.

The dry deposition velocity then has the form of

$$\begin{array}{}\text{(18)}& V_{\mathrm{d}}=V_{\mathrm{g}}+{\displaystyle \frac{\mathrm{1}}{R_{\mathrm{a}}+R_{\mathrm{s}}}}.\end{array}$$

With more observations available to constrain these coefficients, the default Z01 (Zhang et al., 2001) used in CAM6.1 was found to greatly overestimate dry deposition rates for fine particles (diameter < 1 µm: Aitken and accumulation modes) and slightly underestimate (relative to the large change with fine particles) the rates for coarse particles (diameter around 1 or 2 µm) (Petroff and Zhang, 2010; Wu et al., 2018).

The new dry deposition scheme, PZ10

The new scheme (PZ10; Petroff and Zhang, 2010) uses a quite different formula to calculate the dry deposition scheme as follows:

$$\begin{array}{}\text{(19)}& V_{\mathrm{d}}=V_{\mathrm{drift}}+{\displaystyle \frac{\mathrm{1}}{R_{\mathrm{a}}+\mathrm{1}/R_{\mathrm{ds}}}},\end{array}$$

where the drift velocity,

$$\begin{array}{}\text{(20)}& V_{\mathrm{drift}}=V_{\mathrm{g}}+V_{\mathrm{phor}}.\end{array}$$

Therefore, this new scheme includes the effect (V_phor) of different physical processes (thermophoresis, diffusiophoresis, and electricity) occurring between water, ice, and snow surfaces and the air immediately above them, which can result in a downward flux of particles (the phoretic effect; Petroff et al., 2008). PZ10 accounts for such effects of thermophoresis and diffusiophoresis for particle deposition over the three surface types by assigning constant values of $\mathrm{5}\times {\mathrm{10}}^{-\mathrm{5}}$ m s⁻¹ to water and $\mathrm{2}\times {\mathrm{10}}^{-\mathrm{4}}$ m s⁻¹ to ice and snow surfaces, which allows the scheme to better reproduce the available measurements than Z01 (Petroff and Zhang, 2010). This constant is set to zero for all the other surface types. The phoretic effect tends to dominate deposition of fine particles over Brownian diffusion under low wind conditions (friction velocity less than ∼ 11 cm s⁻¹). Because of the reduced Brownian diffusion efficiency compared to Z01, PZ10 corrects the high bias seen in Z01 for the deposition of fine particles (Emerson et al., 2020; Petroff and Zhang, 2010; Wu et al., 2018).

PZ10 parameterizes the gravitational settling velocity V_g in the same way as Z01 but the aerodynamic resistance (R_a) in a different formula.

For non-vegetated surfaces,

$$\begin{array}{}\text{(21)}& R_{\mathrm{a}}={\displaystyle \frac{\mathrm{1}}{\mathit{\kappa }{u}_{\ast }}}\left[\ln \left({\displaystyle \frac{Z_{\mathrm{R}}-d}{Z_{\mathrm{0}}}}\right)-{\mathrm{\Psi }}_{\mathrm{h}}\left({\displaystyle \frac{Z_{\mathrm{R}}-d}{L_{\mathrm{o}}}}\right)+{\mathrm{\Psi }}_{\mathrm{h}}\left({\displaystyle \frac{Z_{\mathrm{0}}}{L_{\mathrm{o}}}}\right)\right],\end{array}$$

where d and L_o are the displacement height of the canopy and the Obhukov length, respectively; and Ψ_h is the integrated form of the stability function for heat.

The surface dry deposition velocity is expressed as

$$\begin{array}{}\text{(22)}& V_{\mathrm{ds}}=u_{\ast }\left(E_{\mathrm{gb}}+E_{\mathrm{IT}}\right),\end{array}$$

where E_IT is the efficiency of collections by turbulent impaction; E_gb represents Brownian diffusion, written as

$$\begin{array}{}\text{(23)}& \begin{array}{rl}{E}_{\mathrm{gb}}& ={\displaystyle \frac{S{c}^{-\mathrm{2}/\mathrm{3}}}{\mathrm{14.5}}}\left[{\displaystyle \frac{\mathrm{1}}{\mathrm{6}}}\ln {\displaystyle \frac{{\left(\mathrm{1}+\sqrt[\mathrm{3}]{Sc}/\mathrm{2.9}\right)}^{\mathrm{2}}}{\mathrm{1}-\sqrt[\mathrm{3}]{Sc}/\mathrm{2.9}+{\left(\sqrt[\mathrm{3}]{Sc}/\mathrm{2.9}\right)}^{\mathrm{2}}}}\right.\\ & {\left.+{\displaystyle \frac{\mathrm{1}}{\sqrt{\mathrm{3}}}}{\tan }^{-\mathrm{1}}{\displaystyle \frac{\mathrm{2}\left(\sqrt[\mathrm{3}]{Sc}/\mathrm{2.9}\right)-\mathrm{1}}{\sqrt{\mathrm{3}}}}+{\displaystyle \frac{\mathit{\pi }}{\mathrm{6}\sqrt{\mathrm{3}}}}\right]}^{-\mathrm{1}}.\end{array}\end{array}$$

For vegetated surfaces,

$$\begin{array}{}\text{(24)}& \begin{array}{rl}{R}_{\mathrm{a}}& ={\displaystyle \frac{\mathrm{1}}{\mathit{\kappa }{u}_{\ast }}}\left[\ln \left({\displaystyle \frac{Z_{\mathrm{R}}-d}{h-d}}\right)-{\mathrm{\Psi }}_{\mathrm{h}}\left({\displaystyle \frac{Z_{\mathrm{R}}-d}{L_{\mathrm{o}}}}\right)\right.\\ & \left.+{\mathrm{\Psi }}_{\mathrm{h}}\left({\displaystyle \frac{h-d}{L_{\mathrm{o}}}}\right)\right],\end{array}\end{array}$$

where h is the canopy height.

The expression for surface deposition velocity is written as

$$\begin{array}{}\text{(25)}& V_{\mathrm{ds}}=u_{\ast }{E}_{\mathrm{g}}{\displaystyle \frac{\mathrm{1}+\left(\frac{Q}{Q_{\mathrm{g}}}-\frac{\mathit{\delta }}{\mathrm{2}}\right)\frac{\tanh \mathit{\eta }}{\mathit{\eta }}}{\mathrm{1}+\left(Q_{\mathrm{g}}+\frac{\mathit{\delta }}{\mathrm{2}}\right)\frac{\tanh \mathit{\eta }}{\mathit{\eta }}}},\end{array}$$

where E_g is the collection efficiency on the ground below the canopy.

E_g includes Brownian diffusion and the turbulent impaction, E_gt, given as

$$\begin{array}{}\text{(26)}& E_{\mathrm{gt}}=\mathrm{2.5}\times {\mathrm{10}}^{-\mathrm{3}}{C}_{\mathrm{IT}}{\mathit{\tau }}_{\mathrm{ph}}^{\mathrm{2}},\end{array}$$

where C_IT equals 0.14 and τ_ph is the dimensionless particle relaxation time.

Q and Q_g in Eq. (25) are given as

$$\begin{array}{}\text{(27)}& {\displaystyle }Q=\mathrm{LAI}\cdot {\displaystyle \frac{\left(E_{\mathrm{B}}+E_{\mathrm{IN}}+E_{\mathrm{IM}}\right)U_h/u_{\ast }+E_{\mathrm{IT}}}{I_{\mathrm{mp}}\left(h\right)/h}},\text{(28)}& {\displaystyle }{Q}_{\mathrm{g}}={\displaystyle \frac{E_{\mathrm{g}}h}{I_{\mathrm{mp}}\left(h\right)}},\end{array}$$

where LAI is the two-sized leaf area index; E_B, E_IN, and E_IM, are Brownian diffusion, interception, and inertial impaction, respectively.

In Eq. (25),

$$\begin{array}{}\text{(29)}& \mathit{\eta }=\sqrt{{\mathit{\delta }}^{\mathrm{2}}/\mathrm{4}+Q},\end{array}$$

where the aerodynamic extinction coefficient δ is expressed as

$$\begin{array}{}\text{(30)}& \mathit{\delta }={\left({\displaystyle \frac{\mathrm{LAI}\cdot k_x}{\mathrm{12}{\mathit{\kappa }}^{\mathrm{2}}{\left(\mathrm{1}-d/h\right)}^{\mathrm{2}}}}\right)}^{\mathrm{1}/\mathrm{3}}{\mathit{\varphi }}_{\mathrm{m}}^{\mathrm{2}/\mathrm{3}}\left({\displaystyle \frac{h-d}{L_{\mathrm{o}}}}\right),\end{array}$$

where k_x is the inclination coefficient of the canopy elements and ϕ_m is the nondimensional stability function for momentum.

2.5.3 Dust asphericity

To account for the influence of dust asphericity on the gravitational settling velocity, we first calculated the asphericity factor γ (defined as the ratio of the gravitational settling velocity of aspherical dust to that of spherical dust) offline based on a combination of observed dust shape parameters (F_s) previously compiled by Huang et al. (2020) following

$$\begin{array}{}\text{(31)}& \mathit{\gamma }={\displaystyle \frac{\mathrm{2}}{F_{\mathrm{s}}^{\mathrm{1}/\mathrm{3}}+\mathrm{1}/F_{\mathrm{s}}^{\mathrm{1}/\mathrm{3}}}},\end{array}$$

F_s is parameterized by the dust shape parameters,

$$\begin{array}{}\text{(32)}& F_{\mathrm{s}}={\displaystyle \frac{D_{\mathrm{g}}^{\mathrm{3}}}{L^{\mathrm{2.3}}\cdot W^{\mathrm{0.7}}}},\end{array}$$

where D_g is the volume-equivalent diameter of the dust particle defined by three axes (L: length, W: width, and H: height, respectively) of the ellipsoid having a form of

$$\begin{array}{}\text{(33)}& D_{\mathrm{g}}=\sqrt[\mathrm{3}]{\mathrm{LWH}}.\end{array}$$

The orientation of particles during gravitational settling determines the drag coefficient. Equation (31) assumes that during settling, aerosol particles randomly orientate. This assumption is reasonable, since for dust falling in the Earth's atmosphere, (1) the Reynolds number, Re≪1 (Kok et al., 2012), and especially (2) CAM6 does not simulate super coarse dust particles (diameter > 10 µm), for which such an assumption may introduce high errors. A previous study (Bagheri and Bonadonna, 2016) suggests that this approximation of the influence of the dust asphericity on the gravitational settling velocity is accurate and reliable with the mean and the maximum errors of 2.4 %, and 33.9 %, respectively. Equations (31)–(33) indicate a range of γ between 0 and 1. When a dust particle becomes less ellipsoidal, γ is getting closer to 1.

In the Stokes regime (Kok et al., 2012), where the gravitational settling of dust usually occurs, the terminal velocity of spherical (sph) and ellipsoidal (asp) dust is approximated as

$$\begin{array}{}\text{(34)}& V_{\mathrm{g},\mathrm{sph}}={\displaystyle \frac{g{\mathit{\rho }}_{\mathrm{p}}}{\mathrm{18}\mathit{\mu }}}{D}_{\mathrm{g}}^{\mathrm{2}},\end{array}$$

and

$$\begin{array}{}\text{(35)}& V_{\mathrm{g},\mathrm{asp}}=\mathit{\gamma }{V}_{\mathrm{g},\mathrm{sph}},\end{array}$$

respectively, where g is the gravitational constant (∼ 9.8 m s⁻²), ρ_p is the dust density (∼ 2500 kg m⁻³), and μ is the dynamic viscosity of air ($\mathrm{1.81}\times {\mathrm{10}}^{-\mathrm{5}}$ Pa s).

In this calculation, we also assume that the dust shape parameters are independent of the size of dust aerosol particles. Therefore, a constant revision of the dust gravitational settling velocity (the calculated value in the model by default is for spherical aerosols) due to dust asphericity by multiplying the velocity by γ was applied to dust species in the three modes that contains dust aerosols (Aitken, accumulation, and coarse). The size independency assumption of dust asphericity follows the recent observational evidence that there does not exist a statistically significant relationship between the shape parameters (aspect ratio and height-to-width ratio) and dust sizes (Huang et al., 2020).

Measurements made at different locations show that the shape parameters (e.g., aspect ratio; Fig. 3 of Huang et al., 2020), which we used to calculate V_g,asp, change for dust during transport. But, because of highly limited measurements of dust shape parameters, we subjectively divided the dust coverage into “close-to-source”, “short-range”, and “long-range” zones and calculated the asphericity factor γ for each of the zones, the global map of which is shown in Fig. S1 in the Supplement, ranging between 0.82 and 0.93. In the regions where the shape parameter measurement is sparse or unavailable, such as those in the Southern Hemisphere, the shape parameters from the global median were used instead to calculate the asphericity factor yielding a value of 15 %.

We acknowledge limitation of the methodology here to account for the lifetime effect of dust asphericity, anticipating improvements on modeling this effect when more high-quality dust shape measurements become available.

Matching modeled DOD to observations requires the model to account for dust asphericity, which acts to enhance the mass extinction efficiency of particles, particularly in the coarse mode (Kok et al., 2017). This enhancement in the mass extinction efficiency due to the dust asphericity is not included in the current version of CESM2.1 but will be likely incorporated into a future officially released CESM version. According to calculations of Kok et al. (2017), the dust mass extinction efficiency at the visible band due to dust asphericity is approximately 16 % and 28 % higher for non-spherical particles than for spherical particles in the fine (accumulation plus Aitken) and coarse modes, respectively. Consequently, the model requires lower dust emissions to achieve a global DOD of ∼ 0.030 compared to simulations without considering dust asphericity. The shape effect on the mass extinction efficiency may also explain the difference between the global mean DOD in Aerosol Comparisons between Observations and Models (AEROCOM; median: 0.023) (Huneeus et al., 2011) and that in Ridley et al. (2016) (0.03 ± 0.005) near the visible band. We have included the enhanced dust mass extinction efficiency due to dust asphericity in our previous studies (e.g., Li et al., 2021; Kok et al., 2021b), which suggests that the inclusion of this enhanced dust mass extinction efficiency would reduce the overestimation of the surface concentration (Kok et al., 2021b). Here we do not investigate its impact on the simulated dust cycle.

2.6 Experiment design

Table 2 lists the simulations designed for the present study. In all simulations, the CAM6.1 with different modifications was configured as a standalone model where the atmosphere is coupled to active land and sea ice models, and to data ocean and slab glacier models. Each simulation in these sets was performed at the spatial resolution of 1.25^∘ × 0.9^∘ × 56 (longitude by latitude by vertical layers) using a data ocean for years 2006–2011, with the simulated data for the last 5 years used for analysis. In addition, the meteorology field (horizontal wind, air temperature T, and relative humidity) was nudged toward the Modern-Era Retrospective analysis for Research and Applications, version 2 (MERRA-2) at a 6 h relaxation timescale. The anthropogenic emissions were taken from the Climate Model Intercomparison Program (CMIP6) inventory for the year 2000 (Eyring et al., 2016).

The enhancement of the mass extinction efficiency of aerosol particles by dust asphericity is included in all the simulations, since we do not attempt to quantify how this enhancement impacts the simulated dust cycle. An offline sensitivity test (Table S1 in the Supplement) supports the use of unity tuning factor to calculate the threshold gravimetric water content which we employed in the experiments for quantifying influence of each modification (speciated-dust simulations listed in Table 2).

NEW_EMIS serves as the baseline simulation for quantifying the impact of the coarse-mode size change on the dust cycle modeling. NEW_EMIS_SIZE completely and NEW_EMIS_SIZE_WIDTH partly reverted the coarse-mode size distribution to that used in CAM5 (Table 1). The other changes to the width of the accumulation mode and the bounds of the simulated GMD online impose negligible impacts on the dust cycle modeling, thus we did not construct sensitivity tests on them in this study. We investigate how the incorrect dust size distribution influences the dust cycle modeling and the estimate of dust DRE in the bulk-dust model rather than in the speciated-dust model, because this incorrect size distribution has been employed in previous studies using the officially released bulk-dust CAM6 only and not in any study using the speciated-dust CAM. It is also reasonable to make all the quantifications in the model that uses a correct dust size distribution. Therefore, we reverted the dust size distribution in all the speciated-dust runs to that configured in CAM5.

We quantified the impact of each of the modifications (Z01 to PZ10, spherical to aspherical dust, and DEAD to BRIFT) on the simulated dust cycle and DRE by differentiating corresponding results in the paired simulations that contain identical developments except for the targeted modification. Specifically, we quantified the impact of changing the following: (1) Z01 to PZ10 by taking the difference between the simulation with Z01 (MINE_NEW_ EMIS_SHAPE) and that with PZ10 (CAM6.α_MINE), (2) spherical to aspherical dust between the simulation with spherical dust (MINE_NEW_EMIS) and that with aspherical dust (MINE_NEW_EMIS_SHAPE), and (3) DEAD to BRIFT between the simulation using DEAD (MINE_BASE) and that using BRIFT (MINE_NEW_EMIS).

Note that there are many ways to conduct sensitivity studies, which could lead to slightly different results. We added the modification on top of the previous change to understand how the simulated dust cycle evolves while updating the model (MINE_BASE) toward the most advanced version (CAM6.α_MINE). This may not hinder a clean comparison of the effect of each development, since the “interaction” between the existing and the newly introduced parameterizations appears weak (Fig. S2).

With the dust tuning applied toward the similar global mean DOD of ∼ 0.030, the modeled dust cycle (i.e., burdens, concentrations, loadings, and deposition fluxes) would be similar between the bulk- and speciated-dust models that are nudged toward identical offline dynamics and are using the same dust size distribution (see Sect. 6). The quantified effect of each of the modifications would thus be similar if using the bulk-dust model instead (Fig. S2), except that the modeled dust optical properties (e.g., single scattering albedo) by the bulk- and speciated-dust models would differ considerably, resulting in considerably different dust DRE (Scanza et al., 2015) and DRE efficiencies between NEW_EMIS (CAM6.α) and MINE_NEW_EMIS (CAM6.α_MINE). A comparison of the bulk- and speciated-dust models on simulating dust DRE had been previously documented (Scanza et al., 2015). This study includes the speciated-dust runs, because we want to verify as well if the updates help improve the agreement with the observed dust DRE efficiency in the speciated-dust model which could better represent the spatial variation of the dust optical properties.

We tuned CAM6.1, NEW_EMIS, CAM6.α, MINE_BASE, and MINE_NEW_EMIS, following Albani et al. (2014), by modifying a CAM namelist variable, dust_emis_fact, such that the simulated global mean DOD is ∼ 0.030 at the visible band centered at 0.53 µm (hereafter unless stated otherwise), an estimate obtained by an integrated analysis of the Aerosol Robotic Network (AERONET)-based measurements, bias-corrected satellite retrievals, and a model ensemble (Ridley et al., 2016). We prefer to tuning the model to reproduce the global mean DOD, 0.030, because DOD is currently the best estimate of global dust quantities compared to the others (i.e., dust concentrations). It turns out that doing so can also reasonably reproduce the other quantities with no need of a regional tuning. MINE_NEW_EMIS requires the dust tuning to use a much larger tuning parameter (dust_emis_fact = 3.6; Table 2) than MINE_BASE (dust_emis_fact = 1.6), because otherwise, if using the same dust_emis_fact as in DEAD, the dust emissions in BRIFT would lead to an unrealistically high global mean DOD (> ∼ 0.50).

The dust tuning was not applied to NEW_EMIS_SIZE and NEW_EMIS_SIZE_WIDTH, but the emissions in which were kept identical to NEW_EMIS, to see how changes in the transported dust size distribution affects the DOD calculation. Because of the rough linearity among DOD, DRE, and dust burdens (Liao and Seinfeld, 1998; Mahowald et al., 2006b), when comparing surface dust concentrations, dust loadings, and deposition fluxes, we rescaled each of them using the same factor to achieve the global mean DOD ∼ 0.030. For the other cases (MINE_NEW_EMIS_SHAPE and CAM6.α_MINE), as will be seen, the global mean DOD only changes slightly within the uncertainty range (0.025–0.035; Ridley et al., 2016). The model retuning is thus not required.

3 Observational datasets for model evaluations

Tables 3–5 summarize available datasets used to evaluate the model performance, detailed descriptions about each datum, and how they are used in the model–data comparison. Due to limitations in precisely matching the period and locations between model results and data, the evaluations focus on checking if models can capture overall features of the observed dust cycle and the corresponding dust DRE efficiency. We summarize limitations going beyond this mismatch on period and location and common in all the model–data comparisons in Sect. 7.

3.1 Surface dust concentrations and dust aerosol optical depth from AERONET

We used monthly surface dust concentration data that Albani et al. (2014) compiled from measurements made using high-volume filter collectors at the University of Miami Ocean Aerosol Network and station-based data that have been previously compiled on annual averages (Mahowald et al., 2009; Zuidema et al., 2019). Because the model only simulates dust < 10 µm (the cut-off value of aerosol size) in diameter, the Albani et al. (2014) compilation had been processed to estimate the flux of dust below the size cut-off according to reported or assumed parameters (e.g., geometric standard deviation) for the size distribution of transported dust. Simulated DOD is compared to AERONET retrievals, which were subject to data quality control and station selection based on the dust dominance in the reported aerosol optical depth (AOD). The data quality control includes a minimum of 10 d per month that contains valid retrievals, the annually averaged Angstrom Exponent < 1.2 (the larger the value, the smaller the aerosol size), and a full coverage of the data availability through a year within the observation period (Albani et al., 2014).

3.2 Surface dust deposition fluxes

The dust deposition flux data used here include those that Albani et al. (2014) compiled from publications (Tegen et al., 2002; Ginoux et al., 2001; Lawrence and Neff, 2009; Mahowald et al., 2009) for the present-day climate. Since the model only simulates dust < 10 µm in diameter, the Albani et al. (2014) compilation was processed to estimate the surface deposition fluxes of dust below the size cut-off, according to reported or assumed parameters (e.g., geometric standard deviation) for the size distribution of transported dust.

Table 3 Observed/retrieved cycle for dust model evaluations including optical depth, surface mass concentrations, surface deposition fluxes, and wet deposition percentages. AERONET: Aerosol Robotic Network; MODIS: Moderate Resolution Imaging Spectroradiometer; AOD: aerosol optical depth; DOD: dust optical depth.

3.3 Size distributions of dust aerosol

Most of the remotely sensed, size-resolved dust volume retrievals used here were taken from the AERONET level 2.0 almucantar retrievals (version 2), which is reported for 22 size bins with bimodal size distribution and ellipsoid shape of aerosol particles (Dubovik et al., 2000). These data overestimate dust mass in the submicron size range and have possible contamination by non-dust aerosols (Albani et al., 2014; Dubovik et al., 2000; Mahowald et al., 2014). We therefore only retained the supermicron fraction of dust in the comparison, even though AERONET may underestimate the mass of dust between 1–10 µm in diameters (McConnell et al., 2008). The data processing procedure is detailed in Albani et al. (2014). Near North Africa, we also compare the modeled size distribution of dust aerosols with measurements from Otto et al. (2007) taken in the vicinity of the Canary Islands, from Ryder et al. (2013) by aircraft with a track between the Canary Islands and Mauritania/Mali, and from Ryder et al. (2018) near Cape Verde.

Table 4 Measured/retrieved dust size distribution for model evaluations. AERONET: Aerosol Robotic Network; DustCOMM: Dust Constraints from joint Observational–Modelling–experiMental analysis.

3.4 The direct radiative effect efficiency of dust

The modeled dust DRE efficiency (ratio of dust DRE to DOD) is compared to satellite-based observations under clear-sky conditions at the top of the atmosphere. These observations include the following:

1. Longwave dust DRE efficiency derived over North Africa by Zhang and Christopher (2003) for September 2000 based on measured longwave (5–200 µm) fluxes at the top of the atmosphere from the Multi-angle Imaging SpectroRadiometer (MISR), and the Clouds and the Earth's Radiant Energy System (CERES) instrument and AOD at 0.55 µm from Moderate Resolution Imaging Spectroradiometer (MODIS).

2. The shortwave (0.3–5 µm) dust DRE efficiency obtained by Li et al. (2004) with the measured shortwave flux from CERES and AOD at 0.55 µm from MODIS near North Africa (15–25^∘ N, 45–15^∘ W) in the summer to winter months between 2000 and 2001.

3. The shortwave (0.3–5 µm) dust DRE efficiency that Patadia et al. (2009) derived using a 1-D radiative transfer model, radiative fluxes from CERES, and AOD at 0.55 µm from MISR and Ozone Monitoring Instrument (OMI) over the high-reflective regions (surface albedo and 0.55 µm > 0.35) of Saharan desert (15–30^∘ N, 10^∘ W–30^∘ E) for the summer months in 2005 and 2006.

Table 5 Retrieved dust direct radiative effect efficiency (DREE) for model evaluations. CERES: Clouds and the Earth's Radiant Energy System; TOA: top of the atmosphere; JJA: June, July, and August; AOD: aerosol optical depth; MISR: Multi-angle Imaging SpectroRadiometer; OMI: Ozone Monitoring Instrument; NDJ: November, December, and January; MODIS: Moderate Resolution Imaging Spectroradiometer; CALIPSO: Cloud–Aerosol Lidar and Infrared Pathfinder Satellite Observations; MFRSR: MultiFilter Rotating Shadowband Radiometer; SEVIRI: Spinning Enhanced Visible and Infrared Imager; GERB: Geostationary Earth Radiation Budget; AERONET: Aerosol Robotic Network; MPL: micro-pulse lidar; AERI: Atmospheric Emitted Radiance Interferometer; SMART: Surface-sensing Measurements for Atmospheric Radiative Transfer; AMJ: April, May, and June.

3.5 Other datasets

In addition to the abovementioned observations, we compare our results to datasets which combine model simulations and observations. Specifically, we compare (1) the modeled transported dust size distribution with that from the Dust Constraints from joint Observational–Modelling–experiMental analysis (DustCOMM) (Adebiyi et al., 2020) in global average; (2) regional dust deposition fluxes with the semi-observational data that were inverted based on an integration of a global model ensemble and quality-controlled observational constraints on the transported dust size distribution, extinction efficiency, and regional DOD (Kok et al., 2021a); this set of semi-observational data were shown to compare better with the high-quality measurement than model ensemble means or any individual model (Kok et al., 2021a); and (3) regional DOD in different seasons with the estimates of Ridley et al. (2016), who obtained DOD by combining four global climate models with multiple satellite aerosol products that were bias corrected using station-based AEROENT data.

4 Model assessment metrics

Metrics used to evaluate the model performance against observations include the root mean square error (RMSE) and correlation efficient (Kendall's τ or Spearman's correlation). Both the Kendall's τ and Spearman's correlation are non-parametric methods which do not require a distribution of the data, such as Gaussian or normal. For dust deposition and loadings, correlations were calculated to assess how well models reproduce both their regional climatology mean or one-time observation and the seasonal cycles. However, because of a lack of reliable monthly data, assessments for the dust DRE efficiency, DOD from Ridley et al. (2016) and percentages of wet deposition in the total deposition on spatial variability are based on the regional climatology mean or one-time observations. We tested the correlation significance of the metrics at the statistical confidence level of 95 %. For the dust DRE efficiency and percentages of wet deposition, some domains only have a range available, such as the Sahara Desert (15–30^∘ N, 10^∘ W–30^∘ E) in the longwave spectral range. For those domains, a mean of the low and high boundaries of the range was used in the calculation of the Spearman's correlation and the corresponding significance test.

5 Results

Each of the modifications made to CAM6.1 (described in Sect. 2.5) is relevant to the modeled dust cycle, and thus relevant to the estimate of dust climatic impacts (e.g., DRE). The proposed new (CAM6.α) and default model versions (CAM6.1) simulate a similar (Fig. S3a: relative change ∼ 16 %; CAM6.α relative to CAM6.1) global mean dust loading of 24 and 29 Tg, respectively, and DOD of 0.032 (Fig. S3c: relative change ∼ 1.3 %) (Table 6). Comparing to the recent estimates that include very coarse dust which are not included in this model, the dust loadings here are well within the range of 22–30 Tg in Kok et al. (2021a) (Table 1 of their study), and are close to the 30 Tg in Adebiyi and Kok (2020). But globally, CAM6.α shows 54 % more dust deposition than in CAM6.1 (Fig. S3b). The general spatial distributions of the relative change of dust loadings, deposition fluxes, and DOD are similar, though the magnitude of this change differs for some regions (e.g., North Africa and India).

Table 6 Simulated annual emission, loading, surface concentration, deposition, lifetime, direct radiative effect (DRE; all-sky conditions), and DREE (DRE efficiency; clear-sky conditions) of dust speciated by mineralogy and bulk-dust CAM6 with offline dynamics. The longwave DRE by dust was augmented by 51 % (Dufresne et al., 2002) to account for dust scattering which is not represented in CAM by default. Dust optical depth (DOD) and single scattering albedo (SSA) shown are for the CAM6.1 visible band centered at 0.53 µm. The global mean dust SSA was calculated over model pixels where DOD/total AOD (aerosol optical depth) > 0.5 as previously done (Scanza et al., 2015; Li et al., 2021). CAM6.α and CAM6.1 in bold represent the proposed new and default model versions, respectively. NA – not available

5.1 Evaluation of model performance and improvements on the dust cycle modeling

5.1.1 Dust emissions

To achieve the global mean DOD of ∼ 0.030, CAM6.α requires a dust emission of 2891 Tg a⁻¹ (Table 6), which falls below the estimate of 3400–9100 Tg a⁻¹ by Kok et al. (2021a; their Table 1) that accounts for dust between 0.1–20 µm in diameter and above the median, 1123 Tg a⁻¹, reported in AEROCOM phase I (Huneeus et al., 2011). The dust emission in CAM6.1 is also much lower than their estimate: 2421 Tg a⁻¹, which is however higher than the previous estimate (1490 Tg a⁻¹) with the same emission scheme (DEAD) and dust size range (< 10 µm) but using the binned method (Zender et al., 2003a).

There are no dust emission estimates from observations at a global-scale coverage. We thus infer the model performance on simulating dust emissions using model–data comparisons on the surface dust concentration and deposition flux. Note that such an evaluation of emission is probably achievable only when the observation site is close to the dust source. Otherwise, the reasoning would become incorrect, because of probable additional errors from the model representation on processes of dust transport and deposition, and interaction of dust with non-dust aerosols (e.g., sea salt and biomass burning). As will be seen in Sect. 5.1.2–5.1.4, in most of the grid cells containing the observational sites in North Africa, all experiments overestimate the deposition fluxes (Fig. 1g) and the surface dust concentrations (at Bani). This might suggest that when turning the global DOD toward ∼ 0.030, the model with the current settings and modifications probably overestimates dust emissions from North African sources, which is also shown in Kok et al. (2021a) using an integrated model ensemble and observational constraints. The smoother distribution of the dust emission in BRIFT than DEAD is due primarily to the use of the source function in DEAD that shifts dust emissions toward the most erodible soil, while in BRIFT, the near-surface friction velocity frequently exceeds the calculated threshold wind fraction velocity, causing dust to emit at more grid cells.

The locally emitted dust from the high-latitude region (> 50^∘ N and < 40^∘ S) in CAM6.α constitutes ∼ 1.6 % of the global total emitted dust flux, which is below the estimate of ∼ 5 % (2 %–3 % for each hemisphere) derived from field and satellite observations (Bullard et al., 2016; Bullard, 2017). Especially for the northern high-latitude region, where local dust sources may dominate the near surface dust concentrations (Groot Zwaaftink et al., 2016), CAM6.α substantially underestimates its contribution to the global dust (< 0.1 %). This underestimation is what we expected, since the new scheme is designed to simulate dust emissions in low-latitude regions predominantly from the impact of saltators (Kok et al., 2012), and thus may not well capture the high-latitude dust emissions which occur through different physical processes.

In comparison, despite missing dust sources > 60^∘ S (Fig. 2a), CAM6.1 may overestimate the contribution of the high-latitude dust emission to the global dust total emission (8.0 %). We attribute the much higher dust emission in the southern high-latitude region in CAM6.1 primarily to the higher emission from the South American sources (i.e., the Patagonian Desert) than in CAM6.α. This much higher dust emission is not due to local dust emissions from the Antarctic, because the local emission in the Antarctic though exists (Delmonte et al., 2013; Meinander et al., 2022; and Fig. 2b), it is weaker in strength (the contribution percentage < 0.01 %) than Patagonian deserts (Fig. 2b), and the two models (CAM6.1 and CAM6.α) also simulate a percentage contribution of dust emission from the Antarctic sources comparable to each other.

Since both dust emission schemes are far from perfect in reproducing the percentage contribution to the global dust emission and thus probably the high-latitude dust loadings, especially in the Arctic (i.e., Fig. 1e of Shi and Liu, 2019) where dust aerosol could impose big impacts on polar clouds (Shi et al., 2022), a regional tuning of the local emission in the high-latitude regions is needed to better quantify the DRE and indirect radiative effect by dust there.

5.1.2 Climatology annual means of dust optical depth, surface concentrations, and deposition fluxes

In over 90 % of the measurement sites, all models reproduce the climatology of DOD from AERONET retrievals, the surface concentration, and deposition within a factor of 10 (Figs. 1 and S4), with the spatial correlation between the models and observations statistically significant. Analysis of the spatial correlation (Pearson; R) and root mean square error (RMSE) suggests a substantial and statistically significant improvement in simulating DOD close to the source region (Fig. 1c versus Fig. 1b: R²=0.63 versus 0.41 for CAM6.α and CAM6.1, respectively, in log space; RMSE = 0.30 versus 0.41 in log space). Note we obtained the RMSE in log space which removes the dominant influence of stations with high DOD (e.g., sites in North Africa and Middle East). So the reduced bias is because the new model better captures DOD over North Africa and Australia. Compared to the improvement in DOD, the modifications do not notably better improve modeling the surface dust concentrations (Fig. 1f versus Fig. 1e: R²=0.86 versus 0.75 for CAM6.α and CAM6.1, respectively; similar RMSE ≈ 0.70 in both models) and dust deposition (Fig. 1i and h; R²: 0.78 versus 0.69 and RMSE = 0.85 versus 0.97 g m⁻² a⁻¹). This is mainly because the model's ability to simulate DOD, especially close to source regions, is subject to fewer potential errors than for surface dust concentration and deposition, which also requires the model to simulate a correct vertical distribution. Therefore, the model's ability to reproduce DOD close to the source region appears to be improved at most of the sites (33 out of a total of 36 sites; especially in Australia as shown in Fig. 1a), but this improvement does not propagate to simulations of the dust surface concentrations (Fig. 1d: improvement at 24 out of a total of 47 sites) and deposition (Fig. 1g: improvement at 62 out of a total of 108 sites).

Figure 1 Model–observation (AERONET) comparison for DOD (dust optical depth) at the visible band centered at 0.53 µm (a, b, c), dust surface concentrations (d, e, f), and surface deposition fluxes (g, h, i). Colored dots in (a), (d), and (g) show the difference between the proposed new model (CAM6.α) and observations. White symbols indicate the new model CAM6.α improves (plus sign) or worsens (minus sign) the model–observation comparison over that between the default model (CAM6.1) and observations with the metric included in the bottom right-hand corner of the figure. Numbers listed in (a), (d), and (g) are counts of the number of improved or worsened stations. The spatial correlation coefficients between model (CAM6.1: b, e, and h; CAM6.α: c, f, and i), and observations were calculated based on the annual mean values in log space (the log of each model and observational value was taken before calculating the correlation coefficient, since the values span several orders of magnitude except DOD). Dashed lines in the scatter plot show 10:1 or 1:10 lines.

Figure 2 Dust emission flux rate (kg m⁻² s⁻¹; panels a and b) and dust burdens (c, d) simulated in default CAM6.1 (panels a and c; dust emission flux rate rescaled up by 10⁹) and new model CAM6.α with the threshold gravimetric water content calculated following Fécan et al. (1999) using unity tuning factor (b=1 in panels b and d; dust emission flux rate rescaled up by 10⁸ in panel b).

Figure 3 Modeled DOD in CAM6.1 (blue) and CAM6.α (orange) in comparison with that from Ridley et al. (2016) at sub regions as defined in their Fig. 1 and from MODIS retrievals (b) at sub regions (see x axis labels). Both correlations, shown as the Kendall's τ in panel (a), are statistically significant at the 95 % confidence level. Black and gray dashed lines in panel (a) represent a factor of 2 and 4 differences.

Improvements are also seen if the climatologic DOD is compared to regional averages of the observationally constrained DOD in Ridley et al. (2016) (Fig. 3a). The new model CAM6.α substantially improves the modeled DOD, increasing the correlation (Kendall's τ coefficient) from 0.49 to 0.79, and reducing RMSE from 0.088 to 0.077, compared to CAM6.1. Spatially, CAM6.α better captures the regional DOD averaged over Australia and South Africa, which is consistent with comparison to the regional MODIS DOD (Fig. 3b). Over Taklamakan and Gobi deserts, however, the new model greatly underestimates the regional DOD compared to both estimates from Ridley et al. (2016) (Fig. 3a) and MODIS DOD (Fig. 3b; near northern China), whereas the default CAM6.1 works better, due very likely to lower dust emissions in the source regions in CAM6.α than in CAM6.1 (Fig. 2b versus Fig. 2a). Comparing with both datasets suggests that the new model may overestimate the regional DOD over North Africa and the Middle East within a factor of 2. Despite the imperfect match on the period between data and model, this overestimated regional DOD probably results from the retuning method, which provides more credits to dust emissions from North Africa and the Middle East. This overestimated DOD in the model near the source regions resulting from the tuning method may also partly explain the imperfect match between the modeled and AERONET-based DOD (Fig. 1a).

The underestimation in the surface dust concentration and overestimation in deposition occurring at several sites (near the El Djouf; near the Antarctic from our model in all cases; Fig. 1d versus Fig. 1g) is noteworthy. At some sites, such as King George in the Antarctic (62^∘ S, 58^∘ W), this phenomenon had been previously revealed by studies with multiple model ensemble mean or individual models, including an earlier version of CAM, model–data integrated study (Kok et al., 2021a), in the results of models other than CAM6.1, such as GFDL Atmospheric Model (version 2) (Li et al., 2008), and in earlier versions of CAM (Albani et al., 2014).

We suggest that the phenomenon occurs likely in part due to the following. (1) Model errors occur in simulating dust wet and dry deposition which are substantially larger than in simulating DOD and surface concentrations (Kok et al., 2021b). In addition to errors in dust emissions, and the parameterization of the dry and wet deposition schemes, MAM4 in CAM6.1 represents dust transport as an internal mixture with other species (e.g., sea salt) in the accumulation and coarse modes (Liu et al., 2016), which may have unduly increased the particle size and hygroscopicity, and thus the removal rate (dry and wet) of dust during transport to the sites (i.e., King George). (2) Dust sources are possibly misrepresented in the Southern Hemisphere in the model. With current emission sources, the increase of the emission rate with BRIFT from Patagonia compared to DEAD slightly mediates the underestimation of dust surface concentration at King George. A further increase of the dust emission may help reduce the underestimation of dust deposition in land and the surface concentration at King George, but it would then exacerbate the bias in simulating the surface deposition at that site. (3) The limited observation period could result in the climatology representative issue, considering the episodic character of dust events. This limitation due to observation period may be particularly important for observed dust in the Southern Hemisphere where the dust quantities tend to be more episodic than in the Northern Hemisphere (Mahowald et al., 2011b).

Figure 4 Modeling performance for the seasonal cycle of DOD (a) and dust surface concentrations (b) by CAM6.α (new model) against in situ (site names listed in the figure) measurements, relative to the performance of CAM6.1 (default model) against in situ measurements (improvement, degradation, no change indicated by “+”, “−”, and none characters, respectively). Colored left and right semi-circles represent Kendell's τ coefficient between CAM6.α and observations and the ratio of standard deviation (CAM6.α over observation; each normalized by annual mean values), respectively.

As to the relative importance of dry and wet deposition, we find that the dust wet deposition may dominate the total deposition of dust in the remote oceanic area (Fig. S5a; Table 7), and thus affects the long-range dust transport. The models tend to overestimate the observed percentages of the wet deposition (Table 7). This overestimation could be due partly to the internal mixing assumption of dust aerosols with sea salts which increases hygroscopicity of the aerosol mixture during transport. Correcting the coarse-mode distribution, following we suggest (Table 1), does not help improve the model performance (Table 7). With that said, a recent study had shown that the CMIP6 models overestimate the precipitation frequency, particularly for the light precipitation (0.1–20 mm d⁻¹) (Na et al., 2020). It could be the same reason – the unduly simulated precipitation frequency in CAM6.1 – that explains the overestimated importance of wet deposition compared to the observations we have here. Therefore, future model changes on the cloud physics that reduce the light precipitation frequency may help better simulate the transport of dust aerosols across zones where frequent precipitation occurs (e.g., the intertropical convergence zone (ITCZ)). Considering the limited observations on the partitioning of the dust total deposition between dry and wet processes however, we cannot draw a concrete conclusion that CAM6.1 overestimates the wet dust deposition fluxes (Table 7).

Table 7 Percentage (%) of wet deposition. Observations compiled by Mahowald et al. (2011b) from data at Bermuda (Jickells et al., 1998), Amsterdam Island, Cap Ferrat, Enewetak Atoll (Arimoto et al., 1985), Samoa; New Zealand sites (Arimoto et al., 1990); North Pacific sites (Uematsu et al., 1985); Greenland dye 3 (Hillamo et al., 1993), coastal Antarctica (Wagenbach et al., 1998), and Dome C of Antarctica (Wolff et al., 2006). RMSE: root mean square error; R: Spearman's correlation.

^a shown are minimum and maximum of the annual wet percent among the four sites. ^b Non sea salt-sulfate.

5.1.3 Seasonal cycle of climatology dust optical depth and surface dust concentrations

Dust optical depth

Both CAM6.1 and CAM6.α reasonably reproduce the retrieved seasonal cycle at the selected AERONET sites except Ilorin (Fig. S6), where both models greatly underestimate the observed DOD in winter (Fig. S6b). It is possible that non-dust aerosols (e.g., black carbon) transported from South Africa contaminated the observation, leading to an artificially high DOD during the winter season at that site.

The new model CAM6.α improves both the temporal correlation based on the monthly values and standard deviation, compared to CAM6.1, only at three (Ilorin, Dhabi, and Tamanrasset_INM; Fig. 4a) out of a total of the 11 selected AERONET sites where the measurements cover the whole 12 months in a year (this reduced the number of the total sites, compared to that used in the climatology comparison). Significant improvements on the modeled seasonal cycle of DOD occurs at Tamanrasset_INM (25^∘ N, 4^∘ E). CAM6.α increases the temporal correlation coefficient from 0.42 to 0.82 (Fig. S6e). Despite the improvement, the new model continues largely overestimating the observed DOD at this site, especially in the peak month of June (Fig. S6e), resulting in an overestimated annual mean DOD.

Similar results are obtained if the seasonal cycle of DOD is compared to model–data constraints on regional DOD in Ridley et al. (2016) (Fig. S7): spatial correlation analysis on the seasonal mean DOD suggests that the new model CAM6.α substantially improves the modeled DOD in all seasons (Fig. S7) with higher correlations that are statistically significant (Fig. S7), compared to simulations using CAM6.1 (e.g., in JJA CAM6.α: R = 0.71 versus CAM6.1: R = 0.48), though CAM6.α does not reduce the RMSEs in all the seasons.

5.1.4 Surface dust concentrations

In terms of temporal correlation and standard deviation for assessing the seasonal cycle, the modifications do not uniformly improve the model performance on reproducing the surface dust concentration (Fig. 4b). Only at 7 out of the 19 sites in total (a reduced number of the total sites, compared to that used in the climatology comparison, due to the removal of sites where there is no full coverage of the measurement over the 12 months in a year) – Bani and Cinz in North Africa, Mace Head in the North Atlantic, Cape Grim in Australia, Funafuti and Norfolk in the Eastern Pacific, and Hawaii in the North Pacific (Fig. 4b) – the modifications result in improvements using both metrics. Examining a third metric, the difference between modeled and observed surface concentration in specific months, we have 13 of the 19 sites where at least half a year shows improvement. Still, the new model overestimates the surface concentration of dust at many of those 13 and the other sites during most months in the year (Figs. S8–S9). This overestimation is particularly pronounced for Cape Verde, likely mainly because of the strong dust emission in western North Africa using BRIFT compared to DEAD. The new model produces significant improvement in terms of all the three metrics at Bani (14^∘ N, 3^∘ E; Fig. S9i), increasing the temporal correlation from 0.21 (insignificant at the 95 % confidence level) between CAM6.1 and the observations to 0.58 (significant at the same confidence level) between CAM6.α and the observations.

5.1.5 Size distribution of transported dust

Figure 5 shows the simulated size-resolved dust mass compared to AERONET retrievals and in situ measurements. In general, the new model CAM6.α with the mode size distribution from CAM5 better reproduces the retrieved atmospheric size distribution than the default CAM6.1 with the size distribution from CAM6.1 over most sites. At only three sites (La Laguna: 28^∘ N, 17^∘ W; Santa Cruz Tenerife: 29^∘ N, 16^∘ W; and Puerto Rico: 18^∘ N, 67^∘ W) the mass size distribution from CAM6.α becomes evidently worse than from CAM6.1. Compared to CAM6.α, CAM6.1 tends to carry more dust in mass with the diameter > ∼ 5.0 µm, which also overshoots AERONET retrievals in that size range (Fig. 5). This bias in CAM6.1 could be lower for mass of dust > ∼ 5.0 µm, considering that AERONET retrievals might have a bias towards fine dust when compared to in situ measurements (McConnell et al., 2008).

Figure 5 Modeled and observed atmospheric size-resolved dust mass in the geometric diameter range of 1–10 µm at AERONET stations. Numbers in each plot indicate the Kendall's τ coefficient between model and observations (blue bars). The model runs here include the one using the old model with the mode size parameters from CAM6 by default (CAM6.1 in cyan) and the other one using the new model with the mode size parameters from CAM5 (CAM6.α in black). In (a), only results from CAM6.1 can be clearly seen because the two lines are overlapping.

Figure 6 Normalized size distribution of dust between 0.2 and 10 µm diameter in the global average (a), near Canary Island (blue colors in b; dot: 2.5 km; cross: 6–7 km; data for June/July 1997 from Otto et al., 2007), and near Cabo Verde (orange colors in c; dot: 2.5 km; cross: 6–7 km; data for August 2015 taken from Ryder et al., 2018). The default model, CAM6.1: (purple line); the new model, CAM6.α: (red line); semi-observations: DustCOMM (black line) inverted based on an integration of a global model ensemble and quality-controlled observational constraints on the transported dust size distribution, extinction efficiency, and regional DOD with data taken from Adebiyi et al. (2020). We chose the model layers and grid cells that are closest to the location and atmospheric height, as well as the months, where and when the measurements were made for comparison.

When comparing global mean model results to those from DustCOMM (Adebiyi et al., 2020) (Fig. 6a), generally CAM6.α better reproduces the atmospheric size distribution (d$V/$dlnD) than most of the other models (e.g., WRF-Chem: Weather Research Forecasting-Chemistry) (Adebiyi et al., 2020) in the full size range, and CAM6.α for dust < 2 µm in diameter (Fig. 6a). Like most climate models shown (Adebiyi et al., 2020), CAM6.α tends to underestimate coarse dust with the diameter greater than ∼ 5 µm (our model currently excludes dust with diameter > 10 µm), which is consistent with the finding suggested by the comparison between CM6.1 and AERONET observations (Fig. 5). The size distribution from CAM6.1 compares well with the DustCOMM result for dust between 1–2 µm in diameter. But it greatly underestimates the fine dust fraction (diameter < 2 µm) which CAM6.α can better capture due primarily to the more correct gravitational settling velocity modeled by using the new dry deposition scheme.

We then evaluate the model's performance in reproducing the size distribution measurements at the high-atmosphere levels (2–5 and 6–7 km) near the Canary Island (Fig. 6b) and Cabo Verde (Fig. 6c) by Ryder et al. (2018) and Otto et al. (2007) for transported dust (for the model–observation comparison at other atmospheric levels, see Fig. S10). Overall, CAM6.α better reproduces the size distribution at the higher atmospheric level (6–7 km) than CAM6.1, but CAM6.α substantially overestimates the mass of dust < 5 µm at the lower atmospheric level (2–5 km) compared to the measurements where CAM6.1 performs better. As also suggested in the global size distribution comparison, CAM6.1 simulates more dust > 5 µm and less dust < 5 µm than CAM6.α. However, both models underestimate the observed mass fraction of dust > 5 µm at the high-atmosphere level (6–7 km) near the Cabo Verde. The models also fail to capture the change in the size distribution between the two atmospheric levels that the measurements suggest. It is worth noting that the measurements are from single campaigns or flights that may have representative issues not reflecting the climatological size and vertical distributions of dust aerosols (i.e., limited by the space and time coverage), and the model results do not perfectly match the periods and locations of the measurements.

5.2 Impacts of each modification on the dust cycle modeling

This section details the relative importance of each modification to the modeled dust properties (loading and/or other dust variables). We show in Sect. 5.2.1–5.2.4 the results on the global mean and spatial distribution, and in Sect. 5.2.5 how the modifications affect the dust properties on the regional mean basis.

5.2.1 Dust emission schemes: BRIFT versus DEAD

The dust emission in MINE_NEW_EMIS using BRIFT (2910 Tg a⁻¹) is 35 % lower than in MINE_BASE using DEAD (4456 Tg a⁻¹), consistent with the lower DOD (0.035 versus 0.029), due primarily to higher dust lifetime in the former (3.1 and 2.2 d) (Table 6). The relative strength of dust emissions for different sources also differs between DEAD and BRIFT, as Kok et al. (2011b) documented based on CAM4. The comparison between the two emission schemes in CAM6 here on the spatial distribution of the dust emission largely remains as in CAM4. For example, the preferential source function of Zender et al. (2003b) used in DEAD simulates most of the emission in the central part of North Africa (e.g., the Bodélé depression) (Fig. 2a). In comparison, the dust emission coefficient in BRIFT (Eq. 10) and the new method of calculating the threshold gravimetric water content of the topsoil layer (Eq. 9; see values for the tuning factor “b” in Table 2) shifts the main dust emission in North African source westward and southward into the dust source belt (Kok et al., 2011b). This shifting in BRIFT, compared to DEAD, tends to have the dust emission occur in the wind erodible areas that satellite-based retrievals suggest (Ashpole and Washington, 2013; Ginoux et al., 2012), though the retrieval of dust beneath clouds are unavailable which may lead to a missing of potential dust sources that satellite retrievals cannot detect, for example, dust emissions occurring at the presence of deep convection (Engelstaedter and Washington, 2007; Marsham et al., 2013). The much lower dust emission in Taklamakan and Gobi deserts in China relative to that from North Africa using BRIFT, concerning that comparison using DEAD (Fig. 2b), is likely due to the high soil moisture simulated in CAM6 in both deserts in China.

Another pronounced difference in the modeled dust emission occurs in less erodible areas (i.e., North America, South Africa, Australia), where BRIFT tends to decrease the emission flux compared to using DEAD, an opposite response than that simulated for the North African sources. Such as in Australia, both schemes simulate the maximum in dust emissions from the Great Artesian Basin and the Murray–Darling Basin, but BRIFT reduces the dust emissions there, bringing a better agreement on the climatological DOD with AERONET observations than DEAD. However, BRIFT, using the unity tuning factor to calculate the threshold gravimetric water content, simulates high dust emissions in western Australia instead of central and eastern Australia as previously documented (Ginoux et al., 2012). Sensitivity tests suggest that using inversed clay fraction can likely better capture the spatial emission pattern in Australia (Fig. S11). In Patagonia, as Kok et al. (2014b) found based on CAM4 simulations, using BRIFT in CAM6.1 substantially increases the dust emission compared to DEAD. In addition, BRIFT simulates the dust emission from a source in northern Chile (the Atacama Desert) and the high-latitude area, where no dust emits in DEAD.

Due to the southwestward shifting of dust emissions in BRIFT to the “real” dust belt in North Africa (Sect. 5.2.1), dust aerosol particles experience stronger vertical transport by near-surface convergence that controls the annual cycle of North African dust (Engelstaedter and Washington, 2007). The lifetime of dust thus tends to be higher for aerosol particles experiencing strong convection, which uplifts them high above the surface (Cakmur et al., 2004), increasing the dust lifetime from 2.2 d in MINE_BASE to 3.1 d in MINE_NEW_EMIS. This lifetime-changing mechanism in turn indicates the importance of accurately simulating convergence-related convection (i.e., haboob) (Marsham et al., 2011) and where the dust emission occurs for dust transport modeling, especially the cross-Atlantic/Pacific (Prospero, 1999; Prospero et al., 2020) and cross-equatorial transport (Kok et al., 2021a; Li et al., 2008), which currently the models do not well represent.

In response to the change in dust emissions due to shifting from DEAD to BRIFT, the global annual mean dust deposition and loadings decrease by 35 % and 7 %, respectively (Table 6). Considering the lower global DOD in BRIFT than in DEAD (0.035 versus 0.029), differences between the global annual mean dust deposition in BRIFT and DEAD would become smaller, if we rescaled the global annual mean dust deposition and loadings offline using factors to make the global mean DOD in the two experiments exactly equal 0.030. The change in the total dust loadings (Fig. 7a) and deposition (Fig. 7b) has a similar spatial distribution: a great increase occurs primarily in the Southern Ocean, the Middle East, and the western Atlantic Ocean, western Australia and its downwind areas, and a great decrease in the Pacific Ocean due primarily to reduced dust emissions in East and Central Asia (Fig. 2b versus Fig. 2a); near Greenland, BRIFT simulates more dust deposition and slightly less dust loadings, owing to the local dust emission that occurs in BRIFT (Fig. 2b) but not in DEAD (Fig. 2a), and the ability of transporting further for dust aerosols in BRIFT because of the increased lifetime (Table 6).

Figure 7 Impacts of the dust emission scheme (a and b: ratio of BRIFT to DEAD), aerosol dry deposition scheme (c–f: ratio of PZ10 to Z01), and dust shape (g and h: ratio of ellipsoidal to spherical dust) on the modeled dust deposition (total: b, d, and h; fine mode: c), and dust loading (total: a, f, and g; fine mode: e). The Taylor diagram (i) compares dust loading in 21 sub-regions defined in Fig. S16. In panels (i), DEAD shows comparison between MINE_BASE and MINE_NEW_EMIS; Asp+PZ10 between CAM6.α_MINE and MINE_NEW_EMIS; PZ10 between CAM6.α_MINE and MINE_NEW_EMIS_SHAPE; the Kendall's τ temporal correlation and the standard deviation were obtained based on monthly values with the seasonal cycle removed.

Interestingly, we also find considerable changes to the simulated mass fraction of dust minerals between using BRIFT and DEAD (Fig. S12). These changes are as expected, given the redistributed “hot spots” where dust emission occurs by switching to BRIFT (Fig. 2) and the grid-dependency of the soil mineralogy that we used to initialize the dust speciation (Fig. 1 of Scanza et al., 2015 and Fig. S2 of Li et al., 2021). This change in the simulated mineral mass fractions of dust matters for quantifying the dust shortwave DRE (e.g., hematite) at the top of the atmosphere (Sokolik and Toon, 1999; Balkanski et al., 2007; Li and Sokolik, 2018b; Li et al., 2021), the cloud–aerosol interaction (e.g., feldspar) (Atkinson et al., 2013), and biogeochemistry effect (e.g., irons) (Mahowald et al., 2011a). It thus deserves quantifying how the shift of dust sources changes the simulated mineral content of iron-bearing minerals, including hematite and illite, and feldspar in the dust.

The results suggest that BRIFT simulates 10 times more hematite than DEAD in terms of mass (or volume) fraction in the Northern Hemisphere (BRIFT: 1.0 %; DEAD: 0.098 %) and 25 % less in the Southern Hemisphere (BRIFT: 1.2 %; DEAD: 1.6 %). Such a decreasing of the simulated mass fraction of hematite aerosol in the Southern Hemisphere is due primarily to reduced dust emissions from the Australian deserts, the soil of which enriches iron oxides (Claquin et al., 1999; Journet et al., 2014). BRIFT also shifts the dust emission westward in Australia (Fig. 2b versus Fig. 2a) where the soil abundance of hematite is lower than in the Australian deserts. Similarly, the increased mass fraction of hematite aerosol in the Northern Hemisphere can be partially attributed to the reduced dust emission from East Asia. The change is also evident (increase with the relative change > 30 %) to feldspar in the South Ocean (Fig. S12d) and to calcite in the North Pacific Ocean (decreasing; Fig. S12f), which may have implications for the amount of the ice nucleation by mineral dust since this nucleation could be dominated by feldspar in mixed-phase clouds (Atkinson et al., 2013).

5.2.2 Dust deposition schemes: PZ10 vs. Z01

The comparison of the dry deposition velocity between PZ10 and Z01 is size dependent. Because of the reduced dry deposition velocity in the fine mode, moving to PZ10 from Z01 greatly decreases the dry deposition of fine-mode (Aitken plus accumulation) dust within the low-to-mid latitude regions (between 40^∘ S and 40^∘ N; Fig. 7c; PZ10 : Z01 < 0.3; similar in the accumulation mode only as Fig. S13b suggests). Since most dust mass is in the coarse mode, the small change of dust deposition in this mode, because of the slightly larger dry deposition velocity in the coarse mode in PZ10 than in Z01, results in a slight change in the total dust deposition in the low-to-mid latitude regions (Fig. 7d or Fig. S14b). Even for dust deposition in the fine mode, the increased wet deposition by using PZ10 (such as in the accumulation mode shown in Fig. S15b) offsets the reduced dry deposition in the low-to-mid latitude regions, resulting in a negligible change spatially and on global average. In the South Ocean (downwind of the Patagonian deserts), a decrease of dry deposition fluxes causes more fine-mode dust aerosol particles near the source regions, which then become cloud borne. The increased cloud-borne particles in turn increase the possibility of horizontal transport and release of particles by the cloud droplet evaporation, leading to an increase of the dry deposition flux at the downwind regions (Fig. 7c). But the reduced total (dry plus wet) dust deposition in the coarse mode dominates over the increased dry deposition flux at the downwind regions, leading to a considerable decrease of total dust deposition by > 30 % (relative change; Fig. 7d).

Compared to Z01, PZ10 increases the global mean dust loading in the fine mode by ∼ 20 % (Fig. 7e). Particularly in the tropics, such an increase in the remote areas can be over 60 %, though the dust abundance there is low. The slight decrease of dust in the coarse mode dominates the change in the total dust loading, resulting in a slight decrease of the global mean total dust loading by 6 % (Fig. 7f). Correspondingly, the global mean DOD remains almost the same between the simulations using PZ10 and Z01 (Table 6).

5.2.3 Dust asphericity

The overall change to the spatial distribution of dry deposition induced by dust asphericity is not as important as the change induced by changing to the dry deposition scheme PZ10. The models simulate similar overall spatial distributions of dust deposition at the surface between modeling dust as spherical and ellipsoid shaped particles (Fig. 7h). The lower gravitational settling velocity when modeling dust as ellipsoids induces a considerable change to dust deposition only locally within remote areas: the South Pacific, western and eastern equatorial Pacific, and downwind of Patagonia has an increase of the dry deposition by up to 30 % (MINE_NEW_ EMIS_SHAPE versus MINE_NEW_EMIS). In comparison, little change to the dust deposition by dust asphericity occurs near/over major dust source regions. This contrast in the changes in the dry deposition flux between close-to-source and remote areas suggests that including dust asphericity could potentially mediate the overestimated dust emission from source regions (e.g., North Africa), because dust asphericity could enlengthen the lifetime in the atmosphere, and thus it takes less amount of dust to have the same amount of dust loadings and DOD as spherical shape assumption does.

5.2.4 Dust size representation

The removal rates of dust aerosol particles by both dry and wet deposition highly depends on their size (Mahowald et al., 2014). Since most of dust loadings are in the coarse mode, changing parameters of the coarse-mode size distribution (σ, initialized GMD, and the prescribed minimal and maximum boundaries within which the modeled GMD can vary, Table 1) from σ=1.2 to 1.8 halves the lifetime of dust (lifetime = 4.9 d versus 2.4 d; Table 6). This reduction of dust lifetime is primarily due to the change in σ of the coarse mode rather than the initialized GMD and its boundaries, as we obtain almost the same dust lifetime (∼ 2.4 d) between experiments with different parameters for dust size distribution but identical σ=1.8 (NEW_EMIS_SIZE versus NEW_EMIS_SIZE_WIDTH; Table 6).

We also notice a different DOD simulated by NEW_EMIS_SIZE (DOD = 0.013) and NEW_EMIS_SIZE_WIDTH (DOD = 0.019). The prescribed GMD boundaries do not affect the simulated dust loadings and DOD, because the predicted GMD in the model varies little. We can therefore derive that the initialized GMD itself is also relevant to simulated DOD, but its influence (relative change = 20 %) is second to that of changing the coarse-mode σ. Thus, it is the increased σ of the coarse mode that explains the reduced dust loadings (22 versus 11 Tg in NEW_EMIS and NEW_EMIS_SIZE, respectively; Table 6; Fig. 8b) and DOD (0.030 versus 0.013 in NEW_EMIS and NEW_EMIS_SIZE, respectively; Table 6). This impact of changing the coarse-mode σ is also greater than that of the other modifications (e.g., speciating dust or changing the dust emission scheme from DEAD to BRIFT) on the simulated dust lifetime which appears trivial (e.g., dust lifetime increased by 0.6 d only by changing to the new emission scheme). Correspondingly, given a similar emission rate, changing the coarse-mode σ affects DOD most, compared to the other modifications we made.

Figure 8 Impact of changing the coarse-mode geometric standard deviation (σ) for transported dust aerosol on the modeled dust surface deposition fluxes and column loading: ratio of NEW_EMIS_SIZE (S5: σ=1.8) to NEW_EMIS (S6: σ=1.2) (see Table 2 for case names). Numbers on the top of the plot show ratios on global average.

5.2.5 Impacts of the modifications on the regional mean basis

The regional analysis over 21 selected sub-regions (Fig. S16 for definition) suggests that over most of those sub-regions, the simulated dust loading/deposition flux using the model under different modifications (PZ10 and/or dust asphericity) except to the dust emission scheme closely correlates (temporal correlation coefficient > 0.85 based on monthly values) with that in the reference case MINE_NEW_EMIS (Fig. 7i). In addition to slightly increasing dust loading, introducing dust asphericity to the model slightly increases the temporal variability of the modeled dust loading, while replacing Z01 with PZ10 slightly decreases the variability of the simulated dust loading with respect to the reference case generally in nearly all the 21 sub-regions (Fig. 7i). The combined effect of the two modifications on this temporal variability is more determined by the choice between PZ10 and Z01 than dust asphericity.

Using different dust emission schemes changes the regional dust loading/deposition flux the most among those modifications in terms of the standard deviation (the BRIFT $/$ DEAD ratio > 1.25 or < 0.75 in many regions) and temporal correlation (low-to-moderate temporal correlation between 0.15 and 0.85) (Fig. 7i). Particularly, the strong regional contrast on the dust loading/deposition exists in northwest Asia (region 5: the BRIFT $/$ DEAD ratio < 0.5 and temporal correlation < 0.5), Central Asia (region 6: the BRIFT $/$ DEAD ratio < 0.6 and temporal correlation < 0.5), southeastern Pacific Ocean (region 4: the BRIFT $/$ DEAD ratio > 1.5 and temporal correlation = ∼ 0.2), and southern America (region 21: the BRIFT $/$ DEAD ratio > 1.5 and temporal correlation < 0.5).

5.3 Dust direct radiative effect

CAM6.α yields a global mean net dust DRE of ∼ −0.26 W m⁻² (shortwave plus longwave; longwave has been augmented by 51 % to include dust scattering), which is slightly less cooling than in CAM6.1 (∼ −0.30 W m⁻²). But the net dust DRE can strongly differ between the two model versions at regional scales (Fig. S17b). For example, CAM6.α suggests more warming (difference > 2 W m⁻² in amplitude) near Australia due to reduced dust loadings (thus DOD) (Fig. 7a) and hematite mass fraction (Fig. S12a), and more cooling (difference > 2 W m⁻² in amplitude) in downwind regions of North Africa primarily due to increased dust loadings (Fig. 7a). The opposite change in one region relative to another however cancels out, resulting in a negligible net DRE change at the global scale (−0.04 W m⁻²). The following subsections evaluate the model performance on reproducing the observed dust DRE efficiency (Sect. 5.3.1), and quantify the impact of each modification on the estimate of dust DRE (Sect. 5.3.2 and 5.3.3).

5.3.1 Dust direct radiative effect efficiency

All model versions as shown in Fig. 9 have difficulty in reproducing the dust DRE efficiency under clear-sky conditions. In the shortwave spectral range (Fig. 9a), the new model, CAM6.α, does not show improvement in general. It works better in reproducing the retrievals only in the Atlantic Ocean (tropical Atlantic-3: 10–30^∘ N, 20–45^∘ W) in the summer and at a site in the Mediterranean Basin (33.5^∘ N, 12.6^∘ W) in September. In the longwave spectral range (Fig. 9b), the dust DRE efficiency in the new model, CAM6-α, agrees better with retrievals than that in CAM6.1 with higher correlation coefficients (CAM6.α: R = 0.86 versus CAM6.1: R = 0.60) and reduced RMSEs (CAM6.α: RMSE = 4.0 W m⁻² DOD⁻¹ versus CAM6.1: R = 5.7 W m⁻² DOD⁻¹), likely mainly owing to the improved representation of the dust cycle. It is worth noting that, in addition to uncertainty due to the imperfect representation of the spatial distribution of dust aerosols (Fig. 1), the different spectral ranges in the model and the satellite-based sensors and radiation parameterization in the model (Jones et al., 2017) may also contribute to the difference between dust DRE efficiency from the model and observations.

All the modifications do not change the global mean net DRE efficiency (Table 6; quantified using the mineralogy runs), except that BRIFT yields the global mean net efficiency value that substantially differs in the shortwave spectral range compared to DEAD (MINE_NEW_EMIS: −2.1 W m⁻² DOD⁻¹ versus MINE_BASE: −4.0 W m⁻² DOD⁻¹; Table 6).

Figure 9 Modeled and observed dust direct radiative effect efficiency (DREE) in the shortwave (SW; a) and longwave (LW; b) spectral ranges under clear-sky conditions at the TOA over the sub-domains (see labels of the x axes) in April–June (AMJ), summer (JJA), fall (NDJ), and September (Sep) for the 2000s climate. The DREE is defined as the ratio of DRE to DOD, so has units of W m⁻² DOD⁻¹. Included cases from left are CAM6.1, CAM6.α, MINE_NEW_EMIS_SHAPE, and CAM6.α_MINE. Colored numbers show correlation coefficient (R) and the root mean square error (RMSE) between the model and retrievals in the SW and LW spectral ranges or in both spectral ranges (numbers in parenthesis in panel a).

5.3.2 Impacts of dust asphericity, dry deposition scheme, and dust emission scheme

The dust asphericity introduces negligible (relative change < 10 %) impacts on the global net dust DRE, and PZ10 enhances the net dust cooling by ∼ 18 % relative to that using Z01 (Table 6). Regionally, the slightly higher/lower dust loading or DOD due to dust asphericity only slightly enhances the warming over land (Fig. 10a; e.g., North African land; net DRE: 0.97 and 1.1 W m⁻² for MINE_NEW_EMIS and MINE_NEW_EMIS_SHAPE, respectively; the single scattering albedo at the visible band ∼ 0.90 for both runs, not shown) and enhances the cooling over ocean (e.g., downwind of North Africa). PZ10 simulates a slightly enhanced cooling relative to Z01 almost everywhere (Fig. 10b; e.g., south northern Atlantic Ocean, net DRE: −0.72 and −0.76 W m⁻² for MINE_NEW_EMIS_SHAPE and CAM6.α_MINE, respectively).

Figure 10 Spatial distribution of the net (shortwave plus longwave) direct radiative effect (DRE) difference (unit: W m⁻²) at the top of the atmosphere under all-sky conditions in current climate between model results using non-spherical (asp) and aspherical dust (sph) (a), PZ10 and Z01 (b), Asp+PZ10 and Sph+Z01 (c), and BRIFT and DEAD (d). The longwave DRE was augmented by 51 % to account for the dust scattering. Numbers shown in each panel title represent annual mean difference in global average.

Calculations suggest a regionally strongly contrasted change to net dust DRE when shifting from DEAD to BRIFT (Fig. 10d), but the enhanced cooling in one region (i.e., the downwind Atlantic Ocean of North Africa: BRIFT: −0.76 W m⁻²; DEAD: −0.64 W m⁻²) and warming in another (i.e., western Africa) cancel out, resulting in a weaker global dust cooling, −0.08 W m⁻² (Table 6). These regional dust DRE differences primarily result from the regional changes to DOD/dust loadings in response to the spatial change in dust emissions, especially for North African sources, and are due partially to the change in hematite mass fractions (Fig. S12a).

5.3.3 Sensitivity to the size distribution

In NEW_EMIS_SIZE, the dust DRE at the shortwave bands at the top of the atmosphere under all-sky conditions is ∼ −0.39 W m⁻² (Table 6). In contrast, NEW_EMIS yields approximately 70 % and 62 % stronger cooling effects of −0.66 W m⁻² by dust aerosols. We attribute this strong shortwave cooling in NEW_EMIS primarily to the greatly overestimated mass fraction of fine dust, which is more scattering than coarse dust. The other parameters, such as the GMD bounds of the coarse mode is also relevant to the shortwave dust DRE calculation, inducing a change of 0.12 W m⁻² (NEW_EMIS_SIZE minus NEW_EMIS_SIZE_WIDTH), which is only slightly smaller than 0.15 W m⁻² (NEW_EMIS_SIZE_WIDTH minus NEW_ EMIS) that results from the σ change from 1.2 to 1.8 (Table 6). Compared to its influence at the shortwave bands, the size change only slightly affects the longwave dust DRE calculation (relative change < 30 %).

Spatially, big differences (less cooling; absolute difference > 3.5 W m⁻²) on shortwave dust DRE caused by the size change (from CAM6 size to CAM5 size) mainly appear over areas close to the non-reflective regions (e.g., dust source regions where annual surface albedo at visible band < ∼ 0.2 and oceanic regions adjacent to North Africa and the Middle East) (Fig. S18a). The coarse mode size change from CAM6 size to CAM5 size systematically reduces the longwave warming over all grid cells (Fig. S18c) primarily due to the σ change, as the other parameters enhance the warming effect instead (Fig. S18d).

5.4 Relative importance of each modification

Figure 11 compares the relative importance of each modification on the modeled dust quantities, and the dust DRE at grid cell scales and on the global average. Overall, replacing the size distribution of dust aerosol and the dust emission scheme with new ones are more influential on the modeled quantities of dust (DOD, burden, and deposition) and the DRE estimate, compared to the other modifications. At model grid cell scales, this is especially true for close-to-source regions: the size change dominates over all the others to be the most important factor in modeling the surface dust concentration which occurs everywhere (Fig. 11a), modeling DOD (Fig. 11c), and estimating the dust lifetime (Fig. 11e) and DRE (Fig. 11f) which occurs at most cells; the choice of the dust emission scheme could be most important in modeling the dust burden (Fig. 11b) and deposition (Fig. 11d). Dust asphericity can only dominate the change to the modeled dust burdens or loadings in the South Pacific Ocean (Fig. 11b), where the dust mass is low relative to close-to-source regions. As for the dry deposition scheme, switching to PZ10 can dominate the change to the dust lifetime at the north polar region (Fig. 11e) where the total dust is more in the fine mode for which PZ10 reduces the dry deposition velocity (Petroff and Zhang, 2010).

On the global average (Fig. 11g–l), the size change is most important in modeling most of the dust quantities, except deposition (Fig. 11j) and in estimating the dust DRE at the top of the atmosphere (Fig. 11l) for which the choice of the dust emission scheme becomes more influential.

Figure 11 Summary of the relative importance of the modifications (spherical, sph versus non-spherical dust, asp; default dry deposition scheme, Z01 versus new, PZ10; DEAD versus BRIFT dust emission scheme; and coarse-mode size distribution used in CAM5, new size, versus that used in CAM6.1, old size) at grid cell levels (a–f) and in global average (g–i) on surface concentration (a, g), burden (b, h), DOD (c, i), surface deposition (d, j), lifetime (e, k), and net DRE (f, i) from simulations with offline dynamics.

6 Bulk- versus speciated-dust model

The bulk- (CAM6.α) and speciated-dust (CAM6.α_MINE) models simulate a similar dust cycle with the difference between the two types of models orders of magnitude smaller than the dust cycle itself modeled either by CAM6.α or CAM6.α_MINE (e.g., Figs. 12 and 13). This similarity results from several factors including:

1. tuning the dust cycle to a global mean DOD of 0.03;

2. nudging both models towards the same meteorology dynamics; and

3. conserving the dust mass when speciating the dust aerosols, such that summing the mass fraction of each dust species equals unity. For the same reasons, the influence of each of the modifications on the modeled dust cycle quantified using the bulk-dust model instead of the speciated-dust model, as this study used, would be similarly comparable.

What differs remarkably is the modeled dust optical properties between the speciated- and bulk-dust simulations. For example, the speciated-dust model (CAM6.α_MINE) yields a lower global-mean dust SSA than the bulk-dust model (CAM6.α): 0.896 versus 0.911 (Table 6), at the visible band centered at 0.53 µm. Note that the dust DRE is sensitive to variation of the dust SSA. This lower dust SSA obtained here in the speciated-dust model than in the bulk-dust model is consistent with the finding of a previous study (Scanza et al., 2015) using an earlier model version (CAM5). Correspondingly, CAM6.α_MINE yields a reduced dust cooling (Table 6) and DRE efficiency (Fig. 9) relative to CAM6.α.

For dust DRE efficiency (Fig. 9), speciating dust in CAM6 tends to reduce the RMSE while retaining the horizontal spatial correlation in either shortwave (CAM6.α: RMSE = 11 W m⁻² DOD⁻¹; R = 0.26 versus CAM6.α_MINE: RMSE = 10 W m⁻² DOD⁻¹; R = 0.20) or longwave (CAM6.α: RMSE = 4.0 W m⁻² DOD⁻¹; R = 0.86 versus CAM6.α_MINE: RMSE = 3.0 W m⁻² DOD⁻¹; R = 0.84) or both spectral ranges (CAM6.α: RMSE = 7.0 W m⁻² DOD⁻¹; R = 0.93 versus CAM6.α_MINE: RMSE = 6.0 W m⁻² DOD⁻¹; R = 0.92). This comparison suggests that modeling dust as component minerals with the coarse-mode dust size distribution of MINE_NEW_EMIS_SIZE could help improve the model performance relative to modeling dust as a bulk to reproduce the retrieved dust DRE efficiency (Fig. 9).

It is worth noting that the improvement in modeling the dust DRE efficiency could be artificial because of the combined use of imaginary part of the complex refractive index of hematite (see Fig. 1b of Li et al., 2021) and the volume mixing rule used in the dust speciated model to compute the bulk-dust complex refractive index, leading to artificially more absorptive dust than in the bulk-dust model (Fig. 9a and Table 6).

Figure 12 Surface dust emissions (a; global annual mean = 2891 Tg) and deposition fluxes (b; global annual mean = 2893 Tg) simulated by CAM6.α and their differences (c and d) between CAM6.α_MINE and CAM6.α.

Figure 13 The same as Fig. 12 but for DOD (a: global annual mean = 0.030 and (c) global mean difference = 0.001) and dust burdens (b: global annual mean of dust mass = 24 Tg and (d) global mean difference ≈ 0 Tg), respectively.

7 Limitations in the model–observation comparison

The following are issues which may affect the model–observation comparison when interpreting the comparison:

1. The period when the measurements were made not perfectly matching with when the simulations were performed for.

2. Different representative space volume between the model results and observations. The model results are representative of a colocation in space which is determined by the spatial resolution and often too large compared to the volume that observations represent (Hamilton et al., 2019; Wang et al., 2014). Ground stations measure dust-related quantities using stationary instruments, and aircraft-onboard instrument measures dust along with the flight tracks.

3. Different size cut-offs. Some observations include dust of size > 10 µm in diameter (between 10–20 µm; dust particles in this size range are also present over the source regions and regions downwind of North Africa as Ryder et al. found in 2019, but nearly all the observational constraints used in this study do not include those dust particles) which our models do not simulate. This might be an important error source (Adebiyi and Kok, 2020). On the other hand, the observations of PM₁₀ are likely to include only PM_6.9, because what is measured is in aerodynamic not geometric diameters (Huang et al., 2021; J. S. Reid et al., 2003). Finally, the modeled dust mass is for dust with our own defined mineralogy composition only (Li et al., 2021; Scanza et al., 2015), but the measured mass could likely also include non-dust particles, such as sea salt (Kandler et al., 2011; Zhang et al., 2006), sulfate (Kandler et al., 2007), biomass burning aerosols (Ansmann et al., 2011; Johnson et al., 2008), or other air pollution aerosols (Huang et al., 2010; Yuan et al., 2008). This contamination of non-dust aerosols on the measurement is especially important for dust in the fine-mode size where the instrument cannot distinguish dust from the fine-sized non-dust aerosols.

4. Different representative periods. Some of the observations were not made for a period long enough to be taken as a representative of climatology (see Tables 3–5). Also considering point (1), the model–observation comparison may be subject to change because of interannual variability or the episodic character of dust aerosols (Li and Sokolik, 2018b; Mahowald et al., 2011b).

5. Uncertainty in the measurements. In addition to contamination of non-dust aerosols on the measurement of dust, there is also uncertainty due to assumed dust shape and complex refractive index to derive dust size, particularly for particles > 1 µm (Laskin et al., 2006), and error in AERONET AOD retrievals (e.g., the cloud-screening algorithms; Levy et al., 2010) and in the method used to filter out the contribution of non-dust aerosols. Note that difference exists between clear-sky from observations and all-sky AOD/DOD from the model and aerosol models but the difference is not a considerable error source (tested; not shown).

6. The method of selecting AERONET sites may introduce uncertainty because of the possible mismatch between simulated and observed AOD for both dust and non-dust aerosols.

8 Concluding remarks and outlook

This study has compared how different modeling representations of the dust emission schemes, the aerosol dry deposition schemes, transported dust particle size distributions, and the dust shape treatments affect the modeled dust cycle in CESM2.1–CAM6.1. We have evaluated model performance using different combinations of those modifications using offline dynamics by comparing the modeled dust properties (DOD, dust surface concentrations, dust deposition fluxes, atmospheric size distribution of transported dust, and dust DRE efficiency at the top of the atmosphere) that are related to the dust life cycle with (semi-) observations in the current climate. Since the new more physically based dust emission scheme shows substantial improvements on the model–observation comparison and the updated aerosol dry deposition scheme corrects the overestimated fine-mode deposition velocity, future model developments will be focused on introducing both these features into a future official CAM version for the benefit and use of the whole community. Results of this work therefore inform modelers how well these new features will improve model performance in reproducing the dust cycle in CESM.

Our analysis suggests that reverting the geometric standard deviation of the transported dust size distribution (coarse mode) from the default 1.2 to 1.8 imposes the most important change among what we introduced to CAM6.1 to the modeled dust cycle, though the linear assumption between DOD and the other dust quantities based on which we rescaled up the concentrations, deposition, burdens, and DRE of dust in the size distribution simulations introduces uncertainty. Since the defaulted 1.2 is too narrow to simulate the dust lifetime, in the next released model version, we recommend reverting the geometric standard deviation to 1.8, as in CAM5. This reverse may require a split of representation of dust and the stratospheric aerosols in the coarse mode, for which the narrow coarse-mode size distribution works better (Mills et al., 2016), and some changes to sea salt.

With the global DOD similarly comparable in different cases because of the retuning we applied or slight impacts by the updates on DOD, the modifications on dry deposition and emission schemes, and the gravitational settling due to dust asphericity only slightly changed the simulated global dust loadings/burdens and deposition. However, regionally, large difference among different model results for dust loadings/burden/DOD and deposition are found. These stem either from the choice of the dust emission schemes (BRIFT versus DEAD) or the width of the coarse-mode size distribution. Consequently, it is due primarily to the inclusion of the new dust emission scheme and the correction to the transported dust size contribution but not use of the new dry deposition scheme and accounting for dust asphericity that the new model, CAM6.α, shows improvements. It is worth noting that the results obtained in this study rely on the models with the offline dynamics, which is subject to change while using the predicted meteorology field online.

Overall, with the offline dynamics, the new model, CAM6.α:

1. can better capture the climatology and seasonal variation of DOD at more observational sites than the default model, CAM6.1, bearing in mind the uncertainty in the measurement and in the way that we did the model–data comparison;

2. results in a dust DRE that is regionally substantially different from CAM6.1 (e.g., stronger warming over most land areas except over South America and stronger cooling over the North Atlantic Ocean; Fig. 10d). Though the opposite change to dust DRE in one region to another partially cancels, (Fig. S17), its influence on the global mean dust DRE remains large (relative change > 55 %; Fig. 1l).

Still, there exists large uncertainty in modeling the global and regional dust cycle in comparison with observations. This large uncertainty could partially result from the constants used in the parametrizations that affect the dust emission and transport processes, such as the critical LAI threshold, the hygroscopicity of dust, and the prescribed scavenging coefficient, though the default values in the model have been used during the past decade in CAM of different versions. In addition, further development and studies focusing on the following processes and dust properties, which the current model does not represent well or omits entirely, may be helpful for further improving the simulation of the dust cycle in CESM:

1. For the dust emission parameterization, the threshold friction velocity calculated in both BRIFT and DEAD does not account for the spatiotemporal variability of the soil properties (e.g., soil grain size distribution, aggregate state, and static drag partition due to rocks etc.; Leung et al., 2022; mainly limited by the sparse information; Kok et al., 2014b) in addition to the soil moisture. The current dust module in CAM6.1 also does not consider the roughness effect due to the presence of non-erodible elements (i.e., rocks and pebbles) on the threshold velocity calculation (Marticorena and Bergametti, 1995). Also, crusted surface layer present at the erodible surface can greatly reduce the wind erodibility by increasing the particle cohesion, and thus the dust emission rate, compared to the surface that does not consist of consolidated aggregates (Rice and McEwan, 2001; Rodriguez-Caballero et al., 2022).

2. The models used here did not simulate anthropogenic dust emissions due to human activities (i.e., agricultural practices, such as overgrazing, and fugitive dust from roads and construction), which may constitute a considerable fraction of the total dust emissions (Ginoux et al., 2012). This could likely be a reason for the underestimated dust emission in the northern high-latitude regions (Sect. 5.1.1), for instance, at the Moscow metropolitan area (∼ 56^∘ N, ∼ 37^∘ E), one of the most significant northern high-latitude sources generated on paved roads and roadside soils (Kasimov et al., 2020), which the current model does not include.

3. Comparisons with the constrained global dust size distribution and measurements downwind of North Africa suggest that the model underestimates dust aerosols in the coarse mode with the geometric diameter > 5 µm and misses aerosol particles with the geometric diameter > 10 µm (Fig. 6). The former happens likely due to an underestimate of dust aerosol particles in that size range upon emissions and/or the removal rate of those particles being too high during transport in the model (Adebiyi and Kok, 2020; Meng et al., 2022), the reason for which is still under exploration. For the latter, extending the dust size range to include particles with the geometric diameter > 10 µm in CAM6 is a worthy endeavor, such as in CAM5 by Ke et al. (2022).

4. As previously noted (Wu et al., 2018), some of the variables in the dry deposition parameterizations could vary in different seasons for certain land cover and land use types, such as the roughness length, Z₀, in Z01 and the displacement height of the canopy, h, in PZ10, for which a fixed climatological mean is used in the models. How accounting for the seasonal variation of those variables in the model can affect the dust cycle modeling deserves further exploration.

5. Compared to bulk dust, modeling dust aerosol as component minerals could better reproduce the observed spatiotemporal variability of dust optical properties and thus the dust DRE efficiency (Fig. 9) with the offline dynamics in the current climate. But the current atlas of soil mineralogy and the optical properties of key minerals (i.e., iron oxides) contain large uncertainties which should be better quantified in the future, such as that planned for the Earth Surface Mineral Dust Source Investigation (EMIT) (Green et al., 2020; Thompson et al., 2020) and in our ongoing work, respectively.

The comparison of modeling the global and regional dust cycle with observations itself is limited by the spatial and temporal coverage of observations, especially for high-latitude dust, particularly dust in the Southern Hemisphere. More intensive measurements on concentration, deposition, atmospheric loading, shape parameters, size distribution, and optical properties of dust aerosols at varied spatiotemporal scales would also help better represent dust and project climate changes in the global climate models.