Elucidating Hidden and Enduring Weaknesses in Dust Emission Modeling

Large‐scale classical dust cycle models, developed more than two decades ago, assume for simplicity that the Earth's land surface is devoid of vegetation, reduce dust emission estimates using a vegetation cover complement, and calibrate estimates to observed atmospheric dust optical depth (DOD). Consequently, these models are expected to be valid for use with dust‐climate projections in Earth System Models. We reveal little spatial relation between DOD frequency and satellite observed dust emission from point sources (DPS) and a difference of up to 2 orders of magnitude. We compared DPS data to an exemplar traditional dust emission model (TEM) and the albedo‐based dust emission model (AEM) which represents aerodynamic roughness over space and time. Both models overestimated dust emission probability but showed strong spatial relations to DPS, suitable for calibration. Relative to the AEM calibrated to the DPS, the TEM overestimated large dust emission over vast vegetated areas and produced considerable false change in dust emission. It is difficult to avoid the conclusion that calibrating dust cycle models to DOD has hidden for more than two decades, these TEM modeling weaknesses. The AEM overcomes these weaknesses without using masks or vegetation cover data. Considerable potential therefore exists for ESMs driven by prognostic albedo, to reveal new insights of aerosol effects on, and responses to, contemporary and environmental change projections.


Introduction
Mineral dust is central to many of Earth's systems (Shao et al., 2011).For example, dust can warm or cool regional climate, depending on the absorption and spectral properties of dust and background characteristics (Sokolik & Toon, 1999).These radiative effects depend on the volume of emitted dust and the mineral composition of atmospheric plumes, which differ over geographical source areas because particle size distribution and mineralogy vary spatially and temporally (Kok et al., 2017).Consequently, assessments of dust radiative effects rely on numerical models that simulate the emission, atmospheric transport, and deposition of the dust cycle (Mahowald et al., 2010) hereafter dust cycle models.Comparisons with atmospheric dust observations indicate that global dust cycle models include large uncertainties in simulated dust magnitude and geochemical properties (Huneeus et al., 2011).For example, global climate models used in the Fifth Assessment Report of the IPCC (IPCC, 2013) failed to reproduce observed North African dust over the second half of the 20th century challenging the validity of 21st century dust-climate projections (Evan et al., 2014).Large uncertainties and intermodel diversity remain and are larger than previous generations of models implying that modeled dust processes are becoming more uncertain as models develop (Zhao et al., 2022).It is common for global dust cycle models to be evaluated against atmospheric dust optical depth (DOD) and tuned to fit the observations from regional measurements, typically North Africa (Huneeus et al., 2011, p. 7809).However, this calibration approach does not enable an evaluation of the separate components of the dust cycle and specifically developments in dust emission modeling.Here, we are concerned that the evaluation of dust cycle models against atmospheric DOD has hidden weaknesses in the dust emission modeling, that parameterizations which attempt to improve dust emission are being falsely accepted, and consequently the approach to model development has become biased toward parsimony rather than achieving a balance with fidelity of dust emission processes (Raupach & Lu, 2004).
Dust emission schemes (Joussaume, 1990;Marticorena & Bergametti, 1995;Shao et al., 1996) were developed more than two decades ago and the underpinning, inextricably linked magnitude and frequency of sediment transport (Wolman & Miller, 1960) has not changed since.The magnitude of sediment transport driven by wind friction velocity is adjusted by the frequency of occurrence based on the wind momentum exceeding a critical sediment entrainment threshold (Wolman & Miller, 1960) causing highly dynamic, nonlinear responses over space and time (Raupach & Lu, 2004).When dust emission models were developed, there were few continuously varying global data sets available and simplifying assumptions were made for their implementation.The soil surface wind friction velocity to drive sediment transport (in the presence of large, typically vegetation, canopy roughness) was not available and instead the above canopy wind friction velocity was used.The partition between those drag forces used aerodynamic roughness lengths which were not available everywhere and therefore were set static over time and fixed over space to a bare soil surface condition (Zender, Bian, et al., 2003).Under these bare soil surface conditions, dust emission estimates were maximized and recognized as overestimated in the presence of vegetation.Consequently, dust emission schemes reduced estimates using E, the area of bare, exposed "erodible" soil surface (Marticorena & Bergametti, 1995).Varying with wind speed, drag was shown to be equivalent to shelter (Raupach, 1992), but the drag partition was parameterized using lateral cover, a static, geometric two-dimensional representation of canopy sheltering with the caution that mutual, three-dimensional, aerodynamic sheltering was not included in lateral cover (Raupach, 1992;Raupach et al., 1993).
For implementation of these approaches in dust emission modeling, researchers assumed for simplicity that dryland aerodynamic roughness, including nonphotosynthetic vegetation, was approximated by lateral cover from photosynthetic vegetation indices (VIs) readily available from satellite remote sensing (Evans, Ginoux, et al., 2016).These approaches using VIs therefore assumed that the sheltering effect of the drag was restricted only to that "green" canopy area and that the sheltered area of the soil surface did not vary with wind speed (was not aerodynamic).Furthermore, it was assumed for simplicity that the sediment entrainment was at the grain scale, static over time and a function of dry, loose and erodible spherical particle diameters discretized across sizes (Shao & Lu, 2000) so that soil texture data (typically aggregated over depth) could be used.Similarly, dust emission models assumed for simplicity that the soil surface had an infinite supply of dry, loose erodible sediment, despite soil surfaces, particularly in drylands, being sealed/crusted and/or with loose sediment occurring sporadically over space and intermittently over time (Vos et al., 2020;Webb & Strong, 2011).These simplifying assumptions in dust emission models represent the parsimony of implementation more than the fidelity of the dust emission processes.Although necessary for dust emission model implementation, the parameterizations need to be openly, transparently, and routinely re-evaluated particularly with 10.1029/2023JD038584 3 of 22 new technology, measurements, and different thinking.That dust models are becoming more uncertain (Zhao et al., 2022) indicates model development and/or its evaluation is not working well.By adjusting the magnitude of the dust cycle model estimates to atmospheric DOD, there is no possibility of recognizing weaknesses in dust emission modeling and therefore to recognize the need to change the modeling approach to compensate for any errors in the dust emission parameterization.The dust emission models provide no estimate of uncertainty, and it is therefore difficult to have confidence that dust emission models developed and evaluated in this way, adequately represent the geographical distribution of dust emission magnitude and frequency, particularly given the rather critical simplifying assumptions about sediment entrainment and sediment supply.This lack of confidence very likely explains the dearth of publications on dust emission model outcomes per se (as opposed to dust cycle model outcomes).Given the considerable advances in dust emission modeling over the last two decades, it is important to enable dust emission model outcomes to tackle these simplifying assumptions about sediment entrainment and sediment supply and the way in which dust emission models are evaluated.
Our aim is to show that there is an urgent and important requirement to routinely evaluate dust emission modeling separately and before the routine evaluation of dust cycle modeling against DOD.Our objectives to achieve that aim are to: (a) demonstrate that in a commonly implemented dust emission scheme, weaknesses have been introduced in its modeling which includes its parameterization and choices of data layers (Raupach & Lu, 2004);(b) show that these weaknesses in the dust emission modeling have been hidden in dust cycle models routinely calibrated against atmospheric DOD and that dust emission processes are not adequately represented in that calibration.Consequently, we expect our findings to support the implication elsewhere (Zhao et al., 2022), and our hypothesis that models have to their detriment, drifted from their original process fidelity toward modeled parsimony as forewarned (Raupach & Lu, 2004).For clarity, this work is less about which model is best (although a comparison is inevitable) and more about how dust emission modeling has been misled (weaknesses hidden) by the focus on calibrating dust cycle models to atmospheric DOD.For clarity, we recognize the valuable nature of calibrating dust cycle models to DOD.
It is timely that a new dust emission point source (DPS) database for southwestern North America has been collated and has been used with the albedo-based dust emission model (AEM) to circumvent the simplifying assumptions about sediment entrainment and sediment supply (Hennen et al., 2022(Hennen et al., , 2023)).We follow that established approach and evaluate dust emission modeling against dichotomous DPS observations and also compare those DPS to the frequency of atmospheric DOD.Given that we are investigating the evolved nature of traditional dust emission models (TEMs), and that many of its components are highly interactive, it is unreasonable to disaggregate the model components and/or make a superficial comparison of any one single component.Instead, we have produced an exemplar which represents TEMs in the view of the experienced large-scale dust modelers contributing to this study.We compare that exemplar TEM with the AEM which by design attempts to overcome the issues related to the soil surface wind friction velocity in the sediment transport equation (Chappell & Webb, 2016), and to redress the balance toward the fidelity of process representation (Raupach & Lu, 2004).

Dust Emission Modeling
Vegetation attenuates dust emission by extracting momentum from the wind and sheltering a portion of the downstream soil.By reducing wind speeds (U) at the soil surface, vegetation makes it more difficult for the wind to overcome the entrainment threshold for initiation of streamwise sediment flux (hereafter entrainment threshold) and consequent emission of dust particles by saltation bombardment and abrasion.Notably, the influence of vegetation sheltering is wind speed dependent (aerodynamic roughness) and both aerodynamic drag and partitioning of wind friction velocity between roughness elements and the soil, respond nonlinearly to changes in wind speed because of mutual sheltering.

Exemplar Traditional Sediment Transport
Calculation of the streamwise sediment flux density Q (g m −1 s −1 ) on a smooth soil for a given particle size fraction (d) on the particle size distribution (i) requires the above canopy wind friction velocity u * (m s −1 ) influenced by all scales of roughness at the Earth's surface, the air density ρ a (g m −3 ), the acceleration due to gravity g (m s −2 ), a dimensionless fitting parameter C and the bare, smooth (no roughness elements) entrainment threshold 10.1029/2023JD038584 4 of 22 of sediment flux u *ts (d) (m s −1 ).The original transport equation is typically rewritten in the dust modeling literature with the typographic correction and reformulated ratios which require a cubic term: Earth System Models (ESMs) or reanalysis wind field models over large areas (large pixels), with horizontal resolutions that are typically on the order of >10 km, use modeled wind speed at 10 m (U 10 ) to calculate the available above canopy u * .In recognition that vegetation exerts drag on the wind, u * must then be partitioned between the roughness elements (typically vegetation), and that available for driving flux at the soil surface (u s* ).The u *ts is adjusted by a function H(w) of soil moisture (w; kg 3 kg −3 ) (Fécan et al., 1999) and  = *  * (Raupach et al., 1993), the wind friction velocity ratio representing the roughness-induced drag partition (Marshall, 1971), modifies (Darmenova et al., 2009) the previous equation The u s* is required for sediment flux equations where u s* ≠ u * in the presence of any roughness canopy.In the absence of being able to estimate directly u s* , the u *ts H is divided by R for the model implementation to account for the drag partition making use of u * (Webb et al., 2020).Following this approach, this form (Equation 2) is incomplete because  3 * should be multiplied by R before it is cubed (Webb et al., 2020).However, values of R(z 0 ,z 0s ) are not known for all pixels and all time steps producing dust.
One of the common approaches to modeling dust emission in ESMs uses globally constant values of aerodynamic roughness length (z 0 ) (Tegen et al., 2002;Woodward, 2001;Zender et al., 2003).Here, we focus on the impact for large scale models and our exemplar TEM uses the incomplete formulation for Q TEM (Equation 2).Fixed aerodynamic roughness length for the landscape z 0 = 100 μm and the soil z 0s = 33.3μm fixes R(z 0 ) ≈ 0.91 to an almost bare land surface condition (see Equation S8 in Supporting Information S1).This approach assumes for simplicity that the Earth's land surface is devoid of vegetation roughness everywhere and static over time.With R(z 0 ,z 0s ) fixed,  (0,0s) * =  * is assumed, and u * changes only when the wind speed changes (see Equation S1 in Supporting Information S1).Consequently, the entrainment threshold u *ts H is not adjusted sufficiently by R(z 0 ,z 0s ) in the presence of vegetation and dust emission is overestimated for a given wind speed.To reduce dust emission in these non bare conditions, sediment transport is adjusted using a function E to compensate for the fixed R(z 0 ,z 0s ), hereafter we will drop the dependencies.
The function E was originally defined (Marticorena & Bergametti, 1995) as the ratio of bare exposed "erodible" surface area to total surface area.Only after the entrainment threshold is calculated, is the sediment transport and dust emission from a bare soil surface reduced by E < 1 typically in the presence of vegetation.However, sheltering is nonlinear over space and time since it depends on the mutual sheltering of the roughness (typically vegetation) structure, configuration and is aerodynamic i.e., sheltering changes with wind speed (Chappell et al., 2010).Nevertheless, it has become acceptable in the dust emission modeling literature to assume for simplicity that where A v is the non-aerodynamic planform area covered by vegetation.This E v is used in some ESMs so that leaf area index (LAI) or satellite "greenness" observations e.g., VIs can be used as a surrogate of the land surface fraction occupied by green vegetation (Evans, Ginoux, et al., 2016;Galloza et al., 2018;Sellar et al., 2019;Zender, Bian, et al., 2003).This parameterization crudely represents the aerodynamic nature of the sheltering (Chappell et al., 2010).Furthermore, the E v does not represent "brown" roughness, not readily evident in VIs caused by non-photosynthetic, dormant or dead vegetation, common in drylands which contain the majority of dust sources.The E v does not represent non-erodible stone-covered surfaces without sediment, also common in dryland regions.The inclusion of E in sediment transport models and its implementation using E v are prime examples of emphasizing parsimony over process fidelity (Raupach & Lu, 2004).Here, we will provide new insights for the impact of E and its implementation via E v on dust emission modeling.
Later publications in this approach provided practical implementations over large areas.However, these later approaches are strictly only valid for nonvegetated regions because of photometric volume scattering assumptions (Roujean et al., 1997).These later approaches are also typically fixed over time and therefore do not represent roughness change e.g., due to seasonality, land use change, and invasive species.Therefore, these estimates of z 0 do not overcome the challenge of estimating the spatiotemporal variation of the drag partition.Variability in the drag partition and u s* can be represented stochastically in models using probability distributions, as suggested by Raupach and Lu (2004).However, representation of the true heterogeneity (and uncertainty) of the drag partition has only recently been possible using field monitoring data (Edwards et al., 2022) and the approach is yet to be extended to regional or global scales.To retain the original dust emission scheme widely implemented in ESMs, and to show the outcome of this original approach, we do not apply these more recent approaches to our exemplar TEM.
In the dust emission modeling literature, there is little recognition of the uncertainty due to the inconsistency in model implementation scales.For example, the entrainment threshold (u *ts ) is calculated at the grain scale as a function of grain diameter, density, and interparticle cohesion (Shao et al., 1996).However, the above canopy u * is for an area when measured using a wind velocity profile or when using modeled wind data.Current dust emission modeling compares u * and u *ts (Equation 2) which assumes they are represented over (the same) area, which they are not.The threshold value at the grain (point) scale will not represent the required value over a large area (e.g., 11 km pixel) as established for large scale sediment transport modeling (Raupach & Lu, 2004) and recognized in other disciplines (de Vrese & Hagemann, 2016;Kyriakidis & Yoo, 2005;Raupach & Finnigan, 1995;Van Looy et al., 2017).For clarity, the inability of the grain (point) scale to represent area is, even assuming homogeneity within a pixel, caused by the inability to adequately represent the spatiotemporal variability in entrainment characteristics of the sheltered portion of the soil surface.Similarly, modeled wind speed over large e.g., 11 km pixels does not adequately represent the subgrid (pixel) scale heterogeneity of the aerodynamic roughness of the land surface.
The substantive remaining issues for TEMs that are currently known or evident from the literature include: 1. the incomplete form of  3 * in Q TEM (Equation 2) widely adopted in TEMs, overestimates sediment transport and dust emission, should have u * multiplied by the drag partition R, but the correct values of R are unknown (for every pixel and every time step producing dust); 2. poorly constrained aerodynamic roughness (z 0 ) and drag partition (R), causes R ≈ 0.91, fixed over space and time, to represent a bare soil surface which, in regions with any vegetation, underestimates u *ts and overestimates sediment transport and dust emission; 3. sediment transport in the presence of vegetation uses the "erodible" fraction E, implemented using dynamic vegetation cover (at nadir) which causes overestimated sediment transport and dust emission and E does not compensate for incorrect R; 4. the grain (point) scale of u *ts is incompatible with areal transport models which very likely underestimates sheltering and combined with the unreasonable assumption for drylands of an infinite supply of sediment, causes overestimated sediment transport and dust emission; 5. areal (e.g., >11 km) pixels of wind speed do not adequately represent the subpixel scale heterogeneity of the aerodynamic roughness and soil surface wind friction, which very likely causes underestimated dust emission which is scale dependent (different dust emission with different pixel size).
These weaknesses have existed for more than 20 years since dust emission schemes were first published.The first, incomplete formulation, only became evident recently by comparison with the albedo-based approach (Webb et al., 2020).That finding demonstrates the need for diversity in modeling approaches and provided the impetus for this broader investigation.The enduring influences of these weaknesses on large scale TEMs have been hidden for so long we contend, because dust cycle models are routinely evaluated against atmospheric DOD.

Albedo-Based Sediment Transport
In our AEM, the spatiotemporal variation in u s* is represented using the concept that aerodynamics of vegetation is proportional to sheltering (Raupach, 1992) and shadow (1-normalized albedo) (Chappell & Webb, 2016;Chappell et al., 2010).This albedo-based approximation of drag and its partition between forces was designed to provide an area-weighted value which was scale invariant (Chappell et al., 2018(Chappell et al., , 2019;;Ziegler et al., 2020) suitable for tackling nonlinearity of scaling (Raupach & Lu, 2004).This albedo-based approach enables direct 10.1029/2023JD038584 6 of 22 calculation of u s* given measurements of albedo from satellites (and many other ground to airborne sources), and enables the complete formulation for sediment transport and dust emission This Q AEM does not require E (or VI data), R, z 0 , or z 0s and thereby has four primary parameters less than the exemplar TEM, and removes their associated sources of uncertainty.Soil moisture w was obtained from the same data as used in the TEM and calculated in the same way.The u s* is obtained directly from ω ns , using shadow from the normalized and rescaled albedo which describes the area-weighted land surface aerodynamic structure (partitioned between above canopy and soil surface) independent of waveband, making it highly suitable for the inclusion of dryland non-photosynthetic material The is a coupled parameter which describes how the soil surface wind friction velocity is dependent on the wind speed at a given height (h), where that height is ideally at freestream.This approach enables an AEM (see Equations S10-S13 in Supporting Information S1 for a full description of the implementation not limited to MODIS albedo products) (Chappell & Webb, 2016;Chappell et al., 2010Chappell et al., , 2018;;Ziegler et al., 2020).This approach assumes that the wind profile is logarithmic with neutral buoyancy, in common with the approach typically used in large-scale dust emission modeling (Marticorena & Bergametti, 1995).
Typically, wind fields over large (e.g., >11 km) pixels (e.g., ERA5-Land) are estimated at a blending height (e.g., 40-m height) and then use z 0 values, assumed homogeneous within few large land cover types, to extrapolate wind speed to typically 10-m height (ECMWF Forecast User Guide).Therefore, these modeled 10-m wind speeds do not adequately represent the subgrid scale heterogeneity of aerodynamic roughness (z 0 ) which causes wind speeds to be underestimated.The albedo-based approach improves constraints on the coupled parameter   *  ℎ and scales linearly over space and time.The approach therefore offers considerable potential to make u s* estimates over large areas (upscaled) or to improve the subgrid scale heterogeneity of z 0 and therefore improve the downscaling of 10-m wind speed.However, improving the aerodynamic scaling of wind speed is beyond the scope of this study.
Notably, the AEM inherits the long-established and poorly constrained entrainment threshold u *ts which at the grain (point) scale is inconsistent with the new area-weighted albedo-based approach.The AEM also retains the unreasonable assumption in drylands of an infinite supply of sediment.Consequently, we tackle this scaling inconsistency and source of uncertainty in the AEM, by removing the frequency distribution of the entrainment threshold (right-hand side of Equation 3) and replacing it with the frequency distribution from satellite observed DPS data which retains the fidelity of the dust emission process (Section 2.4).

Dust Emission Modeling
The vertical dust mass flux (F; g m −2 s −1 ) may be calculated from Q using physically based schemes (Kok et al., 2014;Shao et al., 1996).However, one of the common approaches in regional and global applications, and that used here for the exemplar TEM and AEM, calculates F as an empirical function of Q (Marticorena & Bergametti, 1995) where Q differs between TEM and AEM only because of the sediment transport equations above The dust emission parameterization considers the emission flux to be driven by saltation bombardment, with the intensity proportional to Q, and the soil's clay content (clay % typically <2 μm fraction of soil particles at the soil).The mass fraction of clay particles in the parent soil was allowed to vary over space, but was fixed over time.In both models, we used the latest, reliable and spatially varying layer of particle size (Dai et al., 2019) and restricted clay % to a maximum value of 20% consistent with previous work (Marticorena & Bergametti, 1995) which showed reasonable results when applied in a regional model calibrated to DOD (Woodward, 2001).The Q, which produces dust, is adjusted by the emitted dust fraction M for a given particle size fraction with diameter d which we calculated as 1 < d < 10 μm following Zender, Bian, et al. (2003) by using M = 0.87.When the soil is covered by snow it is unable to provide any dust emission.In this situation, it is most effective to use a mask which determines whether snow is present or absent (A s ).Similarly, if the soil is bare but frozen it is unable to release sediment almost regardless of how much wind energy is applied.In this situation, it is most effective to use a mask which determines whether the soil is frozen or not (A f ).
Some models also use geographically preferential dust sources that limit the magnitude of dust emission (Evans, Ginoux, et al., 2016;Ginoux et al., 2001;Mahowald et al., 2010;Tegen et al., 2002;Woodward, 2001;Zender, Bian, et al., 2003).In our exemplar TEM, we do not use preferential dust sources to make clear in our results the cause of differing dust emission magnitude and frequency.Also, the emission in ESMs is typically "tuned" down to match observed atmospheric DOD (Zender, Bian, et al., 2003).Here, we do not apply this final global tuning.

Large-Scale Dust Emission Modeling, Mapping Spatial Variation, and Change Detection
To implement vertical dust emission, we used contemporary (2001-2020) Earth observation data including spatially and temporally varying wind speeds (at 10 m height), soil moisture (0-7 cm depth), and soil surface temperature (to represent frozen ground which inhibits sediment flux) from the latest ERA5-Land data (Muñoz-Sabater et al., 2021) (hourly; 0.1°).The coverage of snow in a given pixel is an areal quantity and therefore ranges between 0 and 1.Consequently, we applied the MODIS Normalized Difference Snow Index (Hall et al., 2016) (MOD10A1 from Terra, daily at 500 m).We used soil surface temperature available in ERA5-Land and set a threshold of 273.15K above which sediment flux can occur.The use of these data does not imply priority over any other data.We recognize that there are different qualities to different model data, as evident in their wind fields (Fan et al., 2021).Where applicable, we used the same data in both the exemplar TEM and AEM to consider the relative differences.We used the exemplar TEM with R(z 0 , z 0s ) ≈ 0.91 fixed over space and static over time.Following the current practice, we calculated u * from the modeled 10 m height wind velocity using the logarithmic layer profile theory and aeolian roughness length (Darmenova et al., 2009) (details are provided in the Supporting Information S1).We allowed soil moisture to vary in the same way in both TEM and AEM.Only in the exemplar TEM was MODIS Normalized Difference Vegetation Index (NDVI; MOD09GA Collection 6) data used to calculate the bare soil fraction E. For comparison, we used the AEM with soil surface wind friction velocity u s* /U h obtained directly from MODIS albedo (MCD43A3; Collection 6) varying daily, every 500-m pixel across the study area.MODIS is aboard polar-orbiting satellites which cause incomplete coverage.However, the variation in roughness at the daily-weekly scale is so small that we were able to smooth the available data to improve the coverage.Soil clay content was represented with a digital soil texture map (Dai et al., 2019) and used in both models (see Section 2 and Text S1 and S2 in Supporting Information S1).
All data are available from the catalog of the Google Earth Engine (GEE; Gorelick et al., 2017) which then required no downloading and reformatting.We used the Java script coding environment to calculate daily dust emission (kg m −2 yr −1 ).Given the availability of DPS validation data at sites in southwestern USA, we restricted our mapping to North America including dust source regions bordering the USA.This regional focus enabled the spatial patterns and changes over time to be readily visualized by contrast to global maps.Testing the code and visualizing the results for smaller time periods across the study area was almost instantaneous in the GEE.Data processing at 500 m and using daily resolution between 2001 and 2020 across North America took typically <12 hr.These data were exported from the GEE for the calibration/validation in a Python coding environment and images (TIF) from the GEE were also exported for presentation using ArcGIS Pro.
At the sites and days when dust was observed using DPS, we compared them with the dust emission produced by the exemplar TEM, AEM, and atmospheric DOD frequency.For the year 2020 (the most recent year of complete data available at the time of the study) and the main dust emission months of March-May (MAM) in North America (Hennen et al., 2022), we analyzed the spatial variation of the main controlling variables (wind and aerodynamic roughness) and dust emission produced by the exemplar TEM and AEM.The dust emission of both models was restricted to wind speeds between 8.5 and 9.5 m s −1 to emphasize the difference in our modeling approaches below (Figures 2 and 5), which would otherwise be hidden by taking the average for all wind speeds.
CHAPPELL ET AL. 10.1029/2023JD038584 8 of 22 Finally, we also map the difference in driving variables during MAM for the year 2001 compared with the year 2020.This large difference between years provides the greatest opportunity to appreciate the impact on dust emission of any change in land surface conditions and which has not been demonstrated previously in dust emission modeling applications.The dust emission on dust days was used to obtain the mean difference.That mean difference is then tested for significance using the minimum detectable change (MDC) framework (Webb et al., 2019;Woodward, 1992) and the results are displayed.The MDC was established using critical values for false acceptance and false rejection (α = 0.05; β = 0.05, respectively) and the change in dust emission which did not exceed the MDC, was set to 0 (not detectable = white).Details of how the MDC was calculated are described in Text S3 in Supporting Information S1.

DPS and DOD Frequency
Commonly, DOD from ground-based or large area Earth observation data are typically used to evaluate the performance and/or calibrate dust cycle model simulations (Meng et al., 2021).This approach assumes for simplicity that: (a) atmospheric DOD adequately represents the dust emission process and (b) the spatial variation in magnitude and frequency of dust emission in the dust cycle model is correct.However, we know a priori that dust in the atmosphere is only partially related to dust emission because dust concentration is controlled by dust emission magnitude and frequency which varies over space and time, by residence time of dust near the surface which itself is dependent on wind speed (Textor et al., 2006), and on dust deposition in the dust source region, a size dependent process (Mahowald et al., 2014).To understand the extent to which DOD estimates the spatial variation in dust emission magnitude and frequency we calculated the probability of dust occurrence retrieved from the DOD (DOD > 0.2) using the criteria established previously (Ginoux et al., 2012).We note the stated limitations of DOD to be largely restricted to bright land surfaces in the visible wavebands which implies reduced performance over areas where vegetation is present (Ginoux et al., 2012).We demonstrated below (Text S4 in Supporting Information S1) that there is little impact of the chosen DOD threshold on the results presented here.To calculate DOD, we used wavebands available from monthly Moderate Resolution Imaging Spectroradiometer (MODIS; MOD08 M3 V6.1 Deep Blue L2 Aerosol Product) at a 1° pixel resolution (Platnick et al., 2015).We used this resolution because it was consistent with the largest DPS data of our larger study.We also believe that this 1° pixel resolution provides the best opportunity for DOD to be associated with dust emission.We note the difference between our use of collection 6.1 and collection 5.1 used by Ginoux et al. (2012).We assume that the criteria we used following Ginoux et al. (2012) are applicable to the later collection.The DOD was retrieved from those pixels in which dust emission was observed from DPS in space and time throughout 2001-2016.
The identification of DPS data is a highly time-consuming and labor-intensive activity (Text S5 in Supporting Information S1).Consequently, there are few (published) studies relative to the large number of global dust source regions.Here, we use DPS data collated from several previous studies in North America (Baddock et al., 2011;Kandakji et al., 2020;Lee et al., 2012).Those studies identify the point source of dust emissions in New Mexico and Texas between 2001-2016 and 2001-2009 and for 2001-2009 in the Chihuahuan Desert and New Mexico, collating a single data set of DPS data from North America.The DPS observations were identified using MODIS data with visible to thermal infrared wavebands (0.4-14.4 mm; see Text S5 in Supporting Information S5).Modeled (AEM and exemplar TEM) and observed frequencies are aggregated by a 1° × 1° grid matrix, normalizing the results to the lowest resolution MODIS DOD data (Figure 1).
This aggregation is performed to tackle the incompatibility of different scales (Gotway & Young, 2002).At the point scale, there is considerable unexplained variance which is likely related to the DPS data location uncertainty of around ±2 km (Kandakji et al., 2020) due to the phase difference between timing of dust emission and availability of the imagery.The unexplained variance and incompatible scales are well-established in the geostatistical literature (Gotway & Young, 2002).We reduced the unexplained spatial variance by aggregating the DPS data into a 1° grid system, matching the horizontal resolution of the DOD data.A binary value is applied to each DPS location for each day during the observation period (length of the DPS study) where 1 = dust observed, 0 = no dust observed.For comparison, the same process is applied to modeled and DOD data, by spatially interpolating raster images at the native resolution of the original outputs (DOD = 1°, AEM/TEM = 0.1).Aggregation of daily probability for each grid box frequency provides either 1 or 0 observations per day for all measurements CHAPPELL ET AL. 10.1029/2023JD038584 9 of 22 (AEM/TEM/DOD/DPS).For a daily grid cell to record dust (frequency = 1), at least one of the DPS locations must record the emission (F > 0)/presence (DOD > 0.2) of dust.Where emission/presence of dust is identified at more than one DPS location (within the same grid cell) the grid value remains one-this is the maximum.Finally, the sum of these frequencies is divided by the number of years of observation, to provide an annual probability of dust emission (as shown in Figure 4).This normalization is necessary due to varying periods of observation across the DPS studies.

Dust Emission Model Evaluation Against DOD and Calibration Against DPS
We compared DPS observed occurrence with modeled dust emission determined by the exemplar TEM and the AEM.
Similarly, during those same DPS observed occurrences we compared the retrieved estimates of DOD frequency.For all of those model estimates of dust frequency (DOD, exemplar TEM, and AEM), separately we fitted log-linear regression models which produced regression model parameter coefficients, R 2 correlation and the square root of the sum of squared difference (SSE) between DPS and model predictions to form the RMSE = √SSE/(N − df) where N number of data are adjusted by the degrees of freedom (df = number of independent dust emission model parameters).
We improve the constraints   * > *  on dust emission model evaluation by calibrating the dust emission magnitude according to modeled emissions during those observed occurrences.We follow the established approach (Hennen et al., 2022) by using observations of dust emission frequency at dust emission sources during satellite observations.We use dichotomous satellite observed DPS data and its probability of occurrence P(DPS > 0) as a first-order approximation of the probability of sediment transport P(Q > 0) leading to the proportion of dust (F) emission P(F > 0) at those points identified to produce dust .However, we know a priori that there are at least two very weak assumptions in the dust emission modeling: that the soil surface is covered homogeneously with an infinite supply of loose erodible material, which when mobilized by sufficient wind friction, causes dust emission.This approach assumes energy-limited dust emission which is rarely justified in dust source regions where the soil surface is rough due to soil aggregates, rocks, or gravel, sealed with biological or physical crusts, or loose sediment is simply unavailable (Vos et al., 2020;Webb & Strong, 2011).Consequently, we follow the recently established approach (Hennen et al., 2022) and bypass those weak assumptions by using observed dust emission frequencies at DPS data locations to parameterize the entrainment threshold frequency distribution using the established for this region (Hennen et al., 2022) Log 10 (AEMcal) =0.88 Log 10 (DPS) −2.02 where AEM cal is the adjustment of modeled F DPS values using the calibration.This approach calibrates our large-scale dust emission model to be consistent with DPS.This approach overcomes the currently poor model constraint   * > *  of the sediment transport threshold.The DPS data are likely biased away from the smallest dust sources which may not appear or are difficult for operators to detect using optical reflectance (Urban et al., 2018).Nevertheless, under these conditions, the AEM cal provides precise and accurate maps of seasonal dust emission, temporal dynamics, and mean regional dust emission.We applied this calibration only to the AEM change over space and time (Section 3.3) to provide valid, calibrated dust emission estimates for comparison with the exemplar TEM.

The Impact of "Erodible" Fraction (E) Implementation on Dust Emission Modeling
We simulated dust emission with wind speed varying between 0 and 12.5 m s −1 (Figure 2a).The exemplar TEM dust emission is shown with a fixed aerodynamic roughness length for the landscape scale z 0 = 100 μm and the soil scale z 0s = 33.3μm following several previous studies e.g., (Zender, Bian, et al., 2003), which fixes R ≈ 0.91 and assumes for simplicity that the land surface is almost devoid of vegetation roughness and static over time.With E = 1, sediment transport and dust emission are unadjusted and increases along the upper (large dashed) curve as wind speed increases (because wind friction is fixed).When the land surface is partially covered in vegetation and E = 0.5 (all other conditions remain the same), sediment transport and dust emission increases as wind speed increases but at a consistently reduced rate (solid line).This separate curve and reduced rate are caused entirely by using E. The open square is the exemplar TEM at 8 m s −1 and the filled square is the exemplar TEM at 9.2 m s −1 .Despite E being implemented to reduce sediment transport in the presence of vegetation, the (unintended) outcome is very similar dust emission for similar wind speeds (from open square to filled square).For clarity, this finding reveals for the first time that the exemplar TEM will emit larger amounts of dust from vegetated surfaces than from bare surfaces (for the same erodibility conditions).This is because E does not adequately reduce dust emission depending on the interplay between wind and vegetation.Consequently, where there are large wind speeds in vegetated regions, the exemplar TEM will incorrectly produce large amounts of dust emission.
The AEM (uncalibrated) for a smooth, unvegetated situation (   * ∕ℎ = 0.035; dotted line) produces larger dust emission than the exemplar TEM for the same 8 m s −1 wind speed (open triangle; Figure 2a).However, in a 10.1029/2023JD038584 11 of 22 rough, vegetated situation (   * ∕ℎ = 0.022) dust emission declines to almost zero, along the same curve.Despite the larger wind speed of 9.2 m s −1 (closed triangle), the rough surface causes the surface wind friction velocity to decrease, barely exceeding the entrainment threshold, and consequently dust emission is considerably reduced.The increase in roughness is sufficient to overcome the increase in wind speed and causes dust emission to be much smaller.The interplay between wind and roughness realistically produces accurate and precise soil surface wind friction velocity essential for reliable and consistent dust emission estimates across complex terrain including vegetation.
These findings are partially expected based on the theory described above in Section 2. However, the impact of E in the exemplar TEM has not previously been recognized and is only evident here as unusual, relative to the more physically based AEM.The exemplar TEM is driven by wind speed attenuated by aerodynamic roughness, but which is here fixed over space and static over time to a bare soil surface, and dust emission is subsequently reduced by E based on vegetation cover.Consequently, wherever and whenever wind speed exceeds the entrainment threshold, the exemplar TEM will produce sediment flux and dust emission.To illustrate this point, Figure 2b shows change in dust emission with change in soil surface wind friction velocity normalized by wind speed (   * ∕ℎ ).In other words, Figure 2b shows how dust emission changes in either space and/or time as roughness changes in the AEM or as E changes in the exemplar TEM.Since the influence of wind speed is removed on the x-axis (and wind friction is fixed in the model), exemplar TEM produces no change for a given wind speed of e.g., 10 m s −1 .The cause of change in the TEM for 10 m s −1 (solid red line) is due solely to the value of E varying.Since E is not aerodynamic (does not change with wind speed) and wind friction is fixed, dust emission does not change except when E changes.Under a scenario with the wind speed reduced from 10 to 8 m s −1 , the exemplar TEM F increases monotonically but at a reduced rate; that rate does not change with roughness (   * ∕ℎ ) it changes only with E. Similarly, when the wind speed increases from 10 to 12 m s −1 , the exemplar TEM F increases monotonically at an increased rate, but that rate does not change with roughness (   * ∕ℎ ).In contrast, for the wind speed of 10 m s −1 , the AEM produced a greater reduction in dust emission than the exemplar TEM for the greatest decrease in   * ∕ℎ (Figure 2b).With the greatest increase in   * ∕ℎ , the AEM produced a larger increase in dust emission than the exemplar TEM.When wind speed is consistently reduced to 8 m s −1 , the change in dust is smaller than that at 10 m s −1 .Notably, there is no change in dust emission between a change of −0.01 <   * ∕ℎ > 0.01 (Figure 2b).When wind speed is consistently increased to 12 m s −1 , the change in dust emission produced by the AEM is large, continuous, and evident as   * ∕ℎ changes.The results of these simulations illustrate how the exemplar TEM does not adequately represent vegetation sheltering dynamics and that E compensates by adjusting the magnitude, not the onset of dust emission.The exemplar TEM will emit similar amounts of dust from vegetated surfaces as from bare surfaces (for the same conditions).These weaknesses in E have not been apparent previously, despite considerable use and application particularly in large-scale dust models, because dust emission has not been isolated in previous evaluations.In contrast, the AEM provides a direct estimate of   * , which modifies dust emission as roughness and/or wind speed changes.Since this direct estimate of   * is available from albedo, from ground measurements, monitored from satellite remote sensing, or modeled prognostically in ESMs, it is available over space and/or time without the need for R or the bare soil fraction E, thereby reducing uncertainty in the model parameterization.We elaborate the exemplar TEM weakness in comparison with the AEM by modeling dust emission change over space and time (see Section 3.3).

Modeled and Observed Dust Emission Frequency at DPS Locations
We retrieved atmospheric DOD probability P(DOD > 0.2) at previously identified DPS locations across areas of southwestern North America to compare with DPS observed probability P(DPS > 0) (Figure 3).The P(DOD > 0.2) showed little resemblance to P(DPS > 0), with a distinctly different spatial pattern and considerably greater probability in some areas.Peak P(DOD > 0.2) occurred across the USA/Mexico border in the Chihuahuan Desert, while P(DPS > 0) peaked over the Southern High Plains in eastern New Mexico and western Texas.The P(DOD > 0.2) probability increases in areas of reduced vegetation roughness (Figure 1) as difficulties in measuring atmospheric dust over dark surfaces (e.g., vegetation), limit the DOD frequency data to only the most arid regions.In areas where the data are comparable (e.g., northern Chihuahuan Desert; 108°-104°W, 29°-32°N), P(DOD > 0.2) is (at least) an order of magnitude greater than DPS.We compared estimated dust emission frequency (AEM and exemplar TEM models with F > 0 or DOD > 0.2) with observed DPS frequency (in days per year) at each DPS grid box (see in Figure 1).For each model comparison, the observed DPS frequency remained the same, with differences in the model described on the x-axis (Figure 4).At most grid boxes, modeled frequency exceeds observation, consistent with the discrepancies between the grain (point) scale of the entrainment threshold, the area-weighted wind friction, and the areal (11 km) scale of wind speed.Both AEM and exemplar TEM overestimate dust emission frequency with RMSE (Log 10 ) = 0.6 and 0.76 (4 and 5.8 days per year), respectively, relative to the 1:1 line (Figure 4) demonstrating slightly improved performance by the AEM.Nevertheless, across all grid box data, the relation between DOD frequency and DPS was very large exceeding DPS frequency by nearly 2 orders of magnitude, with RMSE (Log 10 ) = 2.09 (123 days per year), considerably larger than the relation between DPS and the dust models.Least squares log-linear regression models were fitted to all models, with AEM and exemplar TEM frequencies showing statistically significant correlation with DPS observed frequency, producing a regression slope of 0.74 (AEM), 0.76 (TEM), and R 2 = 0.80 (P << 0.001).The DOD frequency did not show a statistically significant correlation with DPS observed frequency, with a regression slope of −0.12 and R 2 = 0.02, (P = 0.35).

Modeling Dust Emission Change Over Space
The mean albedo-based u * /U h and full range of U 10 for the year 2020 are shown (Figures 5a and 5b).Since the exemplar TEM has fixed aerodynamic roughness, its wind friction velocity is fixed and hence varies with wind speed and E. For consistency with Figure 2, and to isolate the influence of E, the mean dust emission is shown for selected wind speeds (U 10 = 8.5-9.5 m s −1 ) for both AEM cal and exemplar TEM (Figures 5c and 5d).Consistent with Figure 2, the spatial distribution of mean dust emission was very different between AEM cal and exemplar TEM in both magnitude and spatial extent.According to AEM cal , large dust emissions (0.05-0.12 kg m −2 yr −1 ) occurred in discrete areas across the Southern High Plains (104.5°W,33.5°N), northern Chihuahuan Desert (107.5°W,32°N), southwest Colorado Plateau (110.5°W,35°N), and the Great Divide Basin within the Wyoming  (Baddock et al., 2011;Kandakji et al., 2020;Lee et al., 2012).For each point, the y-axis represents the observed number of DPS observations (per grid box) per year during different observation phases of the DPS data sets within the observation time period (2001)(2002)(2003)(2004)(2005)(2006)(2007)(2008)(2009)(2010)(2011)(2012)(2013)(2014)(2015)(2016).For AEM and exemplar TEM, the x-axis describes the number of modeled observations (F > 0) at DPS locations in each grid box per year during the same time period (x-axis).The x-axis describes the frequency that DOD > 0.2 per year for the same period.The least squares logarithm regression of modeled against DPS observations produced the model parameter coefficients, R 2 correlation and the square root of the mean squared difference between DPS, and model predictions (RMSE) adjusted by the degrees of freedom (df) using the number of model parameters (df = 9 for AEM; df = 12 for TEM; df = 6 for DOD frequency).
CHAPPELL ET AL. 10.1029/2023JD038584 14 of 22 Basin (108.5°W,42°N).These areas correspond with small roughness u * /U h , and large wind speed.Furthermore, large u * /U h reduces the influence of large winds and restricts dust emission.Although the exemplar TEM dust emission occurred with similar maximum magnitude to the AEM cal , the exemplar TEM dust emission was distributed over a much larger area, including large parts of New Mexico and Wyoming, while also extending through the Great Plains in northwest Texas, Oklahoma, Colorado, and Nebraska (Figure 5d).This pattern of exemplar TEM dust emission matches closely the spatial distribution of mean wind speed (Figure 5b) controlled by E. However, consistent with the results of Figure 2, these results show that despite the variation in E implemented with dynamic vegetation cover, dust emission remains very large and similar over space, producing relatively homogeneous dust emission (dark tones) despite considerable land surface heterogeneity (Figure 5d).These results demonstrate that the implementation of E using dynamic vegetation cover, does not compensate adequately for poorly constrained, fixed R (which represents bare soil surface with an infinite supply of sediment).

Modeling Dust Emission Change Over Space and Time
Separate differences in albedo-based roughness (u * /U h ) and wind speed at 10-m height (U 10 ) for the year 2001 compared with the year 2020 and greater than the MDC significance (P < 0.05), were produced (Figures 6a  and 6b).Statistically significant change in roughness across North America occurred with a range ±0.01.The greatest reduction (<−0.01) in roughness occurred in Canada and was very likely caused by changes in the duration of snow coverage.Note that snow is removed from u * /U h when calculating dust emission.South of the USA/Canada border, roughness reduced (−0.01)across large areas of Montana, the Wyoming Basin, and eastern parts of the Great Plains (Colorado, Kansas, and Nebraska).Smaller reductions in u * /U h (−0.01 to −0.005)  10) and without (d; exemplar TEM; Equations 2 and 5) varying aerodynamic roughness.The dust emission displayed is for wind speeds restricted to between 8.5 and 9.5 m s −1 (for comparison with Figure 2).The daily maximum wind speed, described in hourly data from ERA5-Land (Source: ECMWF) are used in both models.
occurred in discrete areas of the Southern High Plains, and northern Chihuahuan Desert.The greatest increase in u * /U h (>+0.01)occurred across the American Mid-West, including Minnesota, Iowa, and South Dakota.In dusty areas (Figure 5), the greatest increase in u * /U h (+0.005 to +0.01) occurred as discrete locations within the Chihuahuan and Sonoran Desert, the Great Basin (Nevada), and the southern limit of the Southern High Plains (eastern New Mexico and western Texas).Mean U 10 changed with a range ±1.6 m s −1 , with the largest increase (>1.6 m s −1 ) across southwest USA, including the Great Basin, Mojave and Sonoran Deserts, and the Colorado Plateau.Mean U 10 reduced (<−0.8m s −1 ) in the Mid-West states of Wisconsin and Illinois.
Differences in mean dust emission during the peak dust season (MAM) for the year 2001 compared with the year 2020 and greater than the MDC significance (P < 0.05), were produced for both exemplar TEM and AEM cal (Figures 6c and 6d).Statistically significant change in dust emission comparing AEM cal and exemplar TEM varied across the range ±2 kg m −2 yr −1 .The AEM cal produced a significant decrease in dust emission (−1 to −2 kg m −2 yr −1 ) from several areas, including the Southern High Plains (eastern New Mexico and western Texas), the Colorado Plateau, and the Sonoran Desert (Figure 6c).The AEM cal showed a significant increase in dust emission from the Wyoming Basin, and discrete locations in Montana, and western areas of the Great Plains (west Colorado, Nebraska).In contrast, where no change in the AEM cal was detected, the exemplar TEM produced a significant decrease of dust emission across large areas of the Great Plains (up to −2 kg m −2 yr −1 ), the arid southwest (−1 to −2 kg m −2 yr −1 ), including the Mojave, Sonoran, and Chihuahuan Deserts, and the Mid-West (−1 to −2 kg m −2 yr −1 ).Exemplar TEM dust emission increased significantly across the Wyoming Basin (up to 2 kg m −2 yr −1 ), the Great Basin and northern Mexico (Figure 6d).The exemplar TEM shows considerable false change in dust emission relative to the AEM cal which is calibrated to DPS frequency distributions.These findings have considerable implications for the use of the exemplar TEM in large-scale dust climate projections of ESMs.10) varying aerodynamic roughness (c) and with exemplar TEM (Equations 2 and 5) z 0 fixed and static over time (d).Wind data are from ERA5-Land (Source: ECMWF).See Text S3 in Supporting Information for details on the calculation of the MDC.

Overcoming Dust Emission Model Weaknesses Using the Albedo-Based Approach
Dust emission modeling has struggled to replicate observed dust emission magnitude and frequency, indicating an inability to adequately represent soil wind friction velocities (Evan et al., 2014).Many of the TEMs assume homogenous bare ground, before using the complement of dynamic vegetation cover to reduce emission.Using satellite DPS (Figure 1), we have shown the exemplar TEM overestimates dust emission frequency by 0.76 of an order of magnitude (RMSE = 0.76 using log 10 ) (Figure 4).Using albedo to describe variability in aerodynamic roughness through changes in vegetation structure, the AEM performs theoretically better (Figure 2) at correctly estimating the probability of u s* exceeding the entrainment threshold, and subsequent changes in dust emission timing and magnitude.When compared to observed DPS (Figure 4), AEM performs only moderately better than the exemplar TEM, still overestimating dust emission frequency by 0.6 orders of magnitude (RMSE = 0.6 using log 10 ).However, the AEM does not use z 0m , z 0s , R, or E, and the monitored normalized shadow is calibrated to wind tunnel u * /U h .In contrast, the exemplar TEM is pre-adjusted to values of z 0m and z 0s for bare soil surfaces in R which are fixed over space and static over time and then adjusted by E. Furthermore, most DPS used in this work are from predominantly barren and windy environments, with mean u * /U h = 0.069 and mean U 10 = 6.9 m s −1 , reducing the potential influence of dynamic vegetation.Nevertheless, the overestimation of modeled dust emission relative to the observed frequency, occurs because of one or more of the main priority factors described in Table 1.Those factors are elucidated in terms of potential impact on dust emission modeling.
We use the latest version of ERA5-Land wind (U at 10-m height; U 10 ) data at a reasonably fine (11 km) resolution.It is evident that U 10 is overestimated in some global regions (Fan et al., 2021).However, there appears to be no systematic bias in the global wind fields that would lead to the systematic overestimation of dust emission frequency.The grain scale of u *ts is evidently incompatible with areal dust emission modeling and this factor appears to be the most likely cause of the overestimated model dust emission frequency and should be a priority for future work.Without resolving the scale of u *ts it is not possible to isolate the impact of the assumed infinite supply of loose erodible material (Table 1).However, the scale invariant nature of the albedo-based approach (Ziegler et al., 2020) holds considerable potential for tackling these long-standing and widely omitted scaling issues in dust emission modeling.It is very likely that these two factors explain the first-order differences between the DPS frequency and the dust emission model frequency.Although we have reduced uncertainty by using grid boxes for the DPS frequency, there remains uncertainty over the use of DPS frequency (Urban et al., 2018).However, by comparison with DOD frequency, the use of DPS frequency is up to 2 orders of magnitude smaller indicating that dust emission models should be evaluated against DPS data.
Beyond the observed DPS, vegetation roughness appears influential, constraining dust emission greater than 0.1 kg m −2 yr −1 to areas where u * /U h is no greater than 0.06, even during periods of peak (8.5-9.5 m s −1 ) wind speed in our case.In contrast, the exemplar TEM predicts dust emission >0.1 kg m −2 yr −1 in areas where u * /U h is greater than 0.075 (vegetated), including large areas of the vegetated Great Plains.This difference is emphasized in parts of western Oklahoma (99.5°W, 35.5°N),where mean u * /U h > 0.08 prevents dust emission from the AEM, despite a mean U 10 > 7 m s −1 .However, in those areas, exemplar TEM dust emission exceeds 0.2 kg m −2 yr −1 .These contrasting estimates emphasize exemplar TEM dependency on variability in U 10 , due to the incomplete modeling using  3 * (Webb et al., 2020) and the inability of R adjusted by E, to correctly attenuate wind speeds by aerodynamic roughness (Table 1).The proposed solution is simple for contemporary, satellite era analysis, remove E (and the need for R, z 0 , and z 0s ) and make use of readily available satellite albedo data to directly estimate   * .The limitations of the approaches and their impacts on dust emission modeling have been hidden for more than two decades since dust emission models were developed.The limitations have created two further issues, (a) a requirement for post-process tuning using dust source maps, which limits the model's ability to predict dust without a priori information; (b) large scale uncertainty and scale dependence driven by a large spatial and temporal variability in U 10 .
Our comparison of dust emission between two time periods (2001 and 2020) emphasizes a previously unrealized impact of varying aerodynamic roughness in the temporal variability of dust emission magnitude.Through the calculation of dynamic u s* , the AEM constrains dust emission to relatively small areas, restricting significant variability between time steps, to only dust producing areas (e.g., the arid southwest and semiarid parts of the Great Plains; Figure 6c).In contrast, the exemplar TEM's dependency on U 10 variability shown here, produces significant changes in dust emission over vast vegetated areas, including those which are very unlikely to produce dust (e.g., temperate areas of the Great Plains and the grasslands of northern Mexico; Figure 6d).These results demonstrate for the first time, that the exemplar TEM using E implemented with VIs and LAI data layers in dust-climate models, will considerably overestimate dust emission in vegetated regions with large wind speeds.Our proposed solution to this previously hidden weakness is to make prognostic ESM albedo data available to the dust module coupled to the climate model.

Overcoming Dust Emission Model Tuning to DOD
Comparing dust cycle models with DOD, indicates large errors in simulated dust magnitude and geochemical properties (Evan et al., 2014;Huneeus et al., 2011;IPCC, 2013).Consequently, dust cycle models are calibrated typically to DOD, which forces dust emissions to match dust in the atmosphere, at often unknown distances from dust sources, which hides the correct magnitude and frequency of emission events at source.We have shown here that DOD frequency poorly represents observed dust emission frequency by nearly 2 orders of magnitude, and with no spatial correlation in frequency variability.Previous studies have suggested that this inconsistency is due to the spatial bias between time of emission and downwind observation in sun-synchronous daily observations (Schepanski et al., 2012).While explaining perhaps some of the variability evident in our results, that inconsistency also illustrates the fundamental problem of calibrating dust cycle models to DOD.The inconsistency in modeled dust emission from areas unlikely to produce dust, has previously been filtered out where preferential dust source maps are used (Ginoux et al., 2012).The probability of dust emission is predefined by the topographic setting, constraining emission to drainage basins (Zender, Newman, et al., 2003).These predefined conditions limit the ability to simulate the spatiotemporal dynamics of dust emission in these areas, as well as omitting most small dust sources in other areas of the basin (Urban et al., 2018).
Using extant DPS (Hennen et al., 2022(Hennen et al., , 2023)), our results demonstrate that DOD frequency is limited to areas with highly reflective surfaces e.g., creating a bias over northern areas of the Chihuahuan Desert.The DOD frequency hotspots for the period 2001-2016 were located upwind of the DPS locations.These findings severely undermine the efficacy of dust cycle model calibration to DOD frequency, especially where dust emission occurs in relatively discrete areas surrounded by more densely vegetated areas such as in North America.Overestimation of dust emission in these environments very likely alters the magnitude, frequency, and geographical distribution of global dust emission, which currently considers continental-scale barren environments (e.g., North Africa) as persistently predominant sources of global dust (Engelstaedter et al., 2006).

Factors increasing uncertainty in dust emission modeling (problem).
Impact on dust emission modeling and proposed solution.
Poorly constrained aerodynamic roughness (z 0 ) and drag partition (R), uses R ≈ 0.91 which represents a bare soil surface, intended for E to reduce dust emission in the presence of vegetation, but E functions incorrectly.
In regions with any vegetation, R ≈ 0.91 underestimates u *ts and the inadequate adjustment by E causes overestimated dust emission, particularly in vegetated regions.
Solution: Remove E, z 0 , and z 0s and use albedo-based   * .
The incomplete form of  3 * in Q TEM (Equation 2) widely adopted in TEMs, should have u * multiplied by the drag partition R, but the correct values of are unknown (for every pixel and every time step producing dust).

Implications of Model Weaknesses for Dust Emission Modeling
This study has demonstrated that dust emission modeling can be considerably improved by utilizing a calibrated drag partition with the AEM.It contrasts with the exemplar TEM by avoiding pre-conditioning the model to bare (devoid of vegetation) z 0 and z 0s to produce R ≈ 0.91 which overestimates sediment transport before adjustment by E the bare "erodible" soil fraction.The TEMs were developed more than two decades ago when dynamic data inputs were less available.Many global dust emission studies still use static inputs, such as selective vegetation cover thresholds and bare soil fraction in global dust emission modeling (Albani et al., 2014).Preferences for which regions emit or how much vegetation to allow before dust emission ceases, have contributed to the inability to detect model weaknesses (Zender, Bian, et al., 2003).The ad hoc delineation of source regions and/or tuning of dust cycle models to DOD constrains dust emission to areas with large concentrations of dust in the atmosphere (Huneeus et al., 2011).In contrast, regional differences in magnitude and frequency of dust emission, wind speed, and particle size controlling dust residence times are at best not prioritized and at worst masked out.Furthermore, contemporary atmospheric dust loads do not enable unbiased reconstruction of past trends or to project future shifts in the location or strength of emissions (Mahowald et al., 2010).
There is also a great risk that major scientific advances made in developing dust emission schemes (Marticorena & Bergametti, 1995;Shao et al., 1996) and newly developed data/parameterizations (Prigent et al., 2012) are being under-utilized by an over-reliance on parsimonious assumptions about dust source location and erodibility to implement dust emission models.Model tuning in its various guises, particularly to DOD, makes it hard to routinely evaluate dust emission model implementation and development.Our findings suggest that it is essential to ensure that dust emission modeling is explicitly balanced between the fidelity of the dust emission scheme (processes) and the parsimony of its implementation (parameterization influenced by data availability) (Raupach & Lu, 2004).As new parameterization schemes are developed and new data sources become available, the aeolian research community will benefit from being open to critical re-evaluations and diversifications to ensure that model development balances parsimony and fidelity and avoids enduring model weaknesses.
The exemplar TEM, in common with many dust emission models, uses  3 * to calculate the magnitude of sediment transport, and predicted unreasonably large dust emission particularly in vegetated regions, because u * should be multiplied by R before being cubed and hence is overestimating   * (Webb et al., 2020).Although our exemplar TEM dust emission is adjusted by the bare soil fraction E, we have shown that this is not functioning adequately.Implementations of E using VIs and recent developments with dynamic vegetation in dust emission modeling are flawed by this weakness in E. Many of the limitations in dust emission modeling using the exemplar TEM were evident in the original dust emission scheme (Marticorena & Bergametti, 1995).However, these limitations have been ignored or overlooked to implement global dust emission schemes.The calibration of dust cycle model estimates against DOD have hidden these limitations, caused weaknesses to endure and made it difficult for the community to evaluate model developments.
Despite its multiple parameters, the exemplar TEM operates like other dust emission models explicitly controlled only by wind speed at some height U f and threshold of U ft (Ginoux et al., 2001) (e.g., GOCART).In our study, we did not include these dust emission models based on wind threshold.However, given their similarity with the exemplar TEM, our results suggest that both model types are inadequate for representing dust emission across Earth's dynamic vegetated drylands and over time.Consequently, the model weaknesses identified here most likely explain why, on monthly time scales, the relation between surface wind speed and TEMs could be linearized, and why differences between CMIP5 models appear to be due solely to wind field biases (Evan et al., 2016).Perhaps most significantly, our results explain to a large extent, how and why the use of exemplar TEM lack validity in 21st century dust emission projections (Evan et al., 2014).Large uncertainties and inter-model diversity remain in CMIP6 models and are larger than previous generations of models (Zhao et al., 2022) implying that modeled dust processes are becoming more uncertain as the latest modeling efforts continue to evaluate dust cycle models against atmospheric DOD (Klose et al., 2021).

Conclusion
Improving climate change projections requires dust models that are sensitive to, and accurately represent, dust emission responses to changing environmental conditions (wind speed, precipitation, evapotranspiration), land use and land cover dynamics.Dust cycle models typically calibrated against atmospheric DOD are assumed valid

Figure 1 .
Figure 1.Location and publication source(Baddock et al., 2011;Kandakji et al., 2020;Lee et al., 2012) inventory in New Mexico, Texas, Arizona, Colorado, Kansas, Oklahoma, and Northern Mexico between 2001 and 2016 (Kandakji), 2001-2009 (Lee) and in 2001-2009 in the Chihuahuan Desert and New Mexico (Baddock) using satellite observed dust emission point sources set against a background of average wind friction velocity normalized by wind speed (u * /U h ) derived from MODIS albedo (500 m).

Figure 2 .
Figure 2. Dust emission (kg m −2 yr −1 ) simulations shown with (a) varying soil surface wind friction velocity and (b) with varying soil surface wind friction velocity (b) normalized by wind speed at 10-m height (U 10 ) using fixed entrainment threshold u *ts = 0.2 m s −1 , clay = 10%, soil moisture function H(w) = 1 and the bare soil function E. The exemplar TEM was implemented (Equations 2 and 5) with fixed aerodynamic roughness length (z 0 ) and consequently fixed R ≈ 0.91 and fixed wind friction.The albedo-based dust emission (AEM uncalibrated) was implemented (Equations 3, 4, and 6) with varying wind friction as described in the main text, with details in Text S2 in Supporting Information S1.

Figure 4 .
Figure 4. Modeled and observed frequency at known southwestern North American satellite observed dust emission point sources (DPS), identified in satellite observations from previous studies(Baddock et al., 2011;Kandakji et al., 2020;Lee et al., 2012).For each point, the y-axis represents the observed number of DPS observations (per grid box) per year during different observation phases of the DPS data sets within the observation time period(2001)(2002)(2003)(2004)(2005)(2006)(2007)(2008)(2009)(2010)(2011)(2012)(2013)(2014)(2015)(2016).For AEM and exemplar TEM, the x-axis describes the number of modeled observations (F > 0) at DPS locations in each grid box per year during the same time period (x-axis).The x-axis describes the frequency that DOD > 0.2 per year for the same period.The least squares logarithm regression of modeled against DPS observations produced the model parameter coefficients, R 2 correlation and the square root of the mean squared difference between DPS, and model predictions (RMSE) adjusted by the degrees of freedom (df) using the number of model parameters (df = 9 for AEM; df = 12 for TEM; df = 6 for DOD frequency).

Figure 5 .
Figure 5. Mean conditions for North America during the year 2020 for peak dust season months March-May, including (a) above canopy albedo-based wind friction velocity normalized by wind speed (u * /U h ), (b) wind speed (at 10-m height), and modeled dust emission with (c; AEM cal Equation10) and without (d; exemplar TEM; Equations 2 and 5) varying aerodynamic roughness.The dust emission displayed is for wind speeds restricted to between 8.5 and 9.5 m s −1 (for comparison with Figure2).The daily maximum wind speed, described in hourly data from ERA5-Land (Source: ECMWF) are used in both models.

Figure 6 .
Figure 6.Difference maps between the year 2001 and the year 2020 for the peak dust season months March-May and only dust days (not all days), showing total difference in (a) albedo-based mean wind friction velocity normalized by wind speed (u * /U h ) and (b) wind speed (U 10 ).Minimum detectable change (MDC) in dust emission with significance (P > 0.05) with AEM cal (Equation10) varying aerodynamic roughness (c) and with exemplar TEM (Equations 2 and 5) z 0 fixed and static over time(d).Wind data are from ERA5-Land (Source: ECMWF).See Text S3 in Supporting Information for details on the calculation of the MDC.

Table 1
All physically based dust emission modeling currently assumes for simplicity an infinite supply of dry, loose erodible material is available once *ts has been exceeded.Modeled u *ts at the grain (point) scale is very likely to be much smaller than areal u *ts and incompatible with areal wind speed.Underestimated u *ts relative to u s* will overestimate the frequency of dust emission where sediment is unavailable and/or restricted by rocks and biogeochemical soil crusts.Assessment of the Priority Factors Causing Uncertainty in Dust Emission Modeling, Their Likely Impact on Dust Emission Modeling and Proposed Solutions Based on Our Research Findings Solution: Either upscale albedo-based u s* to wind speed pixels or downscale wind speed to represent aerodynamic roughness heterogeneity.