Modelled direct causes of dust emission change (2001 – 2020) in southwestern USA and implications for management

North American observed atmospheric dust has shown large variability over the last two decades, coinciding with regional patterns of vegetation and wind speed changes. Dust emission models provide the potential to explain how these direct causes of vegetation and wind speed changes are related to changing dust emission. However, those dust models which assume land cover types are homogeneous over vegetation classes and fixed over time, are unlikely to adequately represent changing aerodynamic roughness of herbaceous cover, woody cover, and litter. To overcome these model limitations and explain changing (2001 – 2020) dust emission, we used a new MODIS albedo-based dust emission model calibrated to satellite-observed magnitude and frequency of dust emission point source (DPS) data. We focused our work on four regions of southwestern USA, identified previously as the main dust emission sources. We classified the interplay of controlling factors (wind speed and aerodynamic roughness) which created disturbance regimes with dust emission change consistent with diverse land use and management drivers. Our calibrated model results show that dust emission is increasing or decreasing, in different regions, at different times, for different reasons, consistent with the absence of a secular change of observed atmospheric dust. Our work demonstrates that using this calibrated dust emission model, sensitive to changing vegetation structure and configuration and wind speeds, provides new insights to the contemporary factors controlling dust emission. With this same approach, the prospect is promising for modelling historical and future dust emission responses using prognostic albedo in Earth System Modelling.


Introduction
Measured March fine atmospheric dust concentrations have shown large variability over the past two decades at remote sites across the south-western United States (Hand et al. 2016;Aryal and Evans, 2021).The changes in dustiness across the western US have been partially explained by inter-annual and decadal-scale climate variability (Hand et al., 2016;Achakulwisut et al., 2017;Aryal and Evans, 2021).However, other contributing factors cannot be ruled out including abandonment of cropland due to declining water resources for irrigation, and disturbances including wildfires, invasive plant species, increased population and economic activity (Hand et al., 2016).It is difficult to attribute the cause of changes in atmospheric dust, not least because there is no direct relation between the main factors controlling dust emission and the different factors controlling atmospheric dust dispersion and deposition.Furthermore, land use change, land management and concurrent ecosystem changes have not been examined in detail.Lambert et al. (2020) showed that increasing atmospheric dust observations (2000-2014) occurred across the Great Plains of the USA and that increasing trend coincided with increasing cropland coverage and planting and harvesting seasons of predominant crops.Public and private rangelands of the dry western US are managed for multiple activities including livestock grazing, recreation, energy development, wildlife habitat, and to maintain the quality of air, soil, water and vegetation resources.The diverse land uses and management approaches create disturbance regimes that produce diverse ecosystem responses e.g., spread of invasive plants and changing fire regimes, and potentially accelerate or reduce (or may have no effect on) aeolian sediment transport rates and dust emission.
It remains unclear what are the direct causes of changing observed dust in the atmosphere.Potential direct causes range from dispersion of dust caused by changing wind speeds, to changing dust emission associated with vegetation composition and cover changes.Hydrological drought associated with soil moisture is implicated in the cause of change in vegetation (Lee and Gill, 2015).However, we are focused here on the direct causal effect of dust emission.Furthermore, it is important to recognise that dust in the atmosphere (e.g., dust optical depth) is not equivalent to dust emission.The amount of dust in the atmosphere is a function of the sediment transport magnitude and frequency (Wolman and Miller, 1960), and dust residency time (Textor et al., 2006), which is influenced by distance from source and dust deposition influenced by the particle size of the dust emitted and scavenged by precipitation (Mahowald et al., 2014).Importantly, if any of these factors controlling atmospheric dust changes e.g., the dust source changes and/or the particle size of dust emitted has decreased over time, this could explain the changing amounts of dust in the atmosphere.Ecosystem change affects land cover (vegetation type and organisation characterised by land surface roughness) and influences wind erosivity and indirectly the spatiotemporal distribution of dust sources.The feedback and interactions between dust emission and dust in the atmosphere therefore requires an isolated and nuanced approach to explain the direct causal effect.
North American dust emission between 2001 and 2020 was previously modelled using an albedo-based dust emission model (AEM) and calibrated using dust emission point source (DPS) frequency observations (Hennen et al., 2022;Fig. 1).These data emphasised significant spatial variability in dust emission, occurring predominantly in the biomes of the Great Plains (GP) and North American Deserts (NAD), producing 7.2 Tg y -1 combined.Within these biomes, dust emission mostly occurred in four ecoregions either side of the Rocky Mountains, including the Western Shortgrass Prairies (in GP), the Chihuahuan Desert (in NAD), the Colorado Plateau (NAD), and the Wyoming Basin (NAD).That work described the broad cause of dust emission difference over space and time.Here, we investigate change in dust emission over time for each of these four regions to determine specifically how changes in wind patterns and land cover, characterised by surface roughness, have influenced dust emissions.
Different types of disturbance may change the vegetation species composition across landscapes and the cover and structure of vegetation.Such changes can reduce the aerodynamic resistance of landscapes to erosive winds (Cowie et al., 2013) and can lead to the creation of dust emission hot spots (Gillette, 1999), or affect the erodibility of entire landscapes (Webb and McGowan, 2009).The magnitude and timing of dust emissions following ecosystem change, likely reflect the scales and rates during and/or after disturbance (Belnap et al., 2007;Munson et al., 2011).Disturbances associated with e.g., agriculture, renewable and non-renewable energy development, livestock grazing, and off-road vehicles may produce different patterns and rates of change in surface soil properties (e.g., crusting, disaggregation) and vegetation, that attenuate aeolian sediment transport dynamics (e.g., Belnap et al., 2007;Aubault et al., 2015;Flagg et al., 2014).
To date, modelling studies have focused on potential dust emission responses to broad-scale (global to regional) changes in vegetation cover but have omitted the effects of changing vegetation species composition, heights, and spatial distributions (Chappell and Webb, 2016).The effects of changing roughness structures need to be addressed to understand how ecosystem changes are influencing spatial and temporal patterns of wind erosion (Webb et al., 2014a).Field studies have identified local (plot scale) aeolian transport responses to abrupt disturbances such as vegetation removal and wildfire (Dukes et al., 2018), and vegetation state-transitions associated with woody shrub encroachment in the Chihuahuan Desert (Bergametti and Gillette, 2010;Webb et al. 2014b).However, the specific impacts of disturbance-induced ecosystem changes on North American dust emissions across landscape-to-regional scales is understudied, and unknown (Nauman et al., 2018).Our objective here is to use a dust emission model to explain the dynamics of wind and drag partition induced by changing vegetation to investigate the magnitude and frequency of spatio-temporal variation in changing dust emission across North America.Our approach is to use an established albedo-based dust emission model (AEM) calibrated to satellite observed dust emission point source frequencies to identify and explain, across different ecoregions in North America, changed dust emission timeframes and spatial patterns.

Albedo-based dust emission model (AEM)
Dust emission was estimated daily using the AEM following established approaches (Chappell and Webb, 2016;Ziegler et al., 2020;Hennen et al., 2022).Consequently, we provide only a brief description of the model which identifies the key differences and similarities with the original dust emission scheme on which it is based (Marticorena and Bergametti, 1995).Our approach similarly relies on estimates of saltation flux Q (g/m/s) to simulate dust emission F (kg m − 2 y -1 ).The Q for a given particle diameter (d), soil moisture (w), wind speed at height h (U h ), and albedo (ω) were calculated as where ρ a is air density (1.23 kg m − 3 ), g is gravitational acceleration (9.81 m s − 2 ), c is a dimensionless fitting parameter (set to 1), u *ts (d) is threshold wind friction velocity (m/s).The soil surface wind friction velocity u s* is the momentum remaining after the partition of the above canopy momentum (u * ) by roughness elements at all larger scales (topography, vegetation).The u s* is compared with the soil surface threshold friction velocity (u *ts ), which assumes a smooth, dry and loose soil surface.The u *ts is adjusted by a soil moisture function H(w) which increases the entrainment threshold as soil moisture increases (Fécan et al., 1999).Consequently, the above equation (Eq.( 1)) calculates the magnitude of sediment transport (left hand side) which is then adjusted by the frequency of occurrence (right hand side; 0 or 1) i.e., u s* > u *ts .
The main difference with the original dust emission scheme is that we estimate directly the u s* normalised by wind speed as a function of aerodynamic shelter (R a ) related to shadow: where ω ns is the normalised and rescaled shadow (ω) translated and scaled (ω n ) from a MODIS range (ω n min = 0, ω n 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): Shadow is the complement of albedo 1 − φ dir (0 • , λ) influenced by spectral characteristics e.g., soil moisture, mineralogy and soil organic carbon.This spectral information was removed by normalizing per pixel, the directional reflectance viewed and illuminated at nadir ρ(0 • , λ) (Chappell et al., 2018): To implement this approach we used MODIS black sky albedo (φ dir ) to estimate ω n and normalized shadow by dividing it by the MODIS isotropic parameter f iso (MCD43A1 Collection 6, daily at 500 m): ω n (0 The f iso is a MODIS parameter that contains information only on the spectral composition as distinct from structural information (Chappell et al., 2018).In theory, the structural information is waveband independent (Chappell et al., 2018).We have shown elsewhere that the normalization of MODIS data using this f iso parameter and that of MODIS Nadir BRDF-Adjusted Reflectance (NBAR) are able to remove the spectral content for all bands examined (Chappell et al., 2018).In practice, we calculated ω n using MODIS band 1 (620-670 nm).To retrieve the wind friction velocity as a function of U h , the daily maximum wind speed at h = 10 m above the soil surface is provided by ECMWF Climate Reanalysis, ERA5-Land hourly wind field data at 11 km spatial resolution (Muñoz Sabater, 2019).Within an ERA5 pixel, the wind speed is applied to the MODIS albedo-based u s* U h available daily at 500 m.Dust emission F (kg m − 2 s − 1 ) is calculated as: MQ(d)10 (13.4 %clay− 6.0) with 0% < clay% < 20%. (6) The % clay was allowed to vary spatially and in common with other dust emission models we restricted the max(% clay ) = 20 (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 which we calculated 1 < d < 10 µm following Zender et al. (2003) by using M = 0.87.The coverage of snow (A s ) and whether the soil surface is frozen (A f ) is used to reduce dust emission by using the MODIS normalised difference snow index (MOD10A1 Collection 6, daily at 500 m).
When dust emission models were first developed there were few continuously varying global datasets 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 et al., 2003).Under these conditions, dust emission estimates were maximised and recognised by modellers as over-estimated in the presence of vegetation.Consequently, dust emission schemes reduced estimates using E, the area of bare, exposed 'erodible' soil surface (Marticorena and Bergametti, 1995).For implementation of these approaches in dust emission modelling, researchers assumed for simplicity that dryland roughness did not vary with wind speed (was not aerodynamic) and was approximated by (at nadir) cover from photosynthetic vegetation indices (VIs and more recently derivations of the VIs including leaf area index; LAI) readily available from satellite remote sensing (e.g., Evans et al., 2016).These approaches using VIs therefore assumed that the sheltering effect of the drag was restricted only to that 'green' canopy area.This approach 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.Furthermore, this approach does not represent non-erodible stone covered surfaces without sediment also common in dryland regions.Perhaps of greatest significance for large scale dust models, this implementation of E is sensitive to vegetation reconstruction of past trends or projections of future shifts in vegetation which are crudely represented in common with critique of dust source masks (Mahowald et al., 2010).
Unlike existing dust models, the use of ω ns to dynamically estimate u s* removes the need for vegetation indices and fixed vegetation coefficients to determine effective aerodynamic roughness.Furthermore, we do not use the the ratio E. Since u s* is spatially explicit, it is not necessary to apply contemporary preferential dust source masks to precondition F (i.e., increasing F in areas perceived to have greater erodibility).
M. Hennen et al.We modelled dust emission using MODIS albedo (daily, 500 m) for the period 2001-2020 for four ecoregions in southwest North America (Hennen et al., 2022).Mean dust emission is provided on days when dust occurred (evident from the dust emission point source data).All other days are removed when making mean calculations.
Our assessment of the spatio-temporal variability and classification of dust emission sources, is with respect to the direct controlling environmental variables of only wind speed and land surface roughness which combine to form the wind friction velocity (aerodynamic nature of vegetation cover).Indirect controlling environmental variables are not considered here because they provide a typically second-order influence on the emission process and the modelling.For example, soil moisture is not considered here.Soil moisture is converted to a function which adjusts the threshold of sediment entrainment (Eq.( 1)).Soil moisture is required at the very surface of the soil which is not readily available, provided by large scale models which reduces its realistic representation.Furthermore, the process of sediment entrainment is parameterised crudely, by assuming it is static over time and fixed over space (even if soil moisture and / or vegetation are varied).The sediment is also assumed to be loose and available with an infinite supply of particles comprising segregated sizes rather than being aggregates of particles.The conventional approach to tackling these weaknesses is to calibrate the dust emission model to dust optical depth (DOD) (Hennen et al., 2022;Chappell et al., 2021).
We take a different approach and instead calibrate the dust emission model to dust emission point source (DPS) observations of dust emission frequency, subjectively identified through visual inspection of satellite data by Baddock et al. (2011), Lee et al. (2009Lee et al. ( , 2012) ) and Kukal and Irmak (2020) to calibrate the albedo-based dust emission model (Chappell and Webb, 2016).Human interpretation of satellite data has many benefits compared to automated approaches (Schepanski et al., 2007).Humans are able to identify the point at the head of a dust plume, especially where dust plumes are visually similar to bare ground (Hennen et al., 2019).Limited image sampling frequency and occasional opaque atmosphere due to clouds or atmospheric dust are unavoidable in DPS (and DOD).The DPS data may be biased towards the larger dust emission events in space and time and may not represent the smaller dust emission events.However, in the absence of extensive networks providing dust emission measurements, DPS data provide the best available and most proximal data for calibrating the location and timing of dust emission (Hennen et al., 2022).
First, we calculated the DPS probability of occurrence P(DPS > 0), a first order approximation of the probability of sediment transport P(Q > 0), which is directly proportional to the probability of dust emission P(F > 0) at those locations.Next, we equated this to study durations equal to the frequency u s* which exceeds u *ts adjusted by H(w): During each simulation, the correct response (F > 0) { 1 0 depends on the correct u *ts H(w).Importantly, the traditional dust emission models, like the AEM used here, assume that the soil surface is smooth and covered with an infinite supply of loose erodible material which when mobilised by sufficient u s* causes transport and dust emission.This (energy-limited) assumption is rarely justified in dust source regions, where the soil surface is rough due to soil aggregates, rocks, or gravels, sealed with biogeochemical crusts, or loose sediment is largely unavailable.Consequently, we follow the newly established approach (Hennen et al., 2022) and bypass those weak assumptions by using observed dust emission frequencies at DPS data locations to parameterise the entrainment threshold frequency distribution On this basis, modelled dust emission is calibrated (F cal ) following Hennen et al. (2022) using the coefficients obtained by comparing satellite observed dust emission (F) against dust emission model (Eqs.( 6) and ( 8)) estimates: where F cal is the adjustment of modelled dust emission using satellite derived dust emission point source observations (Baddock et al., 2011;Lee et al., 2009Lee et al., , 2012;;Kukal and Irmak, 2016;Kandakji et al., 2020).We recognise that these dust emission point sources provide information about observed dust emission and not emission events that are hidden by clouds.This limitation also affects satellite observations of dust in the atmosphere (e.g., dust optical depth).However, these dust emission point sources provide a valuable, large area examination of dust emission which is not currently available in other datasets.

Dust source regions
Ecoregions (Dinerstein et al., 2017) were used to describe vegetation geography, important for patterns of dust emission, wind, and aerodynamic roughness (Fig. 1).The US Environmental Protection Agency (EPA) level II ecological classification scheme divides the area into ecoregions, which are here aggregated into four vegetation-type biomes following Hennen et al. (2022).The Wyoming Basin (438) comprises mostly grasslands, with discrete areas of shrubland.The Colorado Plateau (429) and Chihuahuan Desert (428) reside mainly in the semiarid southwest, comprising mostly shrub and barren land cover types.The Western Shortgrass Prairie (402) transitions from more humid conditions in the north and east, to the semi-arid conditions in the south and southwest (Kukal and Irmak, 2016).Land cover types of this biome follow this gradient, with grasslands and croplands in the north and east, transitioning to more arid shrubland in the southwest (Fig. 1).

Analysis and classification of spatiotemporal change in dust emission
Regional statistics of dust emission conditions use daily modelled estimates from only locations (pixels) and days where dust occurred (Table 1).These data describe conditions during dust emission and do not describe the area weighted means of the entire region on all days.To directly attribute the contribution of surface roughness and wind speed (R a = u * /U h and U 10 ) relative to the change in F cal , we plot the anomalies.For each region, monthly anomalies were calculated from normalised data: where Var norm = is the normalised monthly anomaly, Var i = the monthly anomaly, and Var min and Var max represent the range of monthly values across the time series.We calculate the normalised monthly anomalies as the difference against the median (Table 1) of all normalised monthly values between 2001 and 2020.The median rather than the mean was selected to reduce skewness in the anomalies due to outliers.The relative change in F cal compared to the median is expressed absolutely (magnitude of change irrespective of ± sign) by the size of circle representing each of the monthly means.There are six classes of monthly outcomes.The first three classes (class 1-3) represent months where F cal has increased compared to the median, while the subsequent classes (class 4-6) describe months where F cal decreases relative to the long-term median.In each group, the change in each of the environmental variables are either both increasing or both decreasing (in-phase) or one is increasing while the other is decreasing (out-of-phase: class and 5).When the environmental variables are in-phase (class 1, 3, 4, and 6), the direction of their sign (increasing or decreasing) explain the limits of that system.Where both R a and U 10 are increasing (top right quadrant in Fig. 4), the increasing roughness is the limiting factor (class 1 and 4).Where both u * /U h and U 10 are decreasing (bottom left quadrant in Fig. 4), the decreasing wind speeds are the limiting factor (class 3 and 6).The seasonal average contribution (%) of each class, within each region is calculated.
Positive and negative anomalies in F cal are separated across the 1:1 line, where increasing U 10 with decreasing R a (top left triangle) produces nearly all positive anomalies (class 1-3), while decreasing U 10 with increasing R a (bottom right triangle) is predominantly negative F cal anomalies.We use a 1:1 line to define the point where change in each of the environmental variables is in-phase (both increasing or decreasing) and occurs at the same rate (normalised anomalies).During these conditions the magnitude of change in F cal (described by the size of the marker) shows the smallest rate of change, as the net change in one variable is cancelled out by the other.The farther from the 1:1 line, the rate of change in F cal increases as the variables become less in-phase.

Regional variations in dust emission and controlling environmental variables
During dust events, the Wyoming Basin (WB) produced the highest mean (0.016 kg m − 2 y -1 ) and median (0.014 kg m − 2 y -1 ) modelled dust emission (F) calibrated to observed DPS (Eq.7) of all ecoregions, with a vastly larger maximum monthly value of 0.11 kg m − 2 y -1 and a large standard deviation of 0.014 kg m − 2 y -1 (Table 1).The Chihuahuan Desert ecoregion (CD) produced the smallest F cal during dust events, with a mean of 0.006 kg m − 2 y -1 and a maximum of 0.02 kg m − 2 y -1 .The Colorado Plateau (CP) and Western Shortgrass Prairie (WSP) produced similar F cal , with a mean between 0.009 and 0.01 kg m − 2 y -1 and a maximum monthly average of between 0.02 kg m − 2 y -1 and 0.03 kg m − 2 y -1 .Mean surface roughness conditions during dust events were approximately the same across each of the regions (R a = u * /U h = 0.06 -0.07), with WSP producing the largest (roughest) R a = 0.07, while WB had the largest variance and CD had the smallest.Regional mean surface winds (U 10 ) were largest in WSP (9 m s − 1 ), followed by WB (8.8 m s − 1 ), CP (8.4 m s − 1 ), and smallest in CD (8.2 m s − 1 ).WSP had the largest mean surface winds (U 10 = 9 m s − 1 ), and the CD had the smallest (8.1 m s − 1 ).Variability in U 10 changes per region, with WB producing the largest standard deviation (0.8 m s − 1 ) and maximum monthly value (11.95 m s − 1 ).By comparison, variability in all other regions was far less (~0.5 m s − 1 ), with maximum U 10 = 9.4 -10.6 m s − 1 .The long-term (2001-2020) median was used to calculate the anomalies in the subsequent results.

Temporal variation in dust emission and controlling environmental variables
Seasonal trends in modelled dust emission (F cal on dust days), surface roughness (R a = u * /U h ), and surface winds (U 10 ) of the four ecoregions are described by monthly area averages in Fig. 2. For WB, CD, and WSP, a springtime (MAM) peak in F cal (0.01-0.02 kg m − 2 y -1 ) was estimated.These peak periods of F cal correspond to seasonal patterns in U 10 , with maximum winds of 8.8-9.8 m s − 1 during April.Modelled dust emissions from the CP region peaked (0.025 kg m − 2 y -1 ) during the winter months (DJF), corresponding with monthly mean winds of 8.2-8.7 m s − 1 and monthly minimum vegetation roughness (R a = 0.05).All regions produced decreasing F cal (0.001-0.01 kg m − 2 y -1 ) during the summer months (JJA) as U 10 decreased to a minimum and despite decreasing R a .
Daily mean F cal varied by ±0.0025 kg m − 2 y -1 between years, with a general increase in the least squares linear regression gradient (trend) for WSP, CD, and WB between 2001 and 2010 (significantly in WSP), before decreasing through the subsequent decade.Each of these regions produce a significant increase in U 10 during the 2001-2010 period.The CD and WB regions both had a significant increasing trend in R a during the same period, while the WSP had only a small significant change in R a .However, we found no significant trend in U 10 or R a for each of these regions in the period 2011-2020.Both U 10 or R a decreased initially until around 2012-2015, then increased until 2017 before varying about the mean between 2017 and 2020.The CP region had large variability in F cal , with multiple peaks, including 2002/03, 2010, 2012, 2014, and 2018 of similar magnitude to the maximum positive anomalies modelled in other regions.The F cal estimates for CP are distinct from the other regions and show no significant trends in R a or U 10 .We explore the differences between regions and the peculiarities of CP in the Discussion.

Spatiotemporal variability in dust emission and direct controlling environmental variables
The spatial variability in the mean and changing (regression gradient) F cal , U 10 , and R a = u * /U h are shown per pixel in Fig. 3. Mean dust emission during all dust days was largest (>0.06 kg m − 2 y -1 ) in discrete locations within the CP (429 in Fig. 3) and WB ecoregions (438 in Fig. 3).Large dust emission (>0.04 kg m − 2 y -1 ) was predicted throughout southern areas of CP, along western parts of the WSP (402 in Fig. 3) and northern CD (428 in Fig. 3).Small F cal (<0.0001 kg m − 2 y -1 ) was predicted across large parts of the CD predominantly in the south, northern CP and northeast WSP.U 10 are quite consistent across the four ecoregions (average 9.4-10.5 m s − 1 ).The largest winds occurred in areas within WSP (>12 m/s), and decrease towards the west, with a minimum of 8-9 m/s over central regions of CP.Vegetation roughness was smallest over the CP, northern and eastern areas of CD and the WB.In contrast.WSP had increasing u * /U h values from 0.065 in the west and 0.08 in the east.
The gradient of F cal during the two periods shows contrasting trends, with large areas increasing during 2001-2010 and a greater area reducing during the subsequent 2011-2020 period.The largest variations (gradient > 0.1) were predicted in large parts of CP, WSP, and northern areas of CD.Variations in F cal across CP were associated mostly with changes in U 10 , as R a changes infrequently during the two periods.Here, U 10 alternated between an increasing trend (2001)(2002)(2003)(2004)(2005)(2006)(2007)(2008)(2009)(2010) to a reducing trend during 2011-2020 in areas where modelled F ca varied in-kind.In the CD, changes in F cal appear most influenced by changes in both U 10 , and R a , increasing during 2001-2010 mostly in the north, where R a remained the same but U 10 increased.F cal reduced in the same region during 2011-2020, as R a increased, and U 10 remained constant.During 2001-2010, F cal increased in the WSP predominantly in the east where mean U 10 was greatest.This increase in F cal coincided with a reduction in R a , while changes in U 10 varied throughout the same region.Trends in F cal across WB were generally smaller than in other regions (gradient (<0.1), with alternating patterns from east to west.In the east, changes in F cal followed changes in U 10 , increasing initially (2001-2010) before decreasing in the period 2011-2022, while R a remained constant.In the northwest, small changes in F cal occurred concurrently with opposite change in R a , with F cal increasing between 2001 and 2010 (R a decreases) and decreasing between 2011 and 2020 (R a increases).

Classifying changes in dust emission according to environmental conditions
Across all ecoregions, positive and negative anomalies in F cal are separated across the 1:1 line, where increasing U 10 with decreasing R a (top left triangle) produces nearly all positive anomalies (class 1-3), while decreasing U 10 with increasing R a (bottom right triangle) is predominantly negative F cal anomalies (Fig. 4).This coherence suggests a uniform response across the region, where the dominating environmental conditions (U 10 and R a ) in each month dictate the response in F cal .The WB region provides a notable exception, with multiple months with decreasing F cal anomalies occurring above the 1:1 line.This occurs predominantly during wind limited conditions (class 3 and 6) between DJF and MAM (Table 2).In this case, months with similar environmental conditions in U 10 and R a create different patterns of F cal , suggesting a lack of spatial coherence within the ecoregion, where the mean conditions do not represent discrete areas of large dust emission change.
In each of the ecoregions, the largest rate of change typically occurs in Class 1 and 2, where positive anomalies are driven by an increase in U 10 and either decreasing roughness (Class 2) or moderately increasing roughness (Class 1).Across all regions, Class 1 (increasing dust emission, roughness limited) and Class 6 (decreasing dust, wind limited) make up the largest proportion (51 % combined) of monthly conditions (Table 2).Class 3 (increasing dust, wind limited) and Class 4 (decreasing dust, roughness limited) conditions are least frequent (20.2 % combined).
Opposite phase conditions (increasing F cal -Class 2 and decreasing F cal -Class 4) each comprise 14.3 % of the monthly conditions.Regionally, this pattern is replicated in WSP, CP, and CD, while WB has the largest proportion of months during the opposite phase (Class 2 and 4) conditions (Fig. 3).
Seasonally, each region shows a shift in the proportion of each class representing the change in conditions throughout the year (Table 2).In both CD and CP, predominant dust emission controls change between Class1 and Class 6 throughout the year.
In CD, Class 1 conditions dominate DJF (33.3 %) and MAM (78.3 %), as increasing modelled F cal coincides with increasing U 10 , despite increasing R a .Decreasing modelled F cal (Class 6) typically occurs in JJA (45 %) and SON (76.7 %), despite decreasing R a as U 10 decreases faster.This same pattern occurs on the CP, albeit with an increase in modelled  In WSP and CD, increasing F cal during DJF -MAM and decreasing F cal during JJA -DJF are evident.However, the initial increase in F cal during DJF, occurs despite reducing U 10 as this is offset by a greater reduction in R a (i.e., Class 3-41.7 %), before U 10 and R a both increase during MAM (Class 1 -53.3 %).During JJA, modelled F cal typically reduced (61 % of the time), mostly as U 10 decreases (45 % Class 5 and 6 combined) or R a increases faster than increasing U 10 (Class 4 -26.7 %).Decreasing F cal continues in SON (98 % of the time), with decreasing U 10 76.7 % of the time (Class 5 and 6) and increasing R a 58.4 % of the time (Class 4 and 5).
The WB has increasing F cal during DJF with R a decreasing 61 % of the time and U 10 increasing 23.3 % of the time.During MAM, predominantly increasing F cal coincides with an increase in U 10 58.3 % (Class 1 and 2) and reducing R a (Class 2 and 3) 51.7 % of the time.During JJA, F cal predominantly decreases (60 %), yet Class 1 produces the largest frequency of a single condition (40 %), increasing F cal as U 10 increases at a faster rate relative to increasing R a .

Anomaly maps showing change over time
The direct causal change in dust emission is restricted to wind speed and land surface roughness (aerodynamic nature of vegetation cover).The cause of dust emission is different in different classes as explained in the text related to Fig. 3.Here those explanations are placed into their spatial context which is notably spatially coherent (Fig. 5).The main results are divided in two.The first are those with increasing dust emission and increased wind speeds: The magnitude of the dust emission anomalies is shown using the symbol size and the sign of the dust emission anomalies is shown as a positive above the (dashed, 1:1) line and a negative below the line.The seasonal average contribution (%) of each class, within each region is displayed (see Table 2).
• In a few large patches in all regions, but mainly in the Wyoming Basin and the Colorado Plateau (Class 3; yellow), modelled dust emission increased despite wind speed decreasing, because the decrease in roughness produced a relative wind speed increase.• In northern Chihuahuan Desert, southern Colorado Plateau and the Wyoming Basin (Class 2; red), modelled dust emission increased because wind speed has increased, and roughness has decreased.• In other parts of the Western Shortgrass Prairie and the Chihuahuan Desert (Class 1; orange), modelled dust emission increased despite increasing roughness because the increased wind speed was sufficient to overcome the change.
The main results associated with decreasing dust emission and increasing roughness are: • Over large areas in all regions, except the Western Shortgrass Prairie (Class 6; cyan), modelled dust emission decreased despite decreased roughness because winds decreased which produced a relative roughness increase.• In patches throughout all regions (Class 5; blue), modelled dust emission decreased because wind speed has decreased and roughness increased.• In mainly the Western Shortgrass Prairie but also in small areas of other regions (Class 4; purple), modelled dust emission has decreased despite increased wind speeds because the increased roughness is sufficient to overcome the change.

Discussion
Regional statistics of calibrated, modelled dust emission (Table 1) showed that on dust days the Wyoming Basin (WB) was predicted to have produced the largest mean amount of dust (1.84 kg m − 2 y -1 ) of all ecoregions, the largest maximum monthly average (15.75 kg m − 2 y -1 ) and the greatest variability.The WB also had the greatest variability in surface roughness and wind speeds.All other regions had smaller and similar variability.This indicates that the WB ecosystem behaviour and dust response is very different to the other ecoregions.Seasonality of the ecoregions in calibrated dust emission, surface roughness and surface winds are described by monthly area averages (Fig. 2).All regions except the Colorado Plateau (CP) showed a springtime peak in F cal (0.01 -0.02 kg m − 2 y -1 ) corresponding to seasonal patterns in U 10 .The CP region modelled peak dust emission occurred during the winter months (DJF), corresponding with large mean winds and minimum vegetation roughness (snow was eradicated in the analysis).During the summer months, modelled dust emission decreased in all regions when wind speed was at a minimum regardless of roughness.These basic results indicate that patterns and trends in dust emission are the net outcome of spatiotemporal variation which requires a nuanced approach to exploring and explaining the cause of change.
We found modelled dust emissions changed (2001-2020) within all ecoregions (Fig. 1).Increased dust months (mainly centred on New Mexico) were associated with either directly increasing wind speeds or relative wind speed increases (Figs. 3 and 5).Decreased dust emission was predicted to occur in all ecoregions and all States, with notable predominance in Mexico, and was due directly to increased roughness or indirectly to the relative roughness increase in roughness.It is straightforward to understand the physical basis represented in the model, an increase in roughness reduces the wind friction velocity, which reduces the wind momentum flux at the soil surface and reduces sediment transport and dust emission (class 5) and vice versa (class 2).We described these cases as out-of-phase.Although well described in the literature as aerodynamic roughness, the in-phase cases can produce an increase in either wind speed or roughness.For example, cases 1 and 4 both have increasing roughness and increasing wind speed but can produce either an increase or a decrease in dust emission depending on the balance of the factors.Similarly, cases 3 and 6 both have decreasing roughness and decreasing wind speed but can produce either an increase or a decrease in dust emission depending on the balance of the factors.These cases demonstrate an equifinality in sediment transport and dust emission depending on the influence of the relative wind speed or roughness embodied in the wind friction velocity.These new insights and clarity of explanation are enabled by our recent improvements in dust emission modelling based on albedo (Chappell et al., 2010;Chappell and Webb, 2016;Hennen et al., 2022) and specifically the ability to quantify the soil surface wind friction velocity in sediment transport and dust emission models.
Our albedo-based dust emission model (AEM) results calibrated to satellite-observed dust emission point source (DPS) data describe increasing and decreasing dust emission in different regions across southwest USA and changing dust emission over time with no distinct trend.As explained in the previous paragraph, these results are the consequence of the controlling factors, wind speed and roughness, being in-phase or out-of-phase with each other in different regions and at different times.The complex interplay of controlling factors explained here is consistent with the diverse land use and management drivers, which create disturbance regimes producing diverse ecosystem responses and potentially accelerate or reduce (or may have no effect on) aeolian sediment transport rates and dust emission.Although this complex interplay of controlling factors produces an equifinality of outcomes, dust emission is increasing and decreasing in different regions at different times..Given the related but fundamental differences between dust emission occurrences and transported atmospheric dust (attenuated by dispersion, grain-size dependent processes related to deposition, and wind speed), it is reasonable that our work found no detectable temporal trend in dust emission.Furthermore, our work found considerable spatial coherence in the changing dust emission.Our results (Table 2) showed that, among seasons, most increased dust

Table 2
The proportion of months in each season and each segment of the anomaly time series (2001-2020;Fig. 3)  emission occurred mainly in March-May in all regions, except the Wyoming Basin when wind speed and roughness were both increasing (in-phase wind speed increase).However, most decreased dust emission occurred mainly in September-November when the wind speed and roughness were both decreasing (in-phase relative roughness increase).
It is important to note that our analysis is dependent on the data used to drive the modelling.The ERA5-Land wind data are provided at 11 km pixel resolution.These data are downscaled from 40 m to 10 m height based on aerodynamic roughness (z 0 ) which is static over time and classified over space.Consequently, these winds are currently incapable of representing processes below that grid scale.For example, these modelled winds are unlikely to represent thunderstorm outflow events which are important for North American summer time dust emission in the southwest.The quality of the wind data may be important for outcomes associated with the Wyoming Basin.There is considerable scope to improve this downscaling of ERA5-Land using albedo-based estimates of z 0 .
The new AEM's capability to adequately represent surface wind friction velocity (and aerodynamic roughness), enables explanations of dust emission which include roughness.Without including roughness and, therefore, the complete explanation of the wind friction velocity, it is difficult to avoid interpretations of responses which privilege wind speed.Furthermore, roughness is the key to management since change in vegetation is the only tool available for managers to reduce wind speed and to protect the soil surface from sediment transport and dust emission (Mahowald et al., 2007).Our results demonstrate the importance of both timing of land management activities that influence surface roughness and the amount of roughness change relative to seasonal wind speed changes for increasing or effectively reducing dust emissions.A practical implication of our findings is that management (e.g., restoration) objectives and land health assessments that use ground cover benchmarks to assess wind erosion must also account for effects of wind speed (Webb et al., 2020).Maintaining ground cover (surface roughness) above a target or benchmark level may be effective for avoiding wind erosion in one season but not another with stronger winds.This suggests that a conservative approach to setting ground cover targets for wind erosion and dust emission control is needed to avoid under-protecting soil resources and air quality.

Conclusions
The net outcome of spatiotemporal variation in the direct causal factors of wind speed and aerodynamic roughness are varying patterns and trends in dust emission across persistent dust source areas in the US.Although soil moisture is a well-established control on the entrainment threshold our model calibration circumvents its influence and it cannot be separated.Our per pixel analysis explained dust emission responses for the first time as 'in-phase' and 'out-of-phase' wind speed and surface roughness (vegetation) changes, demonstrating a sophisticated interplay between the controlling factors that produce equifinality in dust emission change.These findings are consistent with diverse land use and management drivers changing over space and time, which create diverse disturbance regimes in ecosystem responses.Dust emission is increasing and decreasing in different regions at different times for different reasons consistent with existing interpretations of climate variability used to explain changes in observed dust in the atmosphere.This new generation of calibrated AEM, sensitive to changing vegetation structure and configuration, provides new insights to understand the factors controlling dust emission.Management can only directly control roughness.With these new insights, land managers can mitigate soil erosion effectively by identifying change in dust emission associated with changes in wind and roughness.

Fig. 1 .
Fig. 1.The main dust emitting ecoregions of North America and their land cover types.Land cover data (left) are described by MODIS product MCD12Q1.The mean calibrated dust emission (F cal ; kg m − 2 y -1; right) is described by the MODIS albedo-based dust emission model following the same approach as Hennen et al. (2022).The ecoregions (Dinerstein et al., 2017) include 402 = Western Shortgrass Prairie, 428 = Chihuahuan Desert, 429 = Colorado Plateau and 438 = Wyoming Basin shrub steppe.

Fig. 2 .
Fig. 2. Seasonality (column 1; a, c and e) and annual anomalies from the long-term (2001-2020) median (column 2; b, d and f) of modelled data of the four contributing ecoregions of North America.Dust (F cal ; kg m − 2 y -1 ; <10 µm) calibrated to observed dust frequency during the period 2001-2020 (top row; a and b), wind friction velocity normalised by wind speed (R a = u * /U h ; middle row; c and d) which describes the land surface roughness, and wind speed (U 10 ; bottom row; e and f).The ecoregions are denoted by CD = Chihuahuan Desert (429), CP = Colorado Plateau (438), WSP = Western Shortgrass Prairie (429) and WB = Wyoming Basin (428).

Fig. 4 .
Fig. 4. Monthly anomalies from the long-term (2001-2020) average of wind speed (U 10 ) and of wind friction velocity normalised by wind speed (R a = u * /U h ) which describes the land surface roughness, and dust emission.The anomalies are shown for the four regions: Western shortgrass prairie (a); Colorado Plateau shrublands (b); Wyoming Basin shrub steppe (c); Chihuahuan desert (d).The magnitude of the dust emission anomalies is shown using the symbol size and the sign of the dust emission anomalies is shown as a positive above the (dashed, 1:1) line and a negative below the line.The seasonal average contribution (%) of each class, within each region is displayed (see Table2).

Fig. 5 .
Fig. 5. Map (a) showing classes of annual average dust emission anomalies increase (classes 3, 2 and 1; yellow, red and orange;) and decrease (classes 6, 5 and 4; cyan, blue and purple) relative to the long-term (2001-2020) mean dust emission anomalies (b; F cal anomalies, kg m − 2 y -1 ) caused by change in wind friction velocity normalised by wind speed (c; R a = u * /U h anomalies) and wind speed (d; U h anomalies, m/s) relative to their respective long-term means.(For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.) including Western Shortgrass Prairie (WSP; 402); Wyoming Basin (WB; 438); Colorado Plateau (CP; 428); Chihuahuan Desert (CD; 429).F cal values are calibrated against observed DPS frequency (equation (7)).