The relative importance of wind and hydroclimate drivers in modulating the interannual variability of dust emissions in Earth system models

Windblown dust emissions are controlled by near-surface wind speed and sediment erodibility, the latter modulated by hydroclimate and land-use conditions. Accurate representations of these drivers are critical for reproducing historical dust variability and projecting future dust changes in Earth system models (ESMs). This study examines the discrepancies among 21 ESMs in the relative importance of wind speed versus five hydroclimate drivers in explaining the historical (1980–2014) variability of dust emissions from global drylands. In hyperarid areas, models show poor agreement in the simulated dust variability, with only 9 % out of 210 inter-model comparisons exhibiting significant positive correlations. In contrast, arid and semiarid areas exhibit a dual pattern driven by a “double-edged sword” effect of land surface memory: models with coherent hydroclimate variability show better agreement, whereas those with divergent hydroclimate representations show larger disagreement. While the ESMs capture the dominant role of wind speed in hyperarid areas, they diverge markedly in the relative contributions of wind and hydroclimate drivers in arid and semiarid areas. Replacing the Zender et al. (2003) dust scheme with the Kok et al. (2014) scheme in CESM and E3SM generally strengthens hydroclimate influences while reducing wind speed contributions to simulated dust variability. MERRA-2 reanalysis produces stronger wind influences than most ESMs across all dryland regions. These results underscore the need for improved near-surface wind simulations in hyperarid areas and more realistic land surface and hydroclimate representations in arid and semiarid areas to reduce uncertainties in global dust emission simulations.

1 Introduction

Windblown dust emission is modulated by near-surface wind speed and the supply and erodibility of fine-grain sediments, which collectively determine the timing, location, duration, and intensity of dust events. The most abundant sediment supply is typically found in low-relief regions with thick layers of unconsolidated materials produced by weathering, fluvial, and/or aeolian processes (Bullard and Livingstone, 2002; Bullard et al., 2011). The erodibility of these fine materials depends on environmental conditions such as surface soil moisture and armoring (e.g., vegetation, soil crusts) which determine the minimum or threshold wind velocity that must be reached to initiate the saltation-sandblasting process (e.g., Zender and Kwon, 2005; Shao et al., 2011; Xi and Sokolik, 2015b). The environmental controls of dust emission have been incorporated into Earth system models (ESMs) via parameterizations of dust emission fluxes as a function of various atmospheric, land surface, and soil parameters, many of which are interactively determined within the models. For example, the horizontal saltation dust flux is calculated as a function of the third or fourth power of wind velocity, reflecting the dominant role of infrequent, high-wind events in dust production (e.g., Owen, 1964; Bagnold, 1974; White, 1979; Kok et al., 2012). Early parameterizations in ESMs use static dust source functions to represent the spatially varying sediment supply, with large values associated with low-relief regions where elevated dust burden has been frequently detected by satellite observations (e.g., Ginoux et al., 2001; Prospero et al., 2002; Zender et al., 2003). The sediment supply is typically assumed to be unlimited in ESMs without accounting for depletion or replenishment over time. In addition, the sediment erodibility is closely connected with hydroclimate and land surface processes in ESMs. For example, surface soil moisture, which is simulated by the land model component, is often used to account for increases in erosion thresholds due to enhanced soil particle cohesion under wet conditions (Fécan et al., 1999). Many ESMs also use the bare soil fraction to adjust dust emissions from areas partially covered by snow, ice, or vegetation. Vegetation also acts to increase the aerodynamic surface roughness and reduce the available wind shear stress exerting on erodible surfaces, which can be represented by drag partitioning schemes but is currently not considered in most ESMs (Raupach et al., 1993; Marticorena and Bergametti, 1995; Shao, 2001).

Numerous studies have evaluated the consistency and performance of current ESMs in simulating the global dust cycle under the Aerosol Comparisons between Observations and Models (AeroCom) initiative and Coupled Model Intercomparison Project (CMIP) (Textor et al., 2006; Huneeus et al., 2011; Kim et al., 2014; Wu et al., 2020; Gliß et al., 2021; Zhao et al., 2022; Kim et al., 2024). Overall, these studies suggest that modern-day dust aerosol column burden is reasonably constrained by ground- and satellite-based aerosol optical depth (AOD) retrievals, leading to better inter-model agreement than those in dust emission and deposition estimates. Knippertz and Todd (2012) pointed out that model tuning to match satellite observations, e.g., via the use of dust source functions, induces a compensational effect between dust emission and deposition, both of which lack direct observational constraints at the global scale. Indeed, previous AeroCom and CMIP model intercomparisons reported substantial discrepancies in global dust emission estimates, with differences spanning an order of magnitude, as well as persistent difficulties in reproducing historical dust variability and its relationships with key driving factors (Huneeus et al., 2011; Evan et al., 2014; Evan, 2018; Pu and Ginoux, 2018; Wu et al., 2020; Gliß et al., 2021; Zhao et al., 2022). Kok et al. (2023) further suggested that current ESMs failed to capture the increase of global dust burden since preindustrial times, likely due to inaccurate model representations of the climate and land-use drivers of dust emissions, and/or the dust sensitivity to these driving factors.

The model discrepancies can be partly explained by the choice of dust emission schemes. Earlier schemes rely on static dust source functions to shift emissions towards satellite-observed hot spot regions (e.g., Ginoux et al., 2001; Zender et al., 2003), whereas newer parameterizations replace prescribed dust source functions with more explicit formulations of sediment erodibility that increase the dust emission sensitivity to soil-moisture-dependent erosion thresholds (e.g., Kok et al., 2014b; Leung et al., 2023). Dust schemes also differ in how they represent the sandblasting efficiency, defined as the ratio of the vertical dust flux to the saltation flux: some schemes assume a global constant (e.g., Ginoux et al., 2001; Volodin and Kostrykin, 2016), whereas more sophisticated schemes account for the dependence of sandblasting efficiency on soil properties and wind speed (e.g., Zender et al., 2003; Kok et al., 2014b). The choice of wind speed also varies: while some dust schemes use 10-m winds, others use friction velocity which more accurately represents the wind shear stress acting on erodible surfaces but requires specification of surface roughness length. In general, more sophisticated schemes, which are derived based on small-scale wind tunnel experiments, require more extensive input parameters which are often poorly constrained at climate model grid levels, necessitating assumptions and empirical tunings. For instance, due to limited data availability, surface aerodynamic roughness is often prescribed as a global constant or based on static satellite-derived maps (e.g., Peng et al., 2012; Tegen et al., 2019). ESM-simulated soil water content may lack the accuracy or dynamic range required by dust emission parameterizations. As a result, some models apply additional tunings or alternative treatments of soil moisture effects (e.g., Zender et al., 2003; Volodin and Kostrykin, 2016), while others disable the soil moisture dependence entirely (e.g., Noije et al., 2021; Shevliakova et al., 2024).

Even when using the same dust emission scheme, ESMs can still diverge in dust emission simulations due to differences in model configurations (e.g., horizontal resolution, vertical levels), input datasets, parameter tunings, and coupled physical processes. For instance, the bare soil fraction is determined from land cover type, vegetation fraction, and snow/ice area extent, all of which may differ across ESMs. Vegetation cover itself may be prescribed from satellite climatology or simulated interactively. ESMs also differ in their representations of soil properties (e.g., hydraulic conductivity), soil column structure (e.g., number and thickness of layers), and land surface hydrologic formulations (e.g., precipitation, runoff, evaporation), which collectively determine the surface soil moisture needed by dust emission schemes. The soil moisture effects may be treated inconsistently, e.g., in how models define the residue moisture level below which dust emission is assumed to be independent on soil water content (e.g., Ginoux et al., 2001; Evans et al., 2016; Volodin and Kostrykin, 2016). Moreover, ESMs employ different planetary boundary layer and surface flux parameterizations, which influence the simulation of near-surface winds and extreme wind events pertinent to dust mobilization and transport. Given the inherent differences in dust process representations and the lack of direct observational constraints, it is thus not surprising that ESM-simulated dust emission fluxes exhibit substantial discrepancies, as documented in previous AeroCom and CMIP intercomparison studies. Therefore, ESM-simulated dust emission fluxes are best viewed as an unconstrained, model-specific quantity characterized by a dynamic range defined by the parameterizations, configurations, and parameter tunings of individual models, similar to Koster et al. (2009)'s view on root-zone soil moisture.

While model discrepancies in global dust emission estimates are well documented, a key remaining question is how consistently and accurately current ESMs represent the temporal variability of dust emission and its sensitivity to underlying physical drivers. Addressing this question is essential for understanding and reducing model uncertainties in predicting dust responses to climate and land-use changes. In this study, we focus on the interannual variability of ESM-simulated dust emissions from global drylands and apply a statistical framework to diagnose the physical controls of dust emissions. Based on fully-coupled historical simulations from a suite of ESMs and two aerosol reanalysis products, we examine the extent of inter-model agreement in simulating the interannual dust variability across different climate aridity regimes and quantify the relative importance of a common set of physical drivers in explaining the simulated dust variability within individual models. Although satellite-derived dust AOD and in-situ dust measurements provide valuable constraints on dust variability (e.g., Prospero and Lamb, 2003; Voss and Evan, 2020), they integrate the effects of emission, transport, and deposition, making it difficult to isolate the emission process itself. Also, due to lack of global validation data, we focus on diagnosing inter-model inconsistency in representing the dust emission variability and its physical controls, rather than validating individual model performance against observations.

The remainder of this paper is organized as follows. Section 2 describes the ESMs and aerosol reanalysis datasets and the dominance analysis technique. Section 3 presents results on comparing model-simulated dust emission fluxes and the relative influences of wind speed versus hydroclimate drivers. Section 4 summarizes the main findings of this study.

2 Data and Approach

2.1 ESMs and aerosol reanalysis

We consider a total of 21 ESMs as summarized in Table 1, including 18 models from the CMIP6 fully coupled historical experiment (1980–2014). For each model, we use the first ensemble member (r1i1p1f1) unless otherwise stated. Two CESM variants use the same dust emission scheme of Zender et al. (2003) (hereafter the Zender scheme) but different atmospheric components: the Community Atmosphere Model (CESM2-CAM-Zender) versus the Whole Atmosphere Community Climate Model (CESM2-WACCM-Zender). In addition, we conducted a CESM experiment (2004–2013) using the dust emission scheme of Kok et al. (2014b) (hereafter the Kok scheme; CESM2-CAM-Kok) (Li et al., 2022). We further conduct two E3SM model experiments (1980–2014) coupled with the Zender (E3SM2-Zender) and Kok (E3SM3-Kok) schemes, respectively (Feng et al., 2022; Xie et al., 2025). A key difference between the two schemes is that, the Zender scheme relies on a prescribed, time-invariant dust source function that shifts emissions towards contemporary dust source regions, whereas the Kok scheme applies more physically based parameterizations of soil erosion thresholds, thereby improving dust simulations without using prescribed dust source functions (Kok et al., 2014a). The paired CESM and E3SM experiments allow us to examine how the choice of dust emission schemes or host models affects the simulated dust variability and sensitivity to driving factors. Nonetheless, it is important to note that comparing these experiments is complicated by additional model differences. For example, CESM2-CAM-Zender does not account for dust mineralogy, whereas CESM2-CAM-Kok simulates dust as mineral components with observationally constrained mineral optical properties (Li et al., 2024). This may lead to inconsistent radiative feedback on meteorology and dust emissions. Similarly, E3SM3 incorporates extensive model updates relative to E3SM2, which may affect near-surface meteorological and land surface conditions relevant to dust emissions (Xie et al., 2025).

Zender et al. (2003)Gettelman et al. (2019)Zender et al. (2003)Albani et al. (2015)Kok et al. (2014b)Li et al. (2022)Zender et al. (2003)Feng et al. (2022)Kok et al. (2014b)Xie et al. (2025)Peng et al. (2012)Swart et al. (2019)Peng et al. (2012)Sigmond et al. (2023)Tegen et al. (2002)Séférian et al. (2019)Tegen et al. (2002)Noije et al. (2021)Miller et al. (2006)Miller et al. (2021)Miller et al. (2006)Miller et al. (2021)Miller et al. (2006)Rind et al. (2020)Woodward (2011)Roberts et al. (2019)Woodward (2001)Woodward et al. (2022)Volodin and Kostrykin (2016)Volodin (2022)Balkanski et al. (2004)Lurton et al. (2020)Shao et al. (1996)Yukimoto et al. (2019)Takemura et al. (2009)Tatebe et al. (2019)Takemura et al. (2009)Hajima et al. (2020)Tegen et al. (2019)Mauritsen et al. (2019)Zender et al. (2003)Seland et al. (2020)Ginoux et al. (2001)Randles et al. (2017)Shao et al. (1996)Yumimoto et al. (2017)

Table 1 Summary of Earth system models and aerosol reanalysis datasets considered in this study. Dust source function (DSF) indicates whether a static dust source function is used. Leaf area index (LAI) indicates whether LAI is treated as a prognostic variable. D_m is the dust particle diameter upper limit. u_* is friction velocity. u₁₀ is 10 m wind velocity.

Several other model families share common heritage but differ in physics options and configurations. For example, CanESM5.1 incorporated physics and technical changes that improved mean climate and dust simulations relative to CanESM5.0 (Sigmond et al., 2023). Three GISS-E2 models use the same dust scheme from Miller et al. (2006) but differ in model version (2.1 vs. 2.2) and aerosol microphysics scheme: One-Moment Aerosol (OMA; ensemble member r1i1p3f1) versus Multiconfiguration Aerosol TRacker of mIXing state (MATRIX; ensemble member r1i1p5f1) (Rind et al., 2020; Miller et al., 2021). UKESM1.0 is developed based on the HadGEM3-GC3.1 general circulation model. They use the same dust scheme from Woodward (2001) but employ different parameter tunings and dust source representations (Woodward et al., 2022). MIROC-ES2L builds upon the MIROC general circulation model version 5.2 (MIROC5) (Hajima et al., 2020), while MIROC6 incorporates physics updates that improved the mean climate state and internal variability compared to MIROC5 (Tatebe et al., 2019). Both MIROC-ES2L and MIROC6 use the dust scheme from the SPRINTARS model (Takemura et al., 2009).

The ESMs are further compared with two aerosol reanalysis products with dust emission flux output: Modern-Era Retrospective Analysis for Research and Applications Version 2 (MERRA-2, 1980–2014) (Gelaro et al., 2017), and Japanese Reanalysis for Aerosol (JRAero, 2011–2017) (Yumimoto et al., 2017). Dust emission in MERRA-2 is simulated using the Ginoux et al. (2001) parameterization within the GOCART aerosol module of GEOS-5 model. In JRAero, dust emission is simulated using the Shao et al. (1996) energy-based scheme (same as in MRI-ESM2.0) within the Japan Meteorological Agency MASINGAR mk-2 global aerosol transport model (Yumimoto et al., 2017; Yukimoto et al., 2019). In both MERRA-2 and JRAero, the meteorological inputs for dust emission calculations are constrained via data assimilation of in situ and remote sensing observations (including surface and upper-air wind measurements), which improves near-surface wind simulations compared to CMIP6 models (Gelaro et al., 2017; Yumimoto et al., 2017). The surface soil moisture in MERRA-2 also benefits from assimilation of observation-corrected precipitation. While both reanalyses assimilate bias-corrected total AOD from satellites, the AOD assimilation is expected to have limited influence on dust emission simulations.

To facilitate the comparison among the ESMs and reanalysis products, we classify global drylands into three climate zones–hyperarid, arid, and semiarid–based on the aridity index (AI), defined as the ratio of climatological mean precipitation to potential evapotranspiration for 1970–2000 following (Zomer et al., 2022). Dry subhumid areas are not considered due to their negligible dust contributions. The hyperarid climate zone is defined as AI ≤ 0.05, arid zone as 0.05 < AI ≤ 0.2, and semiarid zone as 0.2 < AI ≤ 0.5. As shown in Fig. 1, hyperarid areas primarily cover North Africa, Arabian Peninsula, Iranian Plateau, and Tarim Basin. Arid and semiarid areas cover other major dust sources, including the Sahel (North Africa), Turan Depression (Central Asia), Gobi Desert (East Asia), Thar Desert (South Asia), Kalahari Desert (Southern Africa), Chihuahua Desert (North America), Patagonia steppe (South America), and the Great Sandy and Simpson Deserts (Australia). Generally, hyperarid areas are dominated by permanently dry, barren surfaces with very low hydroclimate variability, such that dust emission is expected to be primarily controlled by wind speed. In contrast, arid and semiarid areas experience greater precipitation and hydroclimate fluctuations, which are expected to exert stronger influence on sediment erodibility and dust emissions.

Figure 1 Definitions of climate aridity zones for model intercomparisons.

2.2 Dominance analysis technique

Previous studies commonly used linear regression coefficients to quantify the dust sensitivity to its physical drivers (e.g., Pu and Ginoux, 2016; Aryal and Evans, 2021; Zhao et al., 2022). In multiple linear regression, a regression coefficient represents the mean change in the response variable (e.g., dust emission flux or AOD) associated with a unit change in a given predictor, while holding all other predictors constant. This interpretation assumes mutual independence among predictors, an assumption that is often violated by strong correlations among dust emission drivers. Consequently, linear regression coefficients can cause misleading inference of the relative importance of predictors. In addition, regression coefficients, standardized or not, may not provide a consistent basis for comparing predictor importance across different ESMs, due to their inconsistent dynamic ranges.

In this study, we use the dominance analysis technique to quantify the relative importance of multiple correlated predictors in explaining the variability of monthly total dust emission fluxes simulated by individual ESMs and the MERRA-2 reanalysis. JRAero is excluded from this analysis due to missing predictor data. Although the ESMs differ in how they partition total emissions into discrete size modes or bins, the size partitioning has minor effects on diagnosing the emission process itself. The predictors considered here operate upstream of the size partitioning, and control the initiation and magnitude of total dust emissions rather than its size-resolved characteristics. Here we consider a common set of six predictors from each model: 10 m wind speed, total precipitation (including liquid and solid phases), water content in the uppermost soil layer, 2 m specific humidity, 2 m air temperature, and leaf area index (LAI). Among them, 10 m wind speed represents the wind shear stress driving dust mobilization, while the remaining variables collectively represent the hydroclimate controls on sediment erodibility. The selected predictors are either directly used in dust emission parameterizations or strongly correlated with dust emission intensity, as shown in numerous studies (e.g., Engelstaedter et al., 2003; Ravi et al., 2006; Zou and Zhai, 2004; Cowie et al., 2015; Kim and Choi, 2015; Xi and Sokolik, 2015b, a). Note that we donot not include all the physical drivers represented in each model because of limited data availability in the CMIP6 online archive, and because some models incorporate additional drivers not used by others. Hence we focus on a common set of six readily available predictors to provide a consistent and fair comparison across the ESMs and MERRA-2 reanalysis.

Dominance analysis is applied to the ESMs and MERRA-2 over grid cells with nonzero dust emissions. The dust emission fluxes and predictors are first deseasonalized by subtracting month-wise climatological means and then normalized to 0–1 range via min-max scaling. Dominance analysis quantifies the marginal contribution of each predictor to the total explained variance (R²) by evaluating all possible subset models (2^p−1 subsets for p predictors) in a multiple linear regression framework (Budescu, 1993; Azen and Budescu, 2003). The approach first calculates the average incremental contribution of each predictor to the total R² across all subset models of the same size (i.e., models with the same number of predictors). These incremental R² values are then averaged to obtain the predictor’s overall contribution to the total R². A key feature of this approach is that the sum of individual predictor contributions equals the total R² of the full model (i.e., with all predictors included), thereby allowing the partitioning of total explained variance among correlated predictors. The resulting grid-level predictor R² values from the ESMs and MERRA-2 are used to assess the spatial variability of predictor influences over different climate zones within each model, as well as the inter-model consistency in representing the relative importance of wind speed versus five hydroclimate drivers.

3 Results

3.1 Climatological distribution

Figure 2 displays the climatological mean annual dust emission fluxes from the 21 ESMs, their ensemble mean, and the MERRA-2 and JRAero reanalyses during 2005–2014 (2004–2013 for CESM2-CAM-Kok and 2011–2017 for JRAero). All datasets capture the global dust belt stretching from West Africa across the Middle East to East Asia, as well as weaker sources in the Americas and Australia. Among the ESMs, E3SM3-Kok and HadGEM2-GC31 simulate the most extensive dust-emitting areas extending to high-latitude and subhumid areas. CESM2-CAM-Zender, CESM2-WACCM-Zender and NorESM2 restrict emissions to regions where the dust source function exceeds 0.1, resulting in discrete and spatially limited emission patterns. Conversely, E3SM2-Zender employs the original dust source function of Zender et al. (2003), producing a more spatially continuous emission pattern (Fig. 2e).

Figure 2 Climatological mean dust emission fluxes from (a–u) 21 Earth system models , (v) model ensemble mean, (w) MERRA-2 reanalysis, and (x) JRAero reanalysis. Global annual total emissions are displayed on each panel.

Global annual dust emissions simulated by the ESMs vary greatly, ranging from 890 to 7727 Tg yr⁻¹ with nearly an order of magnitude difference (Fig. 2a–u). The ensemble mean estimate (Fig. 2v) is 2857 Tg yr⁻¹ with a standard deviation of 1835 Tg yr⁻¹, corresponding to a diversity of 64 % (defined as the ratio of standard deviation to ensemble mean). Based on models with an upper particle size limit of 20 µm, global dust emissions vary from 1061 to 6561 Tg yr⁻¹, with a mean of 3048 Tg yr⁻¹ and diversity of 55 %. The ensemble mean is close to JRAero (2780 Tg yr⁻¹, Fig. 2x), but considerably higher than MERRA-2 (1605 Tg yr⁻¹, Fig. 2w). Also, the ensemble mean exhibits a more spatially homogeneous pattern over North Africa and the Arabian Peninsula, whereas MERRA-2 and JRAero display more heterogeneous and localized emission patterns.

The model discrepancies in dust emission magnitude are consistent with previous assessments. For example, Huneeus et al. (2011) compared 14 AeroCom Phase I models and reported a global dust emission range of 500–4400 Tg yr⁻¹ (diversity = 58 %), of which seven using a 20 µm upper size limit yielded 980–4300 Tg yr⁻¹ (diversity = 46 %). Similarly, Gliß et al. (2021) compared 14 AeroCom Phase III models and reported a range of 850–5650 Tg yr⁻¹ with a diversity of 64 %. Based on 15 CMIP5 models, Wu et al. (2020) reported a range of 740–8200 Tg yr⁻¹ (diversity = 66 %), with seven models using particle diameters up to 20 µm producing 740–3600 Tg yr⁻¹ (diversity = 43 %). More recently, Zhao et al. (2022) examined 15 CMIP6 AMIP models and reported a range of 1400–7600 Tg yr⁻¹ with a diversity of 61 %. Collectively, these studies, along with our results, demonstrate persistent large model uncertainties in global dust emission estimates despite advances in model resolutions and physics parameterizations, which reflects the unconstrained, model-specific nature of dust emission fluxes.

Figure 3 displays the contributions of different climate zones to global dust emissions. Based on the model ensemble mean, global dust emissions are partitioned as 61 % from hyperarid, 27 % from arid, and 5 % from semiarid zones. In comparison, MERRA-2 and JRAero allocate the majority of dust emissions to hyperarid and arid zones, with negligible contributions from the semiarid zone.

Figure 3 Contributions of different climate zones to global annual dust emissions. Numbers indicate percentages above 5 %.

The hyperarid zone accounts for more than half of global emissions in all ESMs except CanESM5.0, CanESM5.1, and INM-CM5.0. These models simulate relatively uniform emission patterns with less than 50 % from hyperarid areas, possibly related to their dust emission parameterizations. As noted in Sigmond et al. (2023), parameter tunings related to the hybridization of dust tracers caused spurious dust events in CanESM5.0. An interpolation error in the bare soil fraction also contributed to dust simulation biases in CanESM5.0 and CanESM5.1 compared to other CMIP6 models and satellite aerosol observations (Sigmond et al., 2023). The newer CanESM5.1 simulates 20 % less dust globally but similar spatial distributions compared to CanESM5.0. In INM-CM5.0, the vertical dust flux is calculated as a function of wind speed alone, without accounting for land surface effects on the threshold wind velocity (Volodin and Kostrykin, 2016; Volodin, 2022). While this simplification may be appropriate for hyperarid areas, it can overestimate emissions over arid and semiarid areas where increased soil wetness and vegetation cover suppress dust mobilization.

Contributions of the arid zone range from 8 % (CESM2-CAM-Kok) to 37 % (UKESM-1.0), indicating substantial model discrepancies compared to the hyperarid zone. The discrepancies become even larger over the semiarid zone, where the emission fraction ranges from less than 1 % to 18 %. Particularly, four models allocate more than 10 % to the semiarid zone: CanESM5.0 (16 %), CanESM5.1 (18 %), INM-CM5.0 (15 %), and UKESM1.0 (12 %). Overall, as the climate regime transitions from hyperarid to semiarid, model-estimated dust source strengths become less consistent, revealing increasing uncertainty in how ESMs represent dust sensitivity to hydroclimate conditions.

Among the ESMs, CESM2-CAM-Zender and CESM2-WACCM-Zender produce nearly identical total emissions and spatial patterns, suggesting that the choice between CAM and WACCM atmospheric components has minimal effect. The paired CESM and E3SM experiments, however, show opposite tendencies: contributions of the hyperarid zone increase from 61 % in CESM2-CAM-Zender to 88 % in CESM2-CAM-Kok, but slightly decrease from 63 % in E3SM2-Zender to 58 % in E3SM3-Kok. The GISS-E2 models produce consistent distributions across different climate zones, although total emissions are about 40 % lower when using the MATRIX aerosol scheme, possibly due to parameter tunings or underrepresentation of coarse dust (diameter > 5 µm) in the MATRIX modal size distribution, as noted in Bauer et al. (2022).

UKESM1.0 emits nearly twice as much dust as HadGEM3-GC3.1, and exhibits slightly more uniform spatial distributions. As described in Woodward et al. (2022), UKESM1.0 is built upon HadGEM3-GC3.1 but applies parameter tunings that enhance friction velocity and suppress soil moisture, effectively increasing the wind gustiness and soil aridity, leading to more emissions in UKESM1.0. The three Japanese models (MRI-ESM2.0, MIROC-ES2L, and MIROC6) also differ markedly in total emissions and, to a lesser extent, spatial distributions. MRI-ESM2.0 produces similar regional fractions to JRAero but nearly doubles the total amount. Despite using the same dust parameterization, MIROC-ES2L emits roughly five times more dust than MIROC6. This discrepancy can be largely explained by stronger winds in MIROC-ES2L, which produces 50 % higher global mean wind speed than MIROC6. Moreover, MIROC6 prescribes non-zero LAI even in hyperarid areas, likely further suppressing dust generation relative to MIROC-ES2L (Hiroaki Tatebe, personal communications, 2025).

3.2 Interannual variability

This section examines the consistency of ESMs in simulating the interannual variability of dust emissions. Monthly dust emission fluxes from all ESMs are first regridded to a common resolution of 0.9° × 1.25° (the native grid of CESM2). To remove the influence of annual cycles, month-wise climatological means are subtracted from each grid cell to calculate deseasonalized dust emission flux anomalies. Spearman’s rank correlation coefficients are then calculated between the monthly anomalies for every possible model pair. With 21 ESMs, this results in 210 pairwise comparisons. To quantify the overall model agreement, we calculate the percentage of model pairs exhibiting statistically significant (p≤0.1), positive correlations. A higher percentage indicates stronger model agreement in simulating the dust variability, and vice versa. The results are displayed in Fig. 4.

Figure 4 Percentage of statistically significant (p≤0.1), positive correlations out of 210 pairwise comparisons of deseasonalized monthly dust emission fluxes from 21 Earth system models. Black contours indicate the model ensemble mean emission flux of 10 and 100 Tg yr⁻¹.

Despite its dominant contribution to global dust emissions, the hyperarid zone exhibits poor model agreement, with generally less than 10 % of pairwise comparisons showing statistically significant positive correlations. Because dust emissions from hyperarid areas are predominantly controlled by near-surface wind speed, this poor agreement reflects inconsistent wind simulations among the ESMs. Indeed, we find that only 10 % of model pairs produce positively correlated monthly mean wind speed anomalies. Similarly, Evan (2018) reported that dust-producing winds over the Sahara Desert are driven by large-scale meteorological processes and that most CMIP5 models failed to capture the near-surface wind variability. These results suggest that improving near-surface wind simulations can potentially reduce discrepancies in simulating the dust emission variability over hyperarid regions.

In contrast to the hyperarid zone, arid and semiarid zones – such as the Sahel, South Asia, East and Central Asia, Australia, and North America – exhibit significantly better agreement. To further assess how model consistency varies with aridity, Fig. 5 presents pairwise correlation matrices based on dust emission flux anomalies averaged over hyperarid, arid, and semiarid zones. The percentage of statistically significant, positively correlated model pairs increases from 9 % in the hyperarid zone to 13 % in the arid zone and 16 % in the semiarid zone. Meanwhile, the semiarid zone shows a larger percentage of negatively correlated model pairs (16 %) than the hyperarid (6 %) and arid (7 %) zones. This dual pattern suggests that as the climate regime transitions from hyperarid to semiarid, the ESMs exhibit both stronger agreement and heightened disagreement in simulating the interannual variability of dust emissions.

Figure 5 Spearman's rank correlation coefficients between dust emission flux anomalies averaged over hyperarid, arid, and semiarid climate zones. Dots indicate statistically significant correlations (p≤0.1). Summary tables are based on Earth system models only (MERRA-2 and JRAero not included).

This behavior can be explained by the influence of antecedent land surface conditions on sediment erodibility in semiarid environment such as temperate grasslands and steppes (Shinoda and Nandintsetseg, 2011; Nandintsetseg and Shinoda, 2015). In these regions, factors such as precipitation, soil moisture, and vegetation growth–decay cycles have lagged and long-lasting effects on sediment erodibility. For example, dry anomalies during the wet season, such as reduced rainfall or earlier snowmelt, can reduce soil inter-particle cohesion and suppress vegetation growth, thereby prolonging bare soil exposure and increasing the wind erosion risk. This delayed response exemplifies the land surface memory effect, in which the slow adjustment of soil and vegetation conditions over weeks to months influences subsequent dust emission potentials long after the initial hydroclimate forcing (e.g., drought). Therefore, we hypothesize that the simultaneous increase in both model agreement and disagreement from hyperarid to semiarid zones reflects a “double-edged sword” effect of land surface memory: models with coherent representations of hydroclimate variability tend to converge in the simulated dust variability (i.e., more positive correlations), whereas those with divergent hydroclimate representations diverge in the dust variability (i.e., more negative correlations).

To verify this hypothesis, we examine the statistical association between pairwise model correlations in dust emissions and those in hydroclimate variability. Specifically, we perform a principal component analysis (PCA) of the five hydroclimate variables (i.e., precipitation, soil moisture, specific humidity, air temperature, LAI) separately for the hyperarid, arid, and semiarid zones. The leading principal component (PC1), which explains at least 40 % of the total variance in all zones, is used as a proxy for the dominant hydroclimate variability. Spearman’s rank correlation coefficients are then computed for all pairwise comparisons of deseasonalized monthly PC1 values, following the same approach as in Fig. 5.

Figure 6 compares the correlation coefficients for model pairs with the same sign (i.e., both positive or both negative) in dust emission fluxes and hydroclimate PC1. The regression slope and coefficient of determination (r²) quantify the degree of statistical association between inter-model correlations in dust emission and hydroclimate variability. The positive relationships across all climate zones suggest that ESMs with stronger consensus in hydroclimate variability also tend to produce more consistent dust variability. More importantly, the number of significantly correlated model pairs (N) and correlation strength (slope and r²) show significant increases from hyperarid to semiarid zones. This finding supports our hypothesis on the dual role of land surface memory: it tends to improve agreement among ESMs with coherent hydroclimate representations, while simultaneously amplifying disagreement among those with divergent hydroclimate variability.

Figure 6 Statistical associations between pairwise model correlation coefficients (p≤0.1 shown in red) in dust emission fluxes and hydroclimate variability over (a) hyperarid, (b) arid, and (c) semiarid climate zones.

3.3 Relative importance of wind speed versus hydroclimate drivers

In this section, we present dominance analysis results on the joint and relative influences of wind speed and hydroclimate drivers on the simulated dust variability within individual ESMs. Figure 7 shows the total R² by the six predictors combined. Results for CESM2-WACCM-Zender and NorESM2 are very similar to those for CESM2-CAM-Zender and thus not shown. The ESMs show large discrepancies in the total R², reflecting inherent differences in the coupling strength between dust emission and the selected predictors. When ranked by the global mean R², CanESM5.0 and CanESM5.1 show the lowest explanatory power of the selected predictors, followed by MPI-ESM1.2, MIROC6, and EC-Earth3-AerChem. The low total R² can be explained by several factors. We only consider six common predictors and may omit other predictors that are specific to some models (such as stem area index). Model biases (see Sect. 3.1) may weaken or distort the relationship between dust emission and its physical drivers. Using simplified parameterizations and/or static land surface input can reduce the dust sensitivity to hydroclimate conditions. In addition, because dust emission is governed by highly nonlinear threshold processes, its dependence on the predictors may deviate from the linear assumption underlying dominance analysis. As shown in Fig. 7, total R² values are generally lower in arid and semiarid areas than in hyperarid areas, likely due to increased nonlinearity between dust emission and hydroclimate variables that diminishes the explanatory power of multilinear regression.

Figure 7 Total explained variance (R²) in dust emission fluxes by six near-surface predictors (wind speed, precipitation, soil moisture, specific humidity, air temperature, and LAI) in Earth system models and MERRA-2. Global mean R² values are shown on each panel.

Despite these limitations, most ESMs produce significant total R² values, especially in hyperarid areas with R² generally above 0.5. Replacing the Zender with the Kok dust scheme generally reduces the total R² in both CESM and E3SM (Fig. 7a–d). GISS-E2 models show little differences between the OMA and MATRIX aerosol schemes, and a modest increase from version 2.1 to 2.2. UKESM1.0 and HadGEM3-GC3.1 show minimal differences, both showing high R² values globally. MIROC6 yields lower R² than MIROC-ES2L, particularly in hyperarid areas. MERRA-2 produces higher R² than most ESMs, especially in arid and semiarid zones, indicating stronger coupling between dust emission and the selected predictors.

Figure 8 presents the ratio of the wind-associated R² to the combined R² of five hydroclimate variables (precipitation, soil moisture, specific humidity, air temperature, and LAI). In all ESMs, the wind-to-hydroclimate R² ratio is well above 1 over hyperarid areas, which is consistent with the parameterization of vertical dust flux as a power-law function of wind speed in all models, and the expectation that dust emission from permanently dry and sparsely vegetated surfaces is primarily controlled by wind speed. In contrast, arid and semiarid areas exhibit much larger inconsistency, with ratios either above or below 1 depending on the model. This behavior may be due to inconsistent representations of hydroclimate controls on sediment erodibility, which in turn leads to disagreement in the relative importance of wind versus hydroclimate drivers. For example, most ESMs incorporate soil moisture as a correction to the erosion threshold velocity albeit using different formulations, while INM-CM5.0 treats soil moisture as a simple threshold above which dust emission is switched off. Such inconsistencies inevitably produce varying coupling strengths between dust emission and hydroclimate variables across models.

Figure 8 The ratio of wind-associated R² to the combined R² of five hydroclimate variables (precipitation, soil moisture, specific humidity, air temperature, and LAI) in Earth system models and MERRA-2.

Based on the wind-to-hydroclimate R² ratios, we classify global dust-emitting areas into three regimes: wind-dominated (ratio > 1.2), hydroclimate-dominated (ratio < 0.8), and equally-important (0.8–1.2). Then we calculate the fractions of dust emissions from these regimes within each model. The results are displayed in Fig. 9. The ESMs show general agreement in the “equally-important” regime, with most models simulating less than 10 % from regions where wind and hydroclimate drivers have nearly equal contributions.

Figure 9 Fractional contributions of wind-dominated, equally-important, and hydroclimate-dominated regimes to global dust emissions in Earth system models and MERRA-2.

The wind-dominated regime accounts for more than 80 % dust emissions in most ESMs, consistent with the dominant contribution of the hyperarid zone (Fig. 3). Two models yield significantly lower estimates: INM-CM5.0 (43 %) and CanESM5.0 (65 %). As shown in Fig. 3, these models simulate relatively homogeneous emission patterns, which consequently diminishes the relative contributions from wind-dominated regions compared with other models. For CESM and E3SM, replacing the Zender with Kok dust scheme slightly reduces the wind-dominated dust fraction: from 85 % to 79 % in CESM, and from 99 % to 96 % in E3SM. The three GISS-E2 models yield similar results, with 87 %–90 % dust from the wind-dominated regime. Likewise, UKESM1.0 and HadGEM3-GC3.1 yield nearly identical estimates, with 90 % of dust from wind-dominated regions. MERRA-2 simulates 98 % emissions from wind-dominated areas, higher than most ESMs. Three models produce significantly higher contributions from semiarid areas than others: CanESM5.0 (29 %), CanESM5.1 (21 %), and INM-CM5.0 (49 %).

To further assess the relative importance of wind speed versus hydroclimate drivers, we compute the fractional contributions of wind speed to the total R² at each model grid cell. The statistical distributions of grid-level wind-associated R² fractions are displayed in Fig. 10. In Fig. 10, if the median wind-associated R² fraction exceeds 50 %, it indicates that wind speed dominates dust variability at more than half of the grid cells within a climate zone. Conversely, if the median falls below 50 %, hydroclimate drivers exert dominant control over the majority of grid cells.

Figure 10 Ridgeline plots for the fractional contributions of wind speed to the total R² over (a) hyperarid, (b) arid, and (c) semiarid climate zones. Black vertical lines indicate median values. Color shading represents mean total R² values.

In the hyperarid zone (Fig. 10a), the majority of ESMs capture the dominant wind control, with median wind-associated R² fractions exceeding 80 %. The three GISS-E2 models display slightly lower wind contributions (67 %–74 %). In contrast, CESM2-CAM-Kok exhibits greater spatial variability and lower wind influence with a median wind-associated R² fraction of 63 %, consistent with the model's elevated hydroclimate influence over West Africa and the Tarim Basin as shown in Fig. 8b. Compared to CESM2-CAM-Kok, CESM2-CAM-Zender captures the expected wind dominance with a median of 87 %. The enhanced hydroclimate influence in CESM2-CAM-Kok relative to CESM2-CAM-Zender persists even when comparing common dust-emitting areas in these two models.

In the arid zone (Fig. 10b), total R² values are generally lower, again reflecting reduced explanatory power of the selected predictors. The ESMs exhibit greater discrepancies in the relative importance of wind speed versus hydroclimate drivers in the arid zone. Specifically, wind speed remains the dominant driver of dust variability in most ESMs, despite increased spatial variability. The three GISS-E2 models simulate nearly equal wind and hydroclimate influences. INM-CM5.0, MIROC-ES2L and MIROC6 produce median wind-associated R² fractions well below 50 %, signifying a transition from wind- to hydroclimate-dominated regimes. CESM2-CAM-Kok also displays this transition, although to a smaller extent with a median of 44 % and large spatial variability. In both CESM and E3SM, replacing the Zender with Kok dust scheme weakens the wind influence and strengthens the hydroclimate influence, with the median wind-associated R² fraction declining from 56 % to 44 % in CESM and from 86 % to 74 % in E3SM. This is somewhat consistent with previous findings that the more physically based sediment erodibility formulations in the Kok scheme enhance the dust sensitivity to climate variability relative to the Zender scheme (Kok et al., 2014a).

Results for the semiarid zone (Fig. 10c) are considered less robust due to significantly smaller dust-emitting areas (see Fig. 1). Overall, the contribution of wind speed further weakens, while hydroclimate drivers become more important. The resulting change of predictor relative importance, however, varies considerably. Specifically, four models (CESM2-CAM-Zender, E3SM3-Kok, EC-Earth3-AerChem and MPI-ESM1.2) retain the wind dominance, albeit with increased spatial variability. Hydroclimate dominance persists and strengthens in CESM2-CAM-Kok, INM-CM5.0, MIROC-ES2L and MIROC6, consistent with their behaviors in the arid zone. In contrast, the following models transition from wind- to hydroclimate-dominated regimes: E3SM2-Zender, CNRM-ESM2.1, CanESM5.0, CanESM5.1, HadGEM3-GC3.1, UKESM1.0 and MRI-ESM2.0. GISS-E2 models and IPSL-CM6A exhibit moderate increases of hydroclimate influence, resulting in roughly equal importance of wind and hydroclimate drivers. Compared to the ESMs, MERRA-2 generally produces dominant wind influence across all three climate zones.

To identify the sources of hydroclimate influence in the ESMs, Fig. 11 presents the median fractional contributions of five hydroclimate variables to the total R² in each model. The contribution attributed to specific humidity can largely be interpreted as a soil moisture effect, given the strong coupling between near-surface humidity and surface soil water content through evapotranspiration. In hyperarid regions, dust variability is expected to be dominated by wind speed, with minimal hydroclimate influence. Thus, anomalously large hydroclimate contributions in some ESMs may be explained by two possible mechanisms: (1) the model overestimates hydroclimate variability, thereby producing spurious correlations with dust emissions regardless of whether the hydroclimate variable is explicitly used in the dust parameterizations; and (2) the model reasonably represents hydroclimate variability but overestimates dust sensitivity to the hydroclimate variable, indicating a potential bias in the parameterization itself. Understanding the causes for the statistically inferred predictor influences would require detailed knowledge of the physical parameterizations and model configurations specific to each ESM.

Figure 11 Median factional contributions of hydroclimate drivers to the total explained variance (R²) in Earth system models and MERRA-2 over (a) hyperarid, (b) arid, and (c) semiarid climate zones. Hydroclimate variables are precipitation (P), soil moisture (SM), specific humidity (SH), air temperature (T), and leaf area index (LAI).

In the hyperarid zone, the hydroclimate influence in CESM2-CAM-Kok is associated with precipitation and specific humidity, which may partly reflect the increased dust sensitivity to soil moisture in the Kok scheme relative to the Zender scheme (Kok et al., 2014a). Another possible reason is the short simulation period for CESM2-CAM-Kok (2004–2013), which may not adequately capture the full range of dust variability and predictor relationships as in CESM2-CAM-Zender (1980–2014). In this regard, the E3SM experiments provide a more robust comparison between the two dust schemes. As shown in Fig. 11a, both E3SM2-Zender and E3SM3-Kok produce the expected negligible hydroclimate influence in the hyperarid zone. In the arid zone, however, E3SM3-Kok exhibits stronger hydroclimate influence than E3SM2-Zender, consistent with previous findings that the Kok scheme amplifies the dust sensitivity to hydroclimate conditions compared to the Zender scheme (Kok et al., 2014a). GISS models exhibit elevated influences from specific humidity and soil moisture, which explains their moderate wind contributions as shown in Fig. 10a.

In the arid zone (Fig. 11b), the enhanced hydroclimate influence is primarily associated with soil moisture and specific humidity in most ESMs, consistent with their well-established roles in modulating sediment erodibility (e.g., Csavina et al., 2014; Ravi et al., 2006; Kim and Choi, 2015). Several models – including CESM and INM-CM5.0 – attribute strong influences to LAI. Unlike other hydroclimate variables, LAI may be prescribed from climatology or interactively simulated by dynamic vegetation model components (Table 1). Models using prescribed LAI typically exhibit limited interannual variability in vegetation cover and therefore weak influence on dust emissions. For CESM and INM-CM5.0, the elevated LAI influence can be explained by their parameterizations of the bare soil fraction as a function of LAI.

4 Conclusions

This study examines discrepancies among 21 ESMs in representing the interannual variability of windblown dust emissions and the relative importance of near-surface wind speed versus five hydroclimate drivers (precipitation, soil moisture, specific humidity, air temperature, and LAI) across different climate aridity zones (hyperarid, arid, and semiarid). Recognizing the unconstrained, model-specific nature of dust emission fluxes, we use the dominance analysis technique to quantify the relative influence of a common set of six physical drivers within each model.

The extent of inter-model agreement in dust variability varies strongly with climate aridity. In the hyperarid zone, the ESMs exhibit poor agreement, with only 9 % out of 210 pairwise comparisons showing statistically significant positive correlations, reflecting large inconsistencies in model-simulated near-surface wind speeds. In arid and semiarid zones, the ESMs exhibit a dual pattern driven by a “double-edged sword” effect of land surface memory: models with coherent representations of hydroclimate variability tend to converge in their simulated dust variability, whereas those with divergent hydroclimate representations diverge in dust emission responses.

The relative importance of wind speed versus hydroclimate drivers also varies with climate aridity. In hyperarid areas, most ESMs capture the expected dominance of wind speed and minimal hydroclimate influence. CESM2-CAM-Kok exhibits elevated influences from precipitation and specific humidity, which may partly result from the more physically based sediment erodibility formulations in the Kok et al. (2014b) scheme relative to the Zender et al. (2003) dust emission scheme. A similar behavior is found in E3SM, where replacing the Zender et al. (2003) scheme with the Kok et al. (2014b) scheme reduces the wind dominance and enhances the hydroclimate influence on dust emission. Due to confounding model changes in the CESM and E3SM experiments, however, targeted experiments will be needed to isolate the effect of dust parameterization choice on the dust sensitivity. In arid and semiarid areas, wind influence generally weakens while hydroclimate influence strengthens in all ESMs. However, the relative importance between the two becomes less consistent, with contrasting model behaviors in retaining wind dominance or shifting toward hydroclimate dominance or near-equal importance. Compared to the ESMs, the MERRA-2 aerosol reanalysis generally produce stronger wind influence and weaker hydroclimate influence across all climate aridity regimes.

Note that the physical drivers considered in this study may not fully represent all the dust emission driving factors for specific emission schemes; instead, we focus on a common set of drivers for all models to provide a fair comparison across the ESMs. Therefore, the inferred relative importance from this analysis is limited to those common drivers considered and their influences on dust emissions in different models. Also, because of the statistical nature of dominance analysis, the predictor importance results shall be interpreted with caution when linking to model parameterizations. Despite these limitations, this study introduces a new framework for model intercomparison and yields new insights into how current ESMs represent the variability and physical controls of dust emissions across different climate aridity regimes. Overall, our findings highlight two promising directions for reducing model uncertainties in dust emission simulations: (1) improving the representation of near-surface wind speeds and gustiness in hyperarid regions, and (2) enhancing the representation of hydroclimate and land surface processes that modulate sediment erodibility in arid and semiarid areas.