Weaknesses in Dust Emission Modelling Hidden by Tuning to Dust in the Atmosphere

Vegetation is a major control on dust emission because it extracts momentum from the wind and shelters the soil surface, protecting dry and loose material from erosion by winds. Many traditional dust emission models (TEMs) assume that the Earth’s land surface is devoid of vegetation, adjust dust emission using a vegetation cover complement, and calibrate the magnitude of modelled emissions to atmospheric dust. We compare this approach with a novel albedo-based dust emission model (AEM) which calibrates Earth’s land surface normalised shadow (1-albedo) to shelter depending on wind speed, to represent aerodynamic roughness spatio-temporal variation. Existing datasets of satellite observed dust emission from point sources (DPS) and dust optical depth (DOD) show little spatial relation and DOD frequency exceeds DPS frequency by up to two orders of magnitude. Relative to DPS frequency, both dust emission models showed strong relations, but over-estimate dust emission frequency, suitable for calibration to observed dust emission. Our results show that TEMs over-estimate large dust emission over vast vegetated areas and produce considerable false change in dust emission, relative to the AEM. It is diﬃcult to avoid the conclusion, raised by other literature, that calibrating dust cycle models to atmospheric dust has hidden for more than two decades, these TEM modelling weaknesses and its poor performance. The AEM overcomes these weaknesses and improves performance without masks or vegetation cover. Considerable potential exists for Earth System Models driven by prognostic albedo, to reveal new insights of aerosol eﬀects on, and responses to, contemporary and environmental change


Abstract (250 / 250 words)
Vegetation is a major control on dust emission because it extracts momentum from the wind and shelters the soil surface, protecting dry and loose material from erosion by winds.Many traditional dust emission models (TEMs) assume that the Earth's land surface is devoid of vegetation, adjust dust emission using a vegetation cover complement, and calibrate the magnitude of modelled emissions to atmospheric dust.We compare this approach with a novel albedo-based dust emission model (AEM) which calibrates Earth's land surface normalised shadow (1-albedo) to shelter depending on wind speed, to represent aerodynamic roughness spatio-temporal variation.Existing datasets of satellite observed dust emission from point sources (DPS) and dust optical depth (DOD) show little spatial relation and DOD frequency exceeds DPS frequency by up to two orders of magnitude.Relative to DPS frequency, both dust emission models showed strong relations, but over-estimate dust emission frequency, suitable for calibration to observed dust emission.Our results show that TEMs over-estimate large dust emission over vast vegetated areas and produce considerable false change in dust emission, relative to the AEM.It is difficult to avoid the conclusion, raised by other literature, that calibrating dust cycle models to atmospheric dust has hidden for more than two decades, these TEM modelling weaknesses and its poor performance.The AEM overcomes these weaknesses and improves performance without masks or vegetation cover.Considerable potential exists for Earth System Models driven by prognostic albedo, to reveal new insights of aerosol effects on, and responses to, contemporary and environmental change projections.
Plain Language Summary (155/200 words) Mineral dust influences Earth's systems, and understanding its impacts relies on numerical models which include large uncertainties.We compared measurements of dust optical depth (DOD) frequency and satellite observed dust emission frequency from point sources (DPS) across North America.We found up to two orders of magnitude difference between DOD and DPS frequency without statistically significant relations.Compared with DPS frequency, we found a traditional dust emission model (TEM) and an albedo-based dust emission model (AEM) over-estimated dust emission frequency by up to one order of magnitude with statistically significant relations.Relative to the AEM, TEMs incompletely formulate wind friction, over-estimate large dust emission over vast vegetated areas and produce considerable false change in dust emission.Tuning dust cycle models to dust in the atmosphere has hidden, for more than two decades, these TEM weaknesses with implications for our understanding of

Introduction
Mineral dust is central to many of Earth's systems (Shao, Wyrwoll 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 and Toon 1999).These radiative effects depend on the volume of emitted dust and the mineral composition of atmospheric plumes which differs over geographical source areas, because particle size distribution and mineralogy vary spatially and temporally (Kok, Ridley 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, Kloster et al. 2010).Comparison with atmospheric dust observations indicate that global dust cycle models include large uncertainties in simulated dust magnitude and geochemical properties (Huneeus, Schulz et al. 2011).For example, global climate models used in the Fifth Assessment Report of the IPCC failed to reproduce observed North African dust emission over the second half of the 20 th century challenging the validity of 21 st century dust-climate projections (Evan, Flamant et al. 2014).It is common for dust models to be evaluated against dust optical depth and tuned to fit the observations from North Africa (p.7809) (Huneeus, Schulz et al. 2011).However, this methodology does not enable an evaluation of the separate components of the dust cycle.Here, we are concerned that this evaluation approach has hidden weaknesses in the dust emission phase, that parameterisations which attempt to improve dust emission are being falsely rejected, and consequently the approach to model development has become biased towards parsimony rather than attempting to achieve a balance with fidelity of dust emission processes (Raupach and Lu 2004).
Dust emission models (Joussaume 1990, Marticorena and Bergametti 1995, Shao, Raupach et al. 1996) were developed more than two decades ago and their underpinning basis of sediment transport (Wolman and Miller 1960) has not changed.The magnitude of sediment transport is adjusted by the frequency of occurrence based on the wind momentum exceeding a critical sediment entrainment threshold (Wolman and Miller 1960) causing highly dynamic, nonlinear responses over space and time (Raupach and Lu 2004).When dust emission models were developed there were few continuously varying global datasets available and simplifying assumptions were made.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 (Zender, Bian et al. 2003).Varying with wind speed, drag was shown to be equivalent to shelter (Raupach 1992) and the drag partition was related to lateral cover with the caution that it only represented two dimensions (Raupach, Gillette et al. 1993).For implementa-tion in dust emission modelling, workers assumed that dryland aerodynamic roughness, including non-photosynthetic vegetation, was approximated by lateral cover from photosynthetic vegetation indices (VIs) readily available from satellite remote sensing (Evans, Ginoux et al. 2016).These later approaches using VIs also assumed that the sheltering effect of the drag was restricted only to that 'green' canopy area and that shelter did not vary with wind speed.Whilst dust emission models were being developed in these ways, it was assumed that the sediment entrainment was at the grain scale, static over time and a function of dry, loose and erodible spherical particle diameters discretised across sizes (Shao and Lu 2000) so that soil texture data (typically aggregated over depth) could be used.Similarly, the models assumed 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.
These dust emission model assumptions represent the parsimony of implementation, more than the fidelity of the dust emission processes.Furthermore, by adjusting the magnitude of the dust model estimates to dust in the atmosphere, there is no possibility of compensating for any errors in the frequency and geography of dust emission e.g., global seasonality.It is therefore difficult to have confidence that the dust emission modelling adequately represents the geographical distribution of dust emission magnitude and frequency, particularly given the last two rather critical 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 (not dust cycle model outcomes).Given the considerable advances in dust emission modelling over the last two decades, it is important to enable dust emission model outcomes which requires tackling the assumptions about sediment entrainment and sediment supply and the way in which dust emission models are evaluated.It is timely that a new dust emission frequency point source (DPS) database is available and has been used with the albedo-based dust emission model (AEM) to circumvent the assumptions about sediment entrainment and sediment supply (Hennen, Chappell et al. 2021).We follow that established approach and evaluate dust emission modelling against DPS in addition to the frequency of atmospheric dust optical depth (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 experienced modellers contributing to this study.We compare that exemplar TEM with the AEM which attempts to overcome the issues related to the soil surface wind friction velocity in the sediment transport equation, and to redress the balance towards the fidelity of process representation (Raupach and Lu 2004).

Exemplar traditional sediment transport
Vegetation attenuates dust emission by extracting momentum from the wind and sheltering a portion of the downstream soil.By reducing wind speeds ( ) at the soil surface, vegetation makes it more difficult 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.Calculation of the stream-wise 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  * (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 of sediment flux  *  (d) (m s -1 ) (Kawamura 1951).The original equation is typically rewritten in the dust modelling literature with the typographic correction and reformulated ratios (White 1979) which require a cubic term: In ESMs or reanalysis wind field models over large areas (large pixels), with horizontal resolutions that are typically on the order of >50 km, modelled wind speed at 10 m (U 10 ) is used to calculate the available above canopy  * .In recognition that vegetation exerts drag on the wind,  * must then be partitioned between the dryland roughness elements (typically vegetation), and that available for driving flux at the soil surface (  * ).The  *  is adjusted by a soil moisture function () e.g., (Fécan, Marticorena et al. 1998) 2) In the absence of being able to estimate directly   * , the  *  is divided by R for the model implementation to account for the drag partition making use of  * (Webb, Chappell et al. 2020).Following this approach, this form (Eq. 2) is incomplete because  3 * should be multiplied by R before it is cubed (Webb, Chappell et al. 2020).There is also little recognition of the uncertainty due to the inconsistency in model implementation scales.The entrainment threshold ( *  ) is calculated at the grain scale as a function of grain diameter, density and inter-particle cohesion (Shao, Raupach et al. 1996).However, the above canopy  * is for an area when measured using a wind velocity profile or when using modelled data.This approach requires  *  to be represented over the same area, which it is not.The arising substantive issues for TEMs are therefore that the incomplete form of Q TEM (Eq.2) has been widely adopted in TEMs in which large area estimates of wind speed are typically used, the correct values of R are not known (for every pixel and every time step), and  *  is not scaled correctly and therefore incompatible with areal estimates of wind speed and wind friction.
One of the common approaches to modelling dust emission in ESMs uses globally constant values of aerodynamic roughness length (z 0 ) (Woodward 2001, Tegen, Harrison et al. 2002, Zender, Bian et al. 2003).Here, we focus on the impact for large scale TEMs and our exemplar TEM uses the incomplete formulation for Q TEM (Eq.2) with fixed aerodynamic roughness length for the landscape z 0 =100 µm and the soil z 0s = 33.3µm, which fixes R(z 0 )0.91 (see Eq. 13 in Supplement), assumes that the Earth's land surface is static over time and devoid of vegetation roughness.An adjustment of the dust emission (not the sediment transport) is applied using a vegetation cover function E (described below in section 2.1.3).The emission in the ESMs is then ultimately 'tuned' down to match observed atmospheric dust (Zender, Bian et al. 2003).Here, we do not apply this global tuning.
A second, more recent approach uses satellite remote sensing to provide spatially heterogeneous estimates of z 0 only for arid and semi-arid regions, and fixed over time (Greeley, Blumberg et al. 1997, Roujean, Tanré et al. 1997, Marticorena, Chazette et al. 2004, Prigent, Tegen et al. 2005, Prigent, Jiménez et al. 2012).With this second approach it is still challenging to estimate R. We do not apply this approach to out exemplar TEM.In practice, some models use geographically preferential dust sources that limit the magnitude of dust emission (Ginoux, Chin et al. 2001, Woodward 2001, Tegen, Harrison et al. 2002, Zender, Bian et al. 2003, Mahowald, Kloster et al. 2010, Evans, Ginoux et al. 2016).In our exemplar TEM we do not use this approach to make clear in our results the cause of differing dust emission magnitude and frequency.

Albedo-based sediment transport
In our albedo-based dust emission model (AEM), the spatio-temporal variation in   * is simulated using the concept that aerodynamic sheltering of vegetation is proportional to its shadow (1-albedo) (Chappell, Van Pelt et al. 2010, Chappell andWebb 2016).This albedo-based approximation of the drag partition was investigated and tested to provide an area-weighted value, shown to be scale invariant (Chappell, Webb et al. 2018, Chappell, Webb et al. 2019, Ziegler, Webb et al. 2020).This approach enables direct calculation of   * given measurements of albedo from satellites, and enables the complete formulation for sediment flux and dust emission 3) This AEM does not require R, z 0 or z 0s and thereby has three primary parameters less than the exemplar TEM, and removes the associated sources of uncertainty.Instead, the   * is obtained directly from  ns , the normalised and rescaled areal shadow (1-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 would be ideally at freestream.Notably, this approach retains the long-established entrainment threshold  *  which at the grain-scale is inconsistent with the new area-weighted albedo-based approach.The threshold value at the grain-scale is very unlikely to represent the value over a large area (e.g., 500 m pixel) and is probably much smaller than an areal estimate.Consequently, modelled dust emission is expected to be over-estimated.However, this component of the modelling is beyond the scope of this study.

Dust emission modelling
The vertical dust mass flux (F; g m -2 s -1 ) may be calculated from Q using physically-based schemes (Shao, Raupach et al. 1996, Kok, Albani et al. 2014).One of the common approaches in regional and global applications, and that used here for the exemplar TEM and AEM, F is calculated as an empirical function of Q (Marticorena and Bergametti 1995):  =  ()()10 (0.134clay % −6.0) .(Eq. 5) The dust emission parameterisation 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.We used the latest, reliable spatially varying layer of particle size (Dai, Shangguan et al. 2019) and restricted clay % to a value of 20% consistent with previous work showing reasonable results when applied in regional models (Woodward 2001).
The proportion of emitted dust in the atmosphere M for a given particle size fraction (d) is dependent on the particle size distribution.We calculated the relative particle size surface area (Marticorena and Bergametti 1995) (M ).The vegetation cover function E was originally defined (Marticorena and Bergametti 1995) as the ratio of bare exposed surface area to total surface area when viewed from directly above (at nadir).It is used to adjust linearly the amount of dust emission by the bare soil fraction.However, sheltering is non-linear since it depends on the mutual sheltering of the roughness (typically vegetation) structure, configuration and wind speed (Chappell, Van Pelt et al. 2010).Theoretically, R in the equations above already accounts for the soil area which is exposed to the wind friction velocity relative to that sheltered by upwind roughness elements.Therefore, E is theoretically redundant in the exemplar TEM (Webb, Chappell et al. 2020).Nevertheless, it has become acceptable to assume E=1-A v where A v is the area covered by roughness elements, typically vegetation.This E is used in some ESMs so that leaf area index (LAI) or satellite 'greenness' observations e.g., vegetation indices (VIs) can be used as a surrogate of the land surface fraction occupied by green vegetation (Zender, Bian et al. 2003, Evans, Ginoux et al. 2016, Galloza, Webb et al. 2018, Sellar, Jones et al. 2019).After the sediment flux is calculated, only then is E used to adjust dust emission using the area covered by green vegetation.This planform area covered by vegetation is much smaller than the additional sheltered area downwind of the vegetation (Chappell, Van Pelt et al. 2010).In addition, E does not represent 'brown' roughness caused by dormant or dead vegetation not readily evident in VIs, or non-erodible stone covered surfaces without sediment in dryland regions where most sediment flux and dust emission occurs.Model representation without aerodynamic sheltering is crude and the use of E is unnecessary and increases parameter uncertainty; it is a prime example of emphasising parsimony in model implementation.When TEMs using LAI or VIs are applied in dust-climate ESMs the parameterization is assumed adequate for climate projections.
The albedo-based (not restricted to MODIS products) scheme for sediment flux and dust emission (AEM; Eqs. 3, 4 & 5) represents the drag partition physics without pre-tuning to a fixed land surface condition i.e., R, z 0 or z 0s and also without the need for E. Thereby, the AEM removes these additional sources of uncertainty.

Dust emission frequency point sources (DPS) and dust optical depth (DOD) frequency
Commonly, aerosol optical depth (AOD) from ground-based or large area Earth observation (EO) data are typically used to evaluate the performance and / or calibrate dust emission model simulations (Meng, Martin et al. 2021).This approach assumes that: i) dust in the atmosphere represents the dust emission process, and ii) the spatial variation in magnitude and frequency of dust emission in the 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, and on dust deposition in the dust source region, a size dependent process.To understand the extent to which AOD estimates the spatial variation in dust emission magnitude and frequency we calculated the probability of dust occurrence modelled by the dust optical depth (DOD>0.2) using the criteria established previously (Ginoux, Prospero 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.We demonstrated below (Supplement S4) that there is little impact of the chosen 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°p ixel resolution (Platnick 2015).The DOD was retrieved from those pixels in which dust emission was observed from point sources (DPS) in space and time throughout 2001-2016.All available MODIS DOD data were used.
We described in the previous section how simplifying assumptions are made in TEMs about the dynamics of vegetation sheltering.We also provided a theoretical basis for TEMs formulation to be incorrect.The correct magnitude and frequency of dust emission per unit area depends on the correct probability that sediment flux occurs causing dust emission, which itself depends on the correct  *  () (and the correct R in the exemplar TEM).However, most dust emission schemes using  *  assume that the soil is smooth and covered with an infinite supply of loose erodible material which when mobilised causes dust emission in proportion to the clay content.This (energy limited) assumption is rarely justified in dust source regions where (i) the soil is rough due to soil aggregates, rocks or gravels, (ii) sealed with biogeochemical crusts, or (iii) loose sediment is simply unavailable (Galloza, Webb et al. 2018).Here, we circumvent these assumptions following an established approach by observing dust emission frequency at dust emission sources (see Supplement S5) during satellite observations (Hennen, Chappell et al. 2021).We improve the constraints on dust emission model evaluation by calibrating the dust emission magnitude according to modelled emissions during those observed occurrences.
We define dichotomous satellite observed dust emission point source (DPS) data and its probability of occurrence P(DPS>0) as a first order approximation of the probability of sediment flux P(Q>0) leading to the proportion of dust (F) emission P(F>0) at those locations (Hennen, Chappell et al. 2021).The identification of DPS data is a highly time-consuming and labour-intensive activity.Consequently, there are few (published) studies relative to the large number of global dust source regions.Here, we collate DPS data from several previous studies in North America (Baddock, Gill et al. 2011, Lee, Baddock et al. 2012, Kandakji, Gill et al. 2020).Those studies identify the point source of dust emissions in New Mexico and Texas between 2001-2016, 2001-2009and in 2001-2009 in the Chihuahuan Desert and New Mexico, collating a single dataset of DPS data from North America.DPS observations were identified using MODIS data with visible to thermal infrared wavebands (0.4-14.4mm; Figure 1a, see Supplement S5).Modelled (AEM and exemplar TEM) and observed frequencies are aggregated by a 1°x1°grid matrix, normalizing the results to the lowest resolution data (MODIS DOD) (Fig. 1).This aggregation is performed to tackle the incompatibility of the different scales (Gotway and 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, Gill et al. 2020) due to the phase difference between timing of dust emission and availability of the imagery.The unexplained variance and incompatible scales is well-established in the geostatistical literature (Gotway and Young 2002).We reduced the unexplained spatial variance by aggregating the DPS data.For each grid box location, the observed frequency is calculated as the number of DPS observations per year during observation period (2001 -2016).The AEM and exemplar TEM modelled dust emission frequency describes F>0 at DPS locations in each grid box per year during the same period.DOD modelled frequency describes DOD>0.2 in each grid pixel per year for the same period.At the locations and across the different studies durations of those DPS data, we calculated the AEM and exemplar TEM dust emission.We compared the model estimates during DPS observed occurrence with modelled dust emission determined by the exemplar TEM and AEM.Similarly, during those same DPS observed occurrence we compared the model estimates of dust in the atmosphere approximated using DOD.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).

Large scale dust emission modelling, mapping spatial variation and change detection
We used contemporary (2001)(2002)(2003)(2004)(2005)(2006)(2007)(2008)(2009)(2010)(2011)(2012)(2013)(2014)(2015)(2016)(2017)(2018)(2019)(2020) Earth observation data including spatially and temporally varying wind speeds (at 10 m), soil moisture (0-7 cm) and soil temperature (to represent frozen ground which inhibits sediment flux) from the latest ERA5-Land data (Muñoz Sabater 2019) (hourly; 0.1°).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, Liu 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 (Eqs. 1 & 5) with R(z 0 , z 0s )0.91 fixed over space and static over time.Following the current practice, we calculated  * from the modeled 10 m wind velocity using the logarithmic layer profile theory and aeolian roughness length (Darmenova, Sokolik et al. 2009) (details are provided in the Supplement).We allowed soil moisture to vary and only in the exemplar TEM used MODIS Normalised Difference Vegetation Index (NDVI) data to calculate the bare soil fraction E. For comparison, we used the AEM (Eqs.3, 4 & 5) with soil wind friction velocity   * / ℎ 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 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, Shangguan et al. 2019) and used in both models (see Methods and Supplement S4).
All data were available from the catalogue of the Google Earth Engine (GEE) (Gorelick, Hancher 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 y -1 ).Given the availability of DPS validation data at sites in south-western USA, we restricted the mapping to North America including dust source regions bordering the USA.Testing the code and visualising the results for smaller time periods across the study area was almost instantaneous in the GEE.Data processing at 500 m and daily resolution between 2001-2020 across North America took typically less than 12 hours.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 manipulation and presentation using ArcGIS Pro.
At the sites and days when dust was observed using dust emission point sources (DPS) we compared it with the dust emission produced by the exemplar TEM, AEM and dust in the atmosphere using DOD.For the year 2020 and the main dust emission months of March-May (MAM) in North America, we analysed 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-9.5 m s -1 to emphasise the difference in the modelling approaches, which would otherwise be hidden by taking the average for all wind speeds.Finally, we also map the difference in driving variables during MAM for the year 2001 compared with the year 2020.The dust emission on dust days is similarly compared to obtain the mean difference.That mean difference is then tested for significance using the minimum detectable change (MDC) framework (Woodward 1992, Webb, Chappell et al. 2019) and the results are displayed.The minimum detectable change (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 the Supplement.

The impact of incorrect formulation and fixed drag partition (R) on dust emission modelling
We simulated dust emission separately for a smooth and rough surface with wind speed varying between 0-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(z 0 )0.91 and assumes that the Earth's land surface is devoid of vegetation roughness and static over time.With E=1, dust emission is unadjusted and increases along the upper (large dashed) curve as wind speed increases.When the land surface is partially covered in vegetation it becomes rough and E=0.5, all other conditions remain the same.In this case, dust emission increases as wind speed increases but at a consistently reduced rate (solid curve).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 .The implication is that the same amount of dust emission is produced for a range of wind speeds regardless of whether the land surface is smooth or rough (from open square to filled square).
In contrast, the albedo-based dust emission model (AEM) for the smooth 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).In a rough situation (  * / ℎ =0.022) dust emission declines along the same curve to almost zero.Despite a 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 implication is that 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 speed and roughness influences surface wind friction velocity which is essential for accurate and precise dust emission estimates.These findings are expected based on the theory described above in the Methods section.The exemplar TEM is driven by wind speed attenuated by aerodynamic roughness which is fixed over space and static over time, and dust emission is subsequently reduced by a bare soil fraction (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 as roughness changes in either space and / or time for the exemplar TEM and AEM.Since the influence of wind speed is removed on the x-axis, exemplar TEM produces no change for a given wind speed of e.g., 10 m s -1 .The cause of change in the TEM at 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), dust emission does not change except when E changes.Under a scenario with the wind speed reduced from 10 m s -1 to 8 m s -1 , the exemplar TEM F increases monotonically but at a reduced rate; that rate does not change with roughness (  * / ℎ ).Similarly, when the wind speed increases from 10 m s -1 to 12 m s -1 , the exemplar TEM F increases monotonically at an increased rate, and 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   * / ℎ (the largest increase in roughness; Figure 2b).With the greatest increase in   * / ℎ (the largest decrease in roughness) 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 merely adjusts the magnitude, not the onset of dust emission.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 modelled 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 parameterisation.

Modelled and observed dust emission frequency at DPS locations.
We reproduced DOD>0.2 probability at previously identified DPS locations across areas of southwestern North America to compare with their observed frequency (Figure 3).The probability of DOD showed little resemblance to DPS, with a distinctly different spatial pattern and considerably greater probability in some areas.Peak DOD occurred across the USA / Mexico border in the Chihuahuan Desert, while DPS peaked over the Southern High Plains in eastern New Mexico and western Texas.DOD 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 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), DOD probability 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, modelled frequency exceeds observation.Both AEM and exemplar TEM over-estimate dust emission frequency with RMSE (Log 10 ) = 0.6 and 0.76 (4 and 5.8 day per year) respectively, relative to the 1:1 line (Figure 4).Nevertheless, across all grid box data, the relation between DOD 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.8 (P<<.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.35.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 datasets within the observation time period (2001 -2016).For AEM and exemplar TEM, the x-axis describes the number of modelled observations (F>0) at DPS locations in each grid box per year during the same time period (x-axis).For DOD, the x-axis describes the frequency that DOD>0.2 per year for the same period.The least squares logarithm regression of modelled 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).

Modelling dust emission change over space and time
The mean  * / ℎ and full range of U 10 for the year 2020 are shown (Figure 5a & b).For consistency with Figure 2, the mean dust emission is shown for selected wind speeds ( ℎ = 8.5 -9.5 m s -1 ) from both AEM and exemplar TEM (Figure 5c & 5d).The spatial distribution of mean dust emission var-ied between AEM and exemplar TEM in both magnitude and spatial extent.According to AEM, large dust emissions (0.05 -0.12 kg m -2 y -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 Basin (108.5°W,42°N).These areas correspond with small  * / ℎ , and large wind speed.TEM dust emission occurred with similar magnitude but over a 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 dust emission matches closely the distribution of mean wind speed (Figure 5b) as predicted from Figure 2. The dust emission displayed is for wind speeds restricted to between 8.5-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.
Differences in mean dust emission during peak dust season (MAM) for years 2001 and 2020 greater than the minimum detectable change (MDC) significance (P<0.05), were produced for both exemplar TEM and AEM (Figure 6c and  6d).These were compared to total mean difference in  * / ℎ and U 10 during the same periods (Figure 6a and 6b).Comparing the change between the two periods, changed  * / ℎ across North America producing a range +/-0.01, with the greatest reduction (< -0.01) associated with decreased roughness in Canada, very likely caused by changes in snow coverage.Note that snow is removed from  * / ℎ 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).Further reductions in  * / ℎ (-0.01 to -0.005) occurred in discrete areas of the Southern High Plains, and northern Chihuahuan Desert.The greatest increase in  * / ℎ (> 0.01) occurred across the American Mid-West, including Minnesota, Iowa, and South Dakota.In dusty areas (Figure 5), the greatest increase in  * / ℎ (0.005 to 0.01) occurred as discrete locations within the Chihuahuan and Sonoran Desert, the Great Basin (Nevada), and the southern extent of the Southern High Plains (eastern New Mexico and western Texas).Mean changed U 10 produced 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.
Between 2001 and 2020, significant change in dust emission (F) comparing AEM and exemplar TEM varied across the range ±2 kg m -2 y -1 .The AEM produced a significant decrease in F (-1 to -2 kg m -2 y -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 showed a significant increase in F 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 was detected, the exemplar TEM produced a significant decrease of F during the 2020 period across large areas of the Great Plains (up to -2 kg m -2 y -1 ), the arid southwest (-1 to -2 kg m -2 y -1 ), including the Mojave, Sonoran, and Chihuahuan Deserts, and the Mid-West (-1 to -2 kg m -2 y -1 ).The exemplar TEM F increased significantly across the Wyoming Basin (up to 2 kg m -2 y -1 ), the Great Basin and northern Mexico (Figure 6d).

Overcoming dust emission model weaknesses using the albedo-based approach
Dust emission modelling has struggled to replicate observed dust emission magnitude and frequency, indicating an inability to adequately represent soil wind friction velocities (Evan, Flamant et al. 2014).Many of the TEMs assume homogenous bare ground, before using the complement of vegetation cover to reduce emission.Using satellite observed dust emission point sources (DPS; Figure 1) we have shown the exemplar TEM overestimates dust emission frequency by nearly 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 (Fig. 2) at correctly estimating the probability of   * 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 over-estimating dust emission frequency by 0.6 orders of magnitude (RMSE = 0.6 using log 10 ).However, the AEM is not tuned using z 0m , z 0s , R or E and the monitored normalised shadow is calibrated to wind tunnel  * / ℎ .In contrast, the exemplar TEM is tuned using selected values of z 0m and z 0s for use in R which are fixed over space and static over time and then dust emission is adjusted by E. Furthermore, most DPS used here are from predominantly barren and windy environments, with mean  * / ℎ = 0.069 and mean U 10 = 6.9 m s -1 , reducing the potential influence of dynamic vegetation.Nevertheless, the over-estimation of modelled dust emission relative to the observed frequency, occurs because of one or more of the factors described in Table 1.Those factors are classified to propose future research priorities based on the results and conclusions reached in this study, and based on our understanding of the process that have arisen during our investigation of the results.
Table 1.Assessment of the factors causing over-estimation of dust emission frequency, their likely impact on dust emission modelling and suggested priority for research investment.

Factors causing over-estimated dust emission frequency
Modelled  *  at the grain scale is very likely to be much smaller in value than that of  *  at 500 m (MODIS Dust emission modelling assumes an infinite supply of dry, loose erodible material is available once  *  has be Modelled wind speed may be too large.Scale-invariant albedo-based approach (Ziegler, Webb et al. 2020) ope The DPS are derived from polar-orbiting satellite observations, which may not accurately and completely iden The wind tunnel data used in the albedo-based drag partition calibration may not represent the complete ran Here, we use the latest version of ERA5-Land wind (at 10 m height) data at a reasonably fine (0.1°) resolution.It is evident that U 10 is over-estimated in some regions (Fan, Liu et al. 2021).However, there appears to be no systematic bias that would lead to the over-estimation of dust emission frequency.The grain scale of  *  is evidently incompatible with areal dust emission modelling and this factor appears to be the most likely cause of the over-estimated model dust emission frequency and should be a priority for future work.Without resolving the scale of  *  it is not possible to determine the impact of the assumed infinite supply of loose erodible material (Table 1).It is very likely that these two factors explain the first-order differences between the DPS frequency and the dust emission model frequency.There remains uncertainty over the use of DPS frequency.However, by comparison with dust in the atmosphere represented by DOD, the use of DPS frequency is up to two orders of magnitude smaller.There is a small, perhaps smaller-order likelihood that the original calibration of the albedo-based approach is not representative and universal, despite recent support for the approach (Ziegler, Webb et al. 2020).
Beyond the observed dust emission point sources, vegetation roughness appears influential, constraining dust emission greater than 0.1 kg m -2 y -1 to areas where  * / ℎ is no greater than 0.06, even during periods of peak (8.5 -9.5 m s -1 ) wind speed.In contrast, the exemplar TEM predicts dust emission >0.1 kg m -2 y -1 in areas where  * / ℎ is greater than 0.075, including large areas of the Great Plains.This difference is emphasized in parts of western Oklahoma (99.5°W, 35.5°N),where mean  * / ℎ > 0.08 prevent dust emission from the AEM, despite a mean U 10 > 7 m s -1 .However, in those areas TEM dust emission exceeds 0.2 kg m -2 y -1 .These contrasting estimates emphasise TEM dependency on variability in U 10 , due to the use of  3 * and the inability of R(z 0 ) fixed over space and time to correctly attenuate wind speeds by aerodynamic roughness.This limitation creates two main issues, 1) a requirement for post-process tuning, which limits the model's ability to predict dust without a priori information; 2) large scale uncertainty driven by a large spatial and temporal variability in U 10 .

Overcoming dust emission model tuning to dust in the atmosphere
Inconsistency in modelled dust emission from areas unlikely to produce dust, has previously been filtered out by utilizing a preferential dust source map (Ginoux, Prospero et al. 2012).The probability of dust emission is pre-defined by the topographic setting, constraining emission to drainage basins (Zender, Newman et al. 2003).These pre-defined conditions limit the ability to simulate the spatio-temporal dynamics of dust emission in these areas, as well as omitting most small dust sources in other areas of the basin (Urban, Goldstein et al. 2018).Furthermore, modelled dust emission frequency is typically several orders of magnitude greater than observation, creating the need for calibration when integrated into ESMs.Currently, a global observed dust emission archive does not exist, thus calibration of dust cycle models is achieved against observed dust in the atmosphere (e.g., DOD).However, we have shown that DOD poorly represents observed dust emission frequency by nearly two orders of magnitude, 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, Tegen et al. 2012).Whilst explaining perhaps some of the variability evident in our results, that inconsistency also illustrates the fundamental problem of calibrating dust emission to dust in the atmosphere.Using extant DPS (Hennen, Chappell et al. 2021), our results demonstrate that DOD is limited to areas with highly reflective surfaces e.g., creating a bias over northern areas of the Chihuahuan Desert.The DOD hotspots for the period 2001-2016 were located upwind of the DPS locations.These findings severely undermine the efficacy of dust emission model calibration to DOD, especially where dust emission occurs in relatively discrete areas surrounded by more densely vegetated areas such as in North America.Over-estimation of dust emission in these environments very likely alters the magnitude and frequency of the global dust distribution, which currently considers continental-scale barren environments (e.g., North Africa) as persistently predominant source of global dust (Engelstaedter, Tegen et al. 2006).
Our comparison of dust emission between two time periods emphasizes a previously unrealised impact of varying aerodynamic roughness in the temporal variability of dust emission magnitude.Through the calculation of dynamic   * , 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 semi-arid 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 unlikely to produce dust (e.g., temperate areas of the Great Plains and the grasslands of northern Mexico; Figure 6d).

Implications of model weaknesses for dust emission modelling
This study has demonstrated that dust emission modelling can be considerably improved by utilising a calibrated drag partition with the AEM.It contrasts with the exemplar TEM by avoiding tuning to bare (devoid of vegetation) conditions z 0 and z 0s to produce R0.91 which over-estimates sediment transport and then dust emission before adjustment by E. 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 modelling (Albani, Mahowald 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 to dust in the atmosphere, constrains dust emission to areas with large concentrations of dust in the atmosphere (Huneeus, Schulz et al. 2011).However, there may be regional differences in magnitude and frequency of dust emission, wind speed and particle size controlling dust residence times.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, Kloster et al. 2010).There is also a great risk that the major scientific advances made in developing dust emission schemes (Marticorena andBergametti 1995, Shao, Raupach et al. 1996) and newly developed data / parameterizations (Prigent, Jiménez et al. 2012) are being overlooked by an over-reliance on parsimonious assumptions about dust source location and erodibility to implement dust emission models.Model tuning in its various guises, makes it hard to routinely evaluate the dust emission implementation.We suggest that it is essential to ensure that dust emission modelling is explicitly balanced between the fidelity of the dust emission scheme (processes) and the parsimony of its implementation (parameterization) (Raupach and 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 to ensure that model development balance and avoid model weaknesses enduring.
Our 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  * should be multiplied by R before being cubed and hence is over-estimating   * (Webb, Chappell et al. 2020).Although our exemplar TEM dust emission is adjusted by the bare soil fraction E, we have shown that this is theoretically redundant and is compensating for non-dynamic R0.91.If dust modelling is concerned primarily with parsimony, our results demonstrate that in the exemplar TEM, E should be removed from the dust function (Eq.5) and should replace R in the sediment transport equation (Eq.2).This small change would make explicit RE; the partition of wind speed dependent drag is approximated by static vegetation cover excluding sheltering.It would make a considerable reduction to sediment transport and rectify the over-estimation of dust in vegetated areas.However, it will send dust emission modelling further along a path of parsimony using convenient vegetation indices, dubious in contemporary changing drylands, and unknown in climate projections.
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, Chin 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 of these model types are inadequate for representing dust emission across Earth's dynamic vegetated drylands and over time.Consequently, these 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, Flamant et al. 2016).Perhaps most significantly, our results explain to a large extent, how / why the use of exemplar TEM lack validity in 21 st century dust emission projections (Evan, Flamant et al. 2014).

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.The exemplar TEM was shown here to over-estimate dust magnitude, frequency and extent and lacks the aerodynamics in dust emission of the albedo-based approach.However, the exemplar TEM performs similarly to the AEM because the bare soil fraction E is doing the job of the static drag partition R, but in the wrong place in the equations.If dust modelling is concerned solely with parsimony, our results demonstrate that in the exemplar TEM, E should be removed from the dust function and should replace R in the sediment transport equation.This would make explicit ER and accelerate dust emission modelling along a path of using convenient contemporary vegetation indices, dubious in drylands and unknown in projections.However, we discourage that temptation, and instead encourage dust emission modelling to tip the balance back towards the fidelity of the process by embracing land surface structure and adopting an energy-based approach.The albedo-based approach can be used across timeframes and because it is areal and albedo scales linearly the approach cuts across scales.Aerodynamics can be retrieved from accurate and precise albedo from ground measurements using net radiometers, from various airborne and satellite platforms most notably MODIS, or prognostic estimates used in ESMs.The availability of prognostic albedo provides the opportunity for the albedo-based approach to be readily adopted in energy-driven Earth System Models (ESMs) more suitable for climate projections.We recognise that there is some work to be done to couple prognostic albedo between components in some ESMs.However, that work has the additional benefits of improving consistency within an ESM and reducing uncertainty and independent tuning of the components.Coupling the albedo-based approach to ESMs is expected to reduce uncertainty in dust emission and transform dust-climate change projections. *  (, ,  0 ,  0 ) =  *  ()() ( 0 , 0 ) , (Eq.7) where the entrainment threshold  *  (d) for a given size fraction d, is modified by functions including the drag partition R(z 0 , z 0s ) and the moisture content H (w). The  *  of a given d (mm): (1 + 0.006    2.5 ) 0.5 , (Eq.10) includes p a =1230 g m 3 fixed air density, p p =2650 g m 3 fixed particle density, g=9.81 m s -2 acceleration due to gravity (Marticorena and Bergametti 1995).dimensionless function H (Fécan, Marticorena et al. 1998) was developed using wind tunnel experiments to account for gravimetric surface soil moisture content w (kg 3 kg -3 ) using the difference between the potential w' based on clay content and near surface w: () = √ 1 + (1.21 ( −  ′ ) 0.68 ) (Eq. 11) where  ′ = 0.0014% 2 + 0.17% , (Eq.12) and clay is the finest fraction (expressed as a percentage) of the soil and typically less than 2 µm.
A discussion of the use of this parameterization in dust emission schemes is included elsewhere (Zender, Bian et al. 2003, Xi andSokolik 2015).We make use of the ERA5-Land volumetric soil moisture data (0-7 cm of soil layer; hourly; 0.1°).To convert from volumetric soil moisture to the required gravimetric soil moisture we divided by the soil bulk density.We assumed that the gravimetric moisture of the uppermost soil layer was 20% of the 7 cm soil layer The R(z 0 , z 0s ) is valid for small wakes (z 0 < 1 cm), and to parameterize solid obstacles only (Marticorena and Bergametti 1995).This poses a problem in applying this approach to partially vegetated surfaces such as mixed grasslands, shrublands, and other vegetation mosaics (Darmenova, Sokolik et al. 2009).
Applying different parameterizations for surfaces with similar roughness values could result in a significant discrepancy in the estimated drag partition (Darmenova, Sokolik et al. 2009).To reduce the impact of this discontinuity on R(z 0 , z 0s ), a modification is used (MacKinnon, Clow et al. 2004) because it includes a wider range of land surface types ) . (Eq.

13)
In the absence of regional and global spatio-temporal dynamics of R and aerodynamic roughness length (z 0 ) data to calculate  * from U 10 , two approaches for representing surface roughness have been developed in regional and global dust emission modelling over the last two decades.One of the older, common approaches uses globally constant values of z 0 , fixed over time (Ginoux, Chin et al. 2001, Woodward 2001, Tegen, Harrison et al. 2002, Zender, Bian et al. 2003, Mahowald, Kloster et al. 2010).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 which implies that the Earth's land surface is devoid of vegetation roughness, and static over time.With R(z 0 ) fixed, R(0) * =   * is assumed which tends to maximise dust emission.We recognize that the use of a constant value for z 0s smooths out the heterogeneity of dust sources.We also know that it is recommended to use a z 0s 1/30 of the coarse mode mass median diameter of the undisturbed soil size distribution instead of setting it to a fixed constant (which assumes that the coarse population of an undisturbed soil is equivalent to the coarse population of the soil texture (Darmenova, Sokolik et al. 2009).Nevertheless, we fixed z 0s to ensure that results were consistent with previous work (Woodward 2001, Tegen, Harrison et al. 2002, Zender, Bian et al. 2003).
A second approach is to use spatially heterogeneous estimates for arid and semiarid regions (Prigent, Jiménez et al. 2012).That work follows continued efforts to use active and passive reflectance obtained from satellite remote sensing to characterize aerodynamic roughness (Marticorena, Kardous et al. 2006).Although this approach provides an observation-based approximation of z 0 , it remains challenging to estimate R to approximate   * necessary to complete the sediment flux equation.
To implement vertical dust emission, we used the standard approach of introducing additional terms to Eq. 5 which are explained below:  =       ()  ()10 (0.134clay % −6.0) (Eq.14) Notably, no global tuning is applied to either the exemplar TEM or the AEM.Also in both models, the mass fraction of clay particles in the parent soil was allowed to vary over space but was fixed over time.We used the latest, reli-able spatially varying layer of particle size (Dai, Shangguan et al. 2019) and restricted clay % to a maximum value of 20% consistent with previous work showing reasonable results when applied in regional models (Woodward 2001).
The proportion of emitted dust in the atmosphere M for a given particle size fraction (d) is dependent on the particle size distribution.We calculated the relative particle size surface area (Marticorena and Bergametti 1995) (M ).The parameter E was defined in the main text assuming E=1-A v so that vegetation indices can be used (Shao, Raupach et al. 1996).To conform with that practice, we calculated   = −22.5 + 150   (Eq.15) from global daily NDVI from MODIS (MOD09GA Collection 6 from Terra at 500 m pixel resolution).
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 (  ).The coverage of snow in a given pixel is an areal quantity and therefore ranges between 0-1.Consequently, we applied the MODIS Normalised Difference Snow Index (Hall 2016) (MOD10A1 from Terra, daily at 500 m).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 ).We used soil temperature available in ERA5-Land and set a threshold of 273.15K above which sediment flux can occur.

S2 New parameterization of u s * by relating shelter to shadow (AEM)
To implement Eq. 1, we assume that  * is the above canopy wind friction velocity.We and its spectral influences due to e.g., soil moisture, mineralogy and soil organic carbon, were removed by normalizing (Chappell, Webb et al. 2018) with the directional reflectance viewed and illuminated at nadir (0 ∘ , ) : . (Eq. 18) This was implemented by making use of the available MODIS black sky albedo (Schaaf 2015) to estimate   , and the shadow is normalized by dividing it by the MODIS isotropic parameter f iso (MCD43A1 Collection 6, daily at 500 m) to remove the spectral influences: The f iso is a MODIS parameter that contains information on spectral composition as distinct from structural information (Chappell, Webb et al. 2018).In theory, the structural information is waveband independent (Chappell, Webb et al. 2018).The normalization of MODIS data using this parameter and that of MODIS Nadir BRDF-Adjusted Reflectance (NBAR) is sufficiently similar to remove the spectral content for all bands examined (Chappell, Webb et al. 2018).
In practice, we calculated   using MODIS band 1 (620-670 nm).This structural approach is therefore highly suited to drylands where vegetation indices perform poorly due to non-photosynthetic vegetation.
To calculate the vertical dust emission, we followed the same approach as above (Eq.14) except for E, which was not used.This new implementation provides a highly dynamic representation of the soil wind friction velocity.To this model, we applied no tuning.
Our target quantity d 2,1 is greater than zero and statistically significant and defined as (Woodward 1992):  0 ∶ ẑ ( 1 ) = ẑ ( 2 ),  1 ∶ ẑ ( 1 ) = ẑ ( 2 ) +  ( ≠ 0).(Eq.22) The alternative hypothesis H 1 is the adjustment due to  = d 2,1 which between sampling periods t 1 and t 2 is the net result of change in the property of interest during an intervening time.The uncertainty due to reaching an incorrect conclusion is the minimum detectable change (MDC) which is related to the probability of the errors on the conclusion.In general, the smaller the MDC, the larger the required sample size for a given probability of false acceptance error (De Gruijter 2006).
Our  0 ∶ d 2,1 = 0 is that the average difference in our property of interest has stayed the same over time.The alternative hypothesis H 1 : d 2,1 ≠ 0 is that the average difference in our property of interest has changed over time.In statistical hypothesis testing two types of errors may be made.We may reject H 0 and conclude that there is a positive effect when in reality there is no effect (false rejection; type-I error).We assigned a probability denoted to this type of error and decide on a value of 5% based on the implications of making a false rejection.The alternative error is that we may accept H 0 and conclude that there is no effect, when in reality there is a positive effect (false acceptance; type-II error, ).The probability 1 −  is referred to as the power of the test and is used as a quality measure with a value set at 5%.First the critical value is calculated for the mean beyond which H 0 is rejected.The power is the probability that one correctly concludes that there is a positive effect, that d 2,1 ≠ 0. The power of the test depends on d 2,1 i.e., the greater d 2,1 , the larger the power.
A two-tailed test (for change without direction) statistic is commonly based on the t-test (Woodward 1992): , (Eq.23) where X is a standard normal distribution.Re-arranging to give an expression for d 2,1 , that is the difference between means which it is possible to detect with the specified power (and size) of test or more usefully, the smallest difference detectable with at least the given power This last equation is our estimate of the difference in means and our MDC for a given set of conditions which were applied to our properties of interest.

S4 Comparison of DOD at various flagging thresholds
To determine the sensitivity of frequency of occurrence maps, MODIS DOD are observed at multiple thresholds, including DOD>0, DOD>0.1,DOD>0.2, and DOD>0.25 (Figure A1).Frequencies were observed daily at 1°resolution for areas where a minimum of one DPS observation was made (grid areas in Figure A1), consistent with the methodology in Section 3.2.These results show minimal difference in P(DOD>threshold) along the USA and Mexico border.At lower thresholds (threshold smaller than 0.2), small differences occur further north and east, where frequency of increases.Notably, where no threshold is applied (Figure A1a), dust occurrence increases in the northern areas of the figure.This type of threshold is not applied in the DOD literature because it is vulnerable to erroneous observations due to atmospheric and surface conditions that would otherwise be screened out with the application of a threshold.

S5 Dust emission point source observations (DPS) data
Satellite observed dust emission frequency point source (DPS) data include the location and moment of dust emission events from three dust emission observation studies in the semi-arid southwest of North America.For each study, satellite-derived earth observation (EO) data were acquired at regular inter- The presence of meteorological cloud or dust emission from sources upwind may prevent observation of the source of emission in a single image.The 250 m spatial resolution offered by MODIS data provides enough detail, to allow the observer to detect individual plume shapes, partially mitigating this overlapping effect (Baddock et al., 2011).During these limiting scenarios, subjective interpretation improves upon non-dynamic automated retrieval algorithms, which are required to work in all surface and atmospheric conditions (Schepanski, Tegen et al. 2012).The shape recognition and decision-making ability of human observation currently exceeds those of automated approaches, alleviating many limitations and preventing false positives observations.For each of these studies, criteria are specified for legitimate observation, including i) observation must take place during an emission event, where the deflation surface is clearly identifiable at the head of emission plume; and ii) the distinct dust source must not be obscured by either meteorological clouds or upwind dust emission plumes.As such, these data represent the cutting-edge of dust emission observations in North America, allowing spatial verification of genuine emission events.

Figure 2 .
Figure 2. Dust emission (kg m -2 y -1 ) simulations shown with varying soil wind friction velocity (a) and with varying soil wind friction velocity (b) normalised by wind speed at 10 m height ( 10 ) using fixed entrainment threshold  *  =0.2 m s -1 , clay=10%, soil moisture function H (w)=1 and the bare soil function E depending on the roughness.The exemplar TEM was implemented (Eqs. 2 & 5) with fixed aerodynamic roughness length (z 0 ) and consequently fixed R(z 0 )0.91.The albedo-based dust emission was implemented (Eq.3, 4 & 5) as described in the main text with details in the Supplement.

Figure 3 .
Figure 3.Comparison between the probability of observed dust emission point sources (DPS > 0) observations (a) and MODIS (b) dust optical depth (DOD > 0.2) during the period of DPS observation (2001-2016).All available MODIS DOD data were used.

Figure 4
Figure 4Modelled and observed frequency at known North American satellite observed dust emission point sources (DPS), identified in satellite observations(Baddock, Gill et al. 2011, Lee, Baddock et al. 2012, Kandakji, Gill et al. 2020).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 datasets within the observation time period(2001 -2016).For AEM and exemplar TEM, the x-axis describes the number of modelled observations (F>0) at DPS locations in each grid box per year during the same time period (x-axis).For DOD, the x-axis describes the frequency that DOD>0.2 per year for the same period.The least squares logarithm regression of modelled 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).

Figure 5 .
Figure 5. Mean conditions for North America during the year 2020 for peak dust season months March-May, including (a) above canopy wind friction velocity normalised by wind speed ( * / ℎ ), (b) wind speed (at 10 m height), and modelled dust emission with (c; AEM) and without (d; exemplar TEM) varying aerodynamic roughness.The dust emission displayed is for wind speeds restricted to between 8.5-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) mean wind friction velocity normalised by wind speed ( * / ℎ ) and (b) wind speed (U 10 ).Minimal detectable change in dust emission with significance (P > 0.05) with AEM varying aerodynamic roughness (c) and with exemplar TEM z 0 fixed and static over time (d).Wind data is from ERA5-Land (Source: ECMWF).See Supplement for details on the calculation of the minimum detectable change.
use a new albedo-based implementation of the sediment flux equation which avoids   * = * R and therefore does not use  * , R or the aerodynamic roughness length of vegetation (z 0 ) or that of the soil (z 0s ).Instead we used a robust direct estimation (Chappell and Webb 2016) for   *   *   = 0.0311 (exp − ns 1.131 0.016 ) + 0.007, (Eq.16) where  ns is the normalised and rescaled albedo () translated and scaled (  ) from a MODIS (500 m resolution;  min =0,  max =35) for a given illumination zenith angle (=0°) to that of the calibration data (a=0.0001 to b=0.1) using the following rescaling equation (Chappell and Webb 2016):  ns = (−)(  ()−  () max ) (  () min −  () max ) +  .(Eq. 17) When using different sources of albedo (e.g., ground-based) with different resolutions the value of  max will need to change as shown elsewhere (Ziegler, Webb et al. 2020).Shadow is the complement of albedo 1 −  dir (0 ∘ , )

Figure A2 :
Figure A2: An observed example of dust point source (DPS) identification from MODIS true colour imagery in northern Texas on February 24 2007.Source (Kandakji et al., 2020).
flectance product; NASA EOSDIS Land Processes Distributed Active Archive Center (LP DAAC), USGS/Earth Resources Observation and Science (EROS) Center, Sioux Falls, South Dakota; ISRIC SoilGrids; The work was produced whilst AC and NPW were funded by a joint grant from the National Science Foundation and Natural Environmental Research Council (EAR-1853853).