Retrieving the global distribution of threshold of wind erosion from satellite data and implementing it into the GFDL AM4.0/LM4.0 model

Abstract. Dust emission is initiated when surface wind velocities exceed the threshold of wind erosion. Most dust models used constant threshold values globally. Here we use satellite products to characterize the frequency of dust events and surface properties. By matching this frequency derived from Moderate Resolution Imaging Spectroradiometer (MODIS) Deep Blue aerosol products with surface winds, we are able to retrieve a climatological monthly global distribution of wind erosion threshold (Vthreshold) over dry and sparsely-vegetated surface. This monthly two-dimensional threshold velocity is then implemented into the Geophysical Fluid Dynamics Laboratory coupled land-atmosphere model (AM4.0/LM4.0). It is found that the climatology of dust optical depth (DOD) and total aerosol optical depth, surface PM10 dust concentrations, and seasonal cycle of DOD are better captured over the dust belt (i.e. North Africa and the Middle East) by simulations with the new wind erosion threshold than those using the default globally constant threshold. The most significant improvement is the frequency distribution of dust events, which is generally ignored in model evaluation. By using monthly rather than annual mean Vthreshold, all comparisons with observations are further improved. The monthly global threshold of wind erosion can be retrieved under different spatial resolutions to match the resolution of dust models and thus can help improve the simulations of dust climatology and seasonal cycle as well as dust forecasting.



Introduction 24
Mineral dust is one of the most abundant aerosols by mass and plays an important 25 role in the climate system. Dust particles absorb and scatter solar and terrestrial radiation, 26 thus modifying local energy budget and consequently atmospheric circulation patterns. 27 Studies have shown that the radiative effect of dust can affect a wide range of 28 environmental processes. Dust is shown to modulate West African (e.g., Miller and 29 Tegen, 1998;Miller et al., 2004;Mahowald et al., 2010;Strong et al., 2015) and Indian 30 (e.g., Jin et al., 2014;Vinoj et al., 2014;Jin et al., 2015;Jin et al., 2016;Solmon et al., 31 2015;Kim et al., 2016;Sharma and Miller, 2017) monsoonal precipitation. During severe 32 droughts in North America, there is a positive feedback between dust and the 33 hydrological cycle (Cook et al., 2008(Cook et al., , 2009. African dust is also found to affect 34 Atlantic tropical cyclone activities (e.g., Dunion and Velden, 2004;Wong and Dessler, 35 2005;Evan et al., 2006;Strong et al., 2018). When deposited on snow and ice, dust 36 reduces the surface reflectivity, enhancing net radiant energy loading and accelerating 37 snow and ice melting, and consequently affecting runoff (e.g., Painter et al., 2010;38 Dumont et al., 2014). Dust can serve as ice nuclei and affect the formation, lifetime, and 39 characteristic of clouds (e.g., Levin et al., 1996;Rosenfield et al., 1997;Wurzler et al., 40 2000;Nakajima et al., 2001;Bangert et al., 2012), perturbing the hydrological cycle. Iron 41 and phosphorus enriched dust is also an important nutrient for the marine and terrestrial 42 ecosystems and thus interacts with the ocean and land biogeochemical cycles (e.g., Fung 43 et al., 2000;Jickells et al., 2005;Shao et al., 2011;Bristow et al., 2010;Yu et al., 2015). 44 Given the importance of mineral dust, many climate models incorporate dust 45 emission schemes to simulate the life cycle of dust (e.g., Donner et al., 2011;Collins et 46 al., 2011;Watanabe et al., 2011;Bentsen et al., 2013). Mineral dust particles are lifted 47 from dry and bare soils into the atmosphere by saltation and sandblasting. This process is 48 initiated when surface winds reach a threshold velocity of wind erosion. The value of this 49 wind erosion threshold depends on soil and surface characteristics, including soil 50 moisture, soil texture and particle size, and presence of pebbles, rocks, and vegetation 51 residue (e.g., Gillette et al., 1980;Gillette and Passi, 1988;Raupach et al., 1993;Fécan et 52 al., 1999;Zender et al., 2003;Mahowald et al., 2005), and thus varies spatially and 53 temporally (Helgren and Prospero, 1987). Due to a lack of in-situ data at global scale and 54 uncertainties on these dependencies, most dust and climate models prescribe a spatially 55 and temporally constant threshold of wind erosion for simplicity. Globally uniform 56 values (e.g., around 6 to 6.5 m s -1 ) are either directly used over dry surfaces (e.g., Tegen 57 and Fung, 1994;Takemura et al., 2000;Uno et al., 2001;Donner et al., 2011) or with 58 modulations related to other factors, such as soil moisture (e.g., Takemura et al., 2000;59 Ginoux et al. 2001;Zender et al., 2003;Kok et al., 2014a). For instance, in the 60 Geophysical Fluid Dynamics Laboratory coupled land-atmosphere model AM4.0/LM4.0 61 (Zhao et al., 2018a, b), a constant threshold of 6 m s -1 is used. 62 The threshold of wind erosion may be approximately inferred using observations. 63 For instance, Chomette et al. (1999) used the Infrared Difference Dust Index (IDDI) and 64 10 m winds reanalysis from the European Centre for Medium-Range Weather Forecasts 65 (ECMWF) between 1990 and 1992 to calculate the threshold of wind erosion over seven 66 sites over the Sahel and Sahara. The IDDI was used to determine whether there was a 67 dust event for subsequently calculating an emission index defined as the number of dust 68 events to the total number of potential events. The distribution of surface wind speed was 69 Atmos. Chem. Phys. Discuss., https://doi.org/10.5194/acp-2019-223 Manuscript under review for journal Atmos. Chem. Phys. Discussion started: 19 March 2019 c Author(s) 2019. CC BY 4.0 License.
Integrated Trajectory (HYSPLIT) model (Draxier and Hess, 1998) to forecast dust 93 surface concentrations. It was found that major observed dust plume events in June and 94 July 2007 were successfully captured by the model. Later, Ginoux and Deroubaix (2017)  For individual dust events, the threshold of friction velocity can also be 98 determined by fitting a second-order Taylor series to dust saltation flux measurements 99 (Barchyn and Hugenholtz, 2011;Kok et al., 2014b). 100 Nonetheless, a global distribution of threshold of wind erosion based on 101 observation that may be implemented in climate models is still lacking. In this study, we 102 propose a method to retrieve monthly global threshold of wind erosion (hereafter, 103 V threshold ) for dry and sparsely-vegetated surface using high-resolution satellite products 104 and reanalysis datasets. This two-dimensional threshold is then implemented into the 105 Geophysical Fluid Dynamics Laboratory (GFDL) coupled land-atmosphere model, 106 AM4.0/LM4.0 (Zhao et al., 2018a, b). The benefits of using this threshold in simulating 107 present-day climatology and seasonal cycles of dust are analyzed by comparing the 108 model results with observations. 109 The data and method used to retrieve the threshold of wind erosion are detailed in 110 section 2. The distribution of the derived V threshold and its implication in the climate model 111 is presented in section 3. Section 4 discusses the uncertainties associated with this 112 method, and major conclusions are summarized in section 5. 113 114 115 Atmos. Chem. Phys. Discuss., https://doi.org/10.5194/acp-2019- Sayer et al., 2013): aerosol optical depth (AOD), single-scattering albedo, and the 122 Ångström exponent. All the daily variables are first interpolated to a 0.1° by 0.1° grid 123 using the algorithm described by Ginoux et al. (2010). We require that the single-124 scattering albedo at 470 nm to be less than 1 for dust due to its absorption of solar 125 radiation. This separates dust from scattering aerosols, such as sea salt. Then a continuous 126 function relating the Ångström exponent, which is highly sensitive to particle size (Eck et 127 al., 1999), to fine-mode AOD established by Anderson et al. (2005;their Eq. 5) is used to 128 separate dust from fine particles. Details about the retrieval process and estimated errors 129 are summarized by Pu and Ginoux (2018b). High-resolution MODIS DOD products have 130 been used to identify and characterize dust sources (Ginoux et al., 2012;Baddock et al., 131 2016) and examine the variations of dustiness in different regions (e.g., Pu and Ginoux, 132 2016, 2018b. 133 Following the recommendation from Baddock et al. (2016), who found the 134 dust sources are better detected using DOD with a low-quality flag (i.e., QA=1) than that 135 with a high-quality flag (i.e., QA=3) as retrieved aerosol products were poorly flagged 136 over dust source regions, we also use DOD with the flag of QA=1. Both daily DOD 137 retrieved from Aqua and Terra platforms are used by averaging the two when both 138 Atmos. Chem. Phys. Discuss., https://doi.org/10.5194/acp-2019- products are available or using either one when only one product is available. Since 139 Terra passes the equator from north to south around 10:30 am local time (LT) and Aqua 140 passes the Equator from south to north around 13:30 pm LT, an average of the two 141 combines the information from both morning and afternoon hours. This process also 142 largely reduces missing data (Pu and Ginoux, 2018b Soil moisture is an important factor that affects dust emission (Fécan et al., 1999). Scanning Radiometer 2 (AMSR2) sensor onboard the JAXA GCOM-W1 satellite (from 153 July 2012 to June 2017) from the University of Montana (Du et al., 2017a;Du et al., 154 2017b) was used. Both AMSR-E and AMSR2 sensors provide global measurements of 155 polarized microwave emissions at six channels, with ascending and descending orbits 156 crossing the equator at around 1:30 pm and 1:30 am LT, respectively. The VSM 157 retrievals are derived from an iterative retrieval algorithm that exploits the variable 158 sensitivity of different microwave frequencies and polarizations, and minimizes the 159 potential influence of atmosphere, vegetation, and surface water cover on the soil signal. 160 The VSM record represents surface (top ~2 cm) soil conditions and shows favorable 161 Atmos. Chem. Phys. Discuss., https://doi.org/10.5194/acp-2019- global accuracy and consistent performance (Du et al. 2017b), particularly over areas 162 with low to moderate vegetation cover that are also more susceptible to wind erosion. 163 The horizontal resolution of the product is about 25 km by 25 km, and the daily product 164 from January 2003 to December 2015 is used. The ascending and descending obit VSM 165 retrievals are averaged to get the mean VSM for each day. 166 167

3) Snow cover 168
Snow cover may affect dust emission in the mid-latitudes during spring, for 169 instance, over northern China (Ginoux and Deroubaix, 2017). The interannual variation 170 of snow cover is also found to affect dust emission in regions, such as Mongolia 171 (Kurosaki and Mikami, 2004). Here monthly snow cover data from MODIS/Terra level 172 3 data (Hall and Riggs, 2015) with a resolution of 0.05° by 0.05° from 2003 to 2015 is 173 used. The high spatial resolution of the product is very suitable for this study. 174 175

Reanalysis 185
Surface wind speed is a critical factor that affects wind erosion. Here 6 hourly 10 186 m wind speed from the NCEP/NCAR reanalysis (Kalnay et al., 1996, hereafter NCEP1) 187 on a T62 Gaussian grid (i.e., 192 longitude grids equally spaced and 94 latitude grids 188 unequally spaced) is used. The NCEP1 is a global reanalysis with relatively long 189 temporal coverage, from 1948 to the present. We chose to use the NCEP1 reanalysis also 190 because surface winds in the GFDL AM4.0 model are nudged toward the NCEP1, and we 191 preferred to use the reanalysis surface wind that is closet to the model climatology. The AErosol RObotic NETwork (AERONET; Holben et al., 1998) provides 202 quality assured cloud-screened (level 2) aerosol measurements from sunphotometer 203 records. In this paper we used the data products of the version 3.0 AERONET processing 204 routine. To examine model simulated DOD, we used coarse mode AOD (COD) at 500 205 nm processed by the Spectral Deconvolution Algorithm (O'Neill et al., 2003;hereafter 206 SDA). SDA COD monthly data is first screened to remove those months with less than 207 Atmos. Chem. Phys. Discuss., https://doi.org/10.5194/acp-2019-223 Manuscript under review for journal Atmos. Chem. Phys. i.e., ! , and Ångström exponent between wavelengths λ A and λ B (α) are used to derive 218 AOD 550 nm using the following equations: 219 (1) 220 (2) 222 223 While this process of extrapolating to 550 nm using a classical Ångström exponent is a 224 bit incoherent with the higher order spectral approach of the SDA, errors due to the 225 choice of spectral order will be negligible in comparison with the types of model versus 226 measurement differences that we will be evaluating in this paper. 227 In a manner similar to the process of screening SDA COD data, monthly AOD 228 550 nm data with less than three days of records in a given month are removed. When 229 calculating the annual means we excluded years having less than five months of records. We also developed a method to derive DOD at 550 nm from AOD 550 nm based 233 on the relationship between Ångström exponent and fine-mode AOD established by 234 Anderson et al. (2005;their Eq. 5). This adds a few more sites over the Sahel than the 235 SDA COD stations. DOD is calculated by subtracting the fine-mode AOD from the total 236 AOD. Due to the large uncertainties of single scattering albedo in AERONET records 237 over regions where AOD is lower than 0.4 (e.g., Dubovik and King, 2000;Holben et al., 238 2006;Andrews et al., 2017), we did not use single scattering albedo to screen AOD to 239 further separate dust from scattering aerosols. Therefore, the derived AERONET DOD 240 over coastal stations may be contaminated by sea salt. 241 242 2) RSMAS surface dust concentration 243 The Rosenstiel School of Marine and Atmospheric Science (hereafter RSMAS 244 dataset) at University of Miami collected mass concentration of dust, sea salt, and sulfate 245 over stations globally, with most of stations on islands (Savoie and Prospero, 1989). The 246 dataset has been widely used for model evaluation (e.g., Ginoux et al., 2001;Huneeus et 247 al., 2011). 248 Only stations with records longer than four years were used and of those stations 249 only those years with at least eight months of data are used for calculating climatological 250 annual means. So, totally 16 stations are used. Station names and locations are listed in 251  (Malm et al., 1994;Hand et al., 2011). IMPROVE stations are located in 260 national parks and wilderness areas, and PM 2.5 sampling is performed twice weekly 261 (Wednesday and Saturday; Malm et al., 1994) prior to 2000 and every third day 262 afterwards. Fine dust (with aerodynamic diameter less than 2.5 µm) concentration is 263 calculated using the concentrations of aluminum (Al), silicon (Si), calcium (Ca), iron 264 (Fe), and titanium (Ti) by assuming oxide norms associated with predominant soil 265 species (Malm et al., 1994;their Eq. 5). This dataset has been widely used to study 266 variations in surface fine dust in the U.S. (e.g., Hand et al., 2016;Hand et al., 2017, Tong 267 et al., 2017Pu and Ginoux, 2018a). Here only monthly data with at least 50% of daily 268 data available in a month (i.e., at least 5 records) are used. Since station coverage over the 269 Step1: Since dust is emitted from the dry and sparsely-vegetated surface, the daily DOD 300 data is first masked out to remove the influences of non-erodible factors and unfavorable 301 environmental conditions that are known to prevent dust emission using criteria as 302 follows: daily VSM less than 0.1 cm 3 cm -3 ; monthly LAI less than 0.3; monthly snow 303 cover less than 0.2%; monthly top-layer soil temperature higher than 273.15 K, i.e., over 304 unfrozen surface; and soil depth thicker than 15 cm. 305 Similar criteria have been used in previous studies to detect or confine dust source 306 regions. For instance, Kim et al. (2013) used NDVI less than 0.15, soil depth greater than 307 10 cm, surface temperature greater than 260 K, and without snow cover to mask 308 topography based dust source function. LAI less than 0.3 has been used as a threshold for 309 dust emission in the Community Land Model (Mahowald et al., 2010;Kok et al., 2014a), 310 while gravimetric soil moisture ranging from 1.01 to 11.2 kg kg -3 depending on soil clay 311 content is recommended to constrain dust emission (Fécan et al., 1999). 312 Step 2: Masked daily DOD from Step 1 is then interpolated to a 0.5° by 0.5° grid using 313 bilinear interpolation. This is close to the horizontal resolution of the GFDL 314 AM4.0/LM4.0 model used in this study. Then the cumulative frequency distribution of 315 daily DOD from 2003 to 2015 is derived at each grid point for each month. 316 Step 3: Daily maximum surface wind speed is first derived from 6-hourly NCEP1 surface 317 winds and then interpolated to a 0.5° by 0.5° grid. The cumulative frequency distribution  Step 4: A minimum value of DOD (DOD thresh ) is used to separate dust events from 321 background dust. The cumulative frequency (in %) of dust events passing this threshold 322 is compared to the cumulative frequency of surface winds. The minimum surface winds 323 with the same frequency correspond to the threshold of wind erosion, V threshold (see a 324 schematic diagram in Figure S1 in the Supplement). This operation is performed for all 325 grid points for each month. Ginoux et al. (2012) used DOD thresh = 0.2 to quantify the FoO 326 of local dust events. Similarly, DOD thresh = 0.2 is used here in major dusty regions (North 327 Africa, Middle East, India, northern China), while for less dusty regions, such as the U.S., 328 South America, South Africa, and Australia, DOD thresh = 0.02 is used. The reason to use a 329 lower DOD thresh for less dusty regions is because: i) the overall dust emission in these 330 regions are at least ten times smaller than major dusty regions, such as North Africa (e.g., 331 Huneeus et al., 2011); ii) the frequency distribution of DOD in these regions also peaks at 332 a much lower DOD band (see discussion in section 3.3). 333 Figures 1a-e show the seasonal and annual mean FoO (days when DOD is greater 334 than DOD thresh ) using the DOD thresh defined here. The shaded area covers major dust 335 sources, and the pattern is very similar to that obtained by Ginoux et al. (2012;their Fig. 336 5), although there are some differences, largely due to the masked DOD (i.e., from Step 337 1) used in this study and a lower threshold in less dusty regions. 338 Note that the selections of masking criteria in Step 1 and DOD thresh in Step 4 are 339 empirical and can add uncertainties to this method. Also, we approximate dust emission 340 using cumulative frequency of DOD, which may overestimate dust emission in regions 341 where the contribution of transported dust is significant and thus underestimate the 342 V threshold in those regions. 343

Simulation design 344
The AM4.0/LM4.0 is a coupled land-atmosphere model newly developed at 345 GFDL (Zhao et al., 2018a,b). It uses the recent version of the GFDL Finite-Volume 346 Cubed-Sphere dynamical core (FV 3 ; Putman and Lin, 2007), which is developed for 347 weather and climate applications with both hydrostatic and non-hydrostatic options. The aerosol physics is based in large part on that of the GFDL AM3.0 (Donner et 356 al., 2011), but with a simplified chemistry where ozone climatology from AM3.0 357 simulation (Naik et al., 2013) is prescribed. AM4.0 simulates the mass distribution of five 358 aerosols: sulfate, black carbon, organic carbon, dust, and sea salt. Dust is partitioned into 359 five size bins based on radius: 0.1~1 µm (bin 1), 1~2 µm (bin 2), 2~3 µm (bin 3), 3~6 µm 360 (bin 4), and 6~10 µm (bin 5). The dust emission scheme follows the parameterization of 361 Ginoux et al. (2001), as shown in the following equation: 362 where F p is flux of dust of particle size class p, C is a scaling factor with a unit of µg s 2 365 m -5 , here C is set to 0.75×10 -9 . S is the source function based on topographic depressions 366  (Ginoux et al., 2001), s p is fraction of each size class, and V 10m is surface 10 m wind 367 speed, and V t = 6 m s -1 is the threshold of wind erosion. 368 Three simulations with prescribed sea surface temperature (SST) and sea ice 369 (Table 2) were conducted from 1999 to 2015, with the first year discarded for spin up. 370 The Atmospheric Model Intercomparison Project (AMIP)-style SST and sea ice data 371 (Taylor et al., 2000) are from the Program for Climate Model Diagnosis and 372 Intercomparison (PCMDI), which combined HadISST (Rayner et al., 2003)

from UK Met 373
Office before 1981 and NCEP Optimum Interpolation (OI) v2 SST (Reynolds et al., 374 2002) afterwards. The surface winds in the simulations are nudged toward the NCEP1 375 reanalysis with a relaxation timescale of 6 hours (Moorthi and Suarez, 1992). Note that 376 the nudged surface winds are actually weaker than the surface wind speed simulated by 377 the standard version of AM4.0/LM4.0 without nudging, so the overall magnitude of dust 378 emission is lower than the standard version. Here we choose not to retune the dust 379 emission scheme but instead test the usage of V threshold , which theoretically provides a 380 more physics-based way to improve dust simulation. 381 In the Control run, the default model setting is used for dust emission, with a 382 prescribed 6 m s -1 threshold of wind erosion (cf. Ginoux et al., 2019). In the V thresh 12mn 383 simulation, the observation based climatological monthly V threshold is used to replace the 384 constant wind erosion threshold. The default source function S in Eq. 3 only allows dust 385 emission over bare ground by masking out regions with vegetation cover. Since LAI 386 masking is already applied in the retrieval of V threshold (i.e., LAI<0.3), we choose to use a 387 source function that is the same as the default source function S but without vegetation  The distributions of V threshold for annual mean (black bars) and dusty seasons 419 (color lines; MAM and JJA for the Northern Hemisphere and SON and DJF for the 420 Southern Hemisphere) for each dusty region (see Fig. 1f and Table 1 for locations) are 421 shown in Figs. 2b-j. In the Sahel, the annual mean V threshold peaks around 4 m s -1 (Fig. 2b). 422 This magnitude is lower than indicated from previous studies based on station 423 observations in the region, e.g., Helgren and Prospero (1987) found the threshold velocity 424 over eight stations in Northwest Africa ranged from 6.5 to 13 m s -1 during summer in 425 1974. Chomette et al. (1999) and Marsham et al. (2013) also reported higher wind 426 erosion thresholds around 6-9 m s -1 at individual stations. On the other hand, Cowie et al. 427 (2014) found that the annual threshold of wind erosion at the 25% level, i.e., when 428 surface condition is favorable for dust emission, can be lower than 6 m s -1 at some sites in 429 the Sahel (their Fig. 5). Several factors may contribute to the discrepancies. First, studies 430 suggest that reanalysis datasets may underestimate surface wind speed in spring and for 431 monsoon days in Africa (e.g., Largeron et al., 2015), and therefore could lead to a lower 432 value of V threshold than that derived from station observations. In fact, Bergametti et al. The annual mean V threshold in the Sahara and Arabian Peninsula is a bit higher, 440 with mean values at 4.5 and 5.2 m s -1 , respectively (Figs. 2c-d). The V threshold over 441 northern China is even higher, with an annual mean of 7.9 m s -1 . This is consistent with 442 the results of Kurosaki and Mikami (2007), who found that under favorable land surface 443 conditions the threshold wind speed ranges from 4.4± 0.6 m s -1 in Taklimakan Desert to 444 6.9± 1.2 m s -1 over the Loess Plateau and around 9.8± 1.6 m s -1 in the Gobi Desert. These 445 values are also consistent with Ginoux and Deroubaix (2017) who found that regional 446 mean wind erosion threshold over northern China ranges from 6.5 to 9.1 m s -1 . In India, 447 the V threshold peaks at about 4.5 m s -1 and 6.5 m s -1 , respectively (Fig. 2f). The second peak 448 is probably related to anthropogenic dust sources over the central Indian subcontinent 449 (Ginoux et al., 2012). We also note that in the Northern Hemisphere, the V threshold in dusty 450 seasons is shifted towards lower values than the annual mean (blue and green lines in 451

Climatology of AOD and DOD 461
In order to compare the model results with observations, we first show the 462 climatology of AERONET AOD and COD from 2003 to 2015. As shown in Figure 3, 463 annual mean global AOD is highest over Africa, the Arabian Peninsula, Indian 464 subcontinent, and Southeast Asia. In the latter two regions, high sulfate concentrations 465 (e.g., Ginoux et al., 2006) and organic carbon from biomass burning in Southeast Asia 466 (e.g., Lin et al., 2014) contribute substantially to the total AOD. The SDA COD shows 467 the optical depth due to coarse aerosols, which includes both dust and sea salt, and sea 468 salt over coastal regions or islands can be a major contributor. Here, high values (>0.2) 469 are largely located over dusty regions such as North Africa, the Arabian Peninsula, and 470 northern India (Fig. 3b). 471 percentage of DOD to total COD in the model is displayed at the bottom (Fig. 4e). The 476 simulated AOD is lower than that from the AERONET over North Africa, the Middle 477 East, and western India, largely due to low values of COD simulated in these regions 478 (Fig. 4d). Besides these regions, the COD over North America, South America, South 479 Africa, and northern Eurasia is also, for the most part, underestimated by the model. Dust The underestimation of COD (and effectively DOD given its dominance in most 483 regions) was improved in the subsequent model run using a prescribed 12-month V threshold . 484 AERONET COD (Fig. 5d). Over the Indian subcontinent, COD is overestimated, while 488 over North America excluding the east coast, northern Eurasia, and part of South 489 America, COD is also better captured than in the Control run. 490 These improvements are largely associated with a better simulation of DOD in the 491 "dust belt" (i.e., North Africa and the Middle East). Figure 6 shows the DOD at 550 nm 492 derived from AERONET AOD (see methodology for details) versus that from the 493 V thresh 12mn simulation. Over most stations in the Sahel, Mediterranean coasts, and 494 central Middle East, the relative differences between modeled and observed DOD is 495 within ± 25%. 496 Figure 7 shows the regional averaged annual mean DOD over nine dusty regions 497 from MODIS and three simulations. The Control run largely underestimates DOD in all 498 regions, while the magnitude of DOD is better captured in the V thresh 12mn and V thresh Ann 499 simulations, although slightly overestimated in the Sahel and greatly overestimated over 500 Australia. In general, DOD simulated by the V thresh Ann run using a constant annual mean 501 V threshold is higher than that simulated by the V thresh 12mn run, consistent with the higher 502 dust emission in the V thresh Ann run (Table S2 in  sites 16 and 13), dust concentrations are overestimated, consistent with the 517 overestimation of DOD over the U.S. and the Sahel (Fig. 7). Dust concentrations in 518 Australia and the east coast of China are also overestimated by more than five-folds. 519 Surface dust concentration is further improved in the V thresh 12mn simulation (Fig. 8,  520 bottom), with eight stations showing a model/observation ratio between 0.5 and 2 and 521 only four stations overestimating or underestimating dust concentrations by more than 522 five times. 523 Simulated surface fine dust concentration (calculated as dust bin 1+0.25×dust bin 524 2) in the U.S. is compared with gridded IMPROVE data (Figure 9). While the Control 525 run largely underestimates surface fine dust concentration, the simulated concentration is 526 Atmos. Chem. Phys. Discuss., https://doi.org/10.5194/acp-2019-223 Manuscript under review for journal Atmos. Chem. Phys. Sahel, the Sahara, and the Arabian Peninsula (Figs. 10a-c), although the spring and 545 summer peak in the Sahel is overestimated and winter minimum in the Sahara is 546 underestimated. The MAM peak of MODIS DOD in northern China is missed by both 547 V thresh 12mn and V thresh Ann simulations (Fig. 10d), while the JJA peak over India is 548 largely overestimated (Fig. 10e). Over the U.S. dusty region, the seasonal cycle in the 549 Atmos. Chem. Phys. Discuss., https://doi.org/10.5194/acp-2019- V thresh 12mn simulation is slightly underestimated compared to MODIS DOD but 550 overestimated from May to August in the V thresh Ann simulation (Fig. 10f). DOD is 551 underestimated in South Africa in all three simulations (Fig. 10g). Over South America, 552 the peak from October to February is roughly captured by the V thresh 12mn run but is 553 overestimated by the V thresh Ann run (Fig. 10h). The seasonal cycles of DOD in Australia 554 are very similar in all three simulations and largely resemble that in the MODIS, although 555 both the V thresh 12mn and V thresh Ann simulations overestimate the DOD by about an order 556 of magnitude. 557 Figure 11 shows the seasonal cycle of COD from 12 AERONET SDA sites over 558 North Africa and nearby islands (see Figure S5 in the Supplement for site locations) 559 along with MODIS DOD and DOD simulated in three runs. The magnitude of 560 AERONET COD and MODIS DOD in these sites are very similar, despite missing values 561 at sites 1, 4, 5, 8, 11, and a smaller value at site 2 in MODIS. Over most of the sites, the 562 seasonal cycle is better captured in the V thresh 12mn and V thresh Ann simulations than the 563 Control run, although the peak over Cairo_EMA_2 (site 12) is slightly underestimated, 564 which is consistent with the underestimation of annual mean DOD in the area (Fig. 6). 565 We also examined the seasonal cycle of PM10 surface concentration at three 566 Sahelian INDAAF stations (see Figure S5  winter that transports Saharan dust southward to the Guinean coast and the scavenging 573 effect of monsoonal rainfall in boreal summer that removes surface dust (Marticorena et 574 al., 2010). While the Control run does not capture the seasonal cycles in these sites, the 575 V thresh 12mn run largely captures the spring peak and summer minimum, although the 576 magnitude is overestimated. In all three sites, the simulated concentration in the 577 V thresh Ann run is larger than that in the V thresh 12mn run, especially in boreal fall to early 578 spring. Such an overestimation is probably due to the prescribed constant annual mean 579 V threshold , which is lower than it would be during the less dusty season (i.e., boreal fall to 580 winter) and thus increases dust emission and surface concentration. surface winds. However, this storm was not predicted by the forecast models, such as the 597 Goddard Earth Observing System version 5 (GEOS-5; Rienecker et al., 2008) and Navy 598 Aerosol Analysis Prediction System (NAAPS; Witek et al., 2007;Reid et al., 2009;599 Westphal et al., 2009). 600 As shown in Figure 13, MODIS DOD also captures this event, with a peak value 601 above 0.5 over southwest Nebraska and northern Kansas on Oct. 18 th , 2012. The 602 V thresh 12mn run also largely captures this event (Fig. 13 bottom panel), although the 603

Control run totally misses it (not shown). In the model, the dust storm appears in South 604
Dakota and Nebraska on Oct. 17 th , 2012, along with the anomalous southwesterly winds. 605 It reaches a maximum on Oct. 18 th , in association with intensified anomalous 606 southwesterly winds at the surface and an anomalous low-pressure system at 850 hPa 607 ( Figure S6 in the Supplement). Note that the modeled dust storm centers a bit 608 northeastward compared to the MODIS DOD pattern and it also has greater magnitude 609 and covers a larger area. On Oct. 19 th , both the anomalous low-pressure system and 610 surface wind speeds weaken and the dust storm dissipates, with slightly elevated DOD 611 levels over a region extending over the lower Mississippi River basin and the Midwest. 612 This is somewhat consistent with MODIS records, which also shows slightly higher DOD 613 levels over Tennessee and northern Alabama on Oct. 19 th , regardless of large area of 614 missing values.  Figure 14 shows the frequency of regional mean DOD during one dusty season 618 (MAM in the Northern Hemisphere and SON in the Southern Hemisphere) for nine 619 regions. Results from MODIS, the Control, and V thresh 12mn runs are shown in black, 620 blue, and orange lines, respectively. In most dusty regions, such as the Sahara, Sahel, 621 Arabian Peninsula, India, and northern China, MODIS DOD frequency largely peaks 622 between 0.2 to 0.4, while DOD frequency peaks at a much lower level between 0.02 to 623 0.08 in less dusty regions, such as the U.S., South America, South Africa and Australia. 624 The DOD distribution in the Control run is biased low and peaks around 0.05 in those 625 dusty regions and between 0 and 0.01 in less dusty regions. The frequency is much better 626 captured in the V thresh 12mn run over the Arabian Peninsula and the Sahel, slightly 627 improved but still biased low over the Sahara, northern China, India, and the U.S. The 628 modeled frequency in the V thresh 12mn run is biased high in Australia (peaks outside the 629 maximum of x-axis, not shown) and shows little improvement over South Africa and 630 South America. The overall improvement of DOD frequency using the time-varying 2D 631 V threshold occurs mostly over major dusty regions, which is consistent with the 632 improvements in DOD climatology and seasonal cycle in the model simulations. MODIS and other satellite products, thus the uncertainties in the satellite products are 641 inherited in the derived DOD frequency distribution. Due to the cloud screening 642 processes of MODIS products, dust activities over cloud-covered regions may be 643 underestimated. Also, DOD frequency is derived based on daily observations over a 13-644 year record, so that some variability of dust emission associated with alluvial sediments 645 deposited by seasonal flooding may be not captured. Diurnal variability of dust emission 646 and short-duration events such as haboobs are also not included. Since DOD is a column 647 integrated variable, it includes both local emitted and remotely transported dust. When 648 using DOD frequency distribution to approximate dust emission, it may overestimate dust 649 emission in regions where transported dust is dominated, e.g., over the southern Sahel, 650 and lead to an underestimation of V threshold . 651 Previous study found that over regions such as North Africa, reanalysis products 652 may underestimate surface wind speed in spring and monsoon seasons but overestimate it 653 during dry nights (e.g., Largeron et al., 2015). This is largely because mechanisms such 654 as density current that can enhance surface wind speed are not parameterized in the 655 atmospheric models to produce the reanalysis products, while coarse spatial and temporal 656 sampling may also contribute to the underestimation of reanalysis wind speeds. These 657 limitations add uncertainties to the V threshold estimates derived here. 658 In addition, V threshold is derived by matching the frequency distribution of DOD at 659 certain levels (0.2 or 0.02) with the frequency distribution of daily maximum wind, and 660 these two values are derived empirically. The influences of soil properties such as soil 661 cohesion, particle size, and particle compositions on the threshold of wind erosion (e.g., 662 Atmos. Chem. Phys. Discuss., https://doi.org/10.5194/acp-2019-223 Manuscript under review for journal Atmos. Chem. Phys. Discussion started: 19 March 2019 c Author(s) 2019. CC BY 4.0 License. Fécan et al., 1999;Alfaro and Gomes, 2001;Shao, 2001;Kok et al., 2014b) are not 663 explicitly examined here and will need further investigation.

Conclusion 681
While dust aerosols play important roles in the Earth's climate system, large 682 uncertainties exist in modeling its lifecycle (e.g., Huneeus et al., 2011;Pu and Ginoux, 683 2018b). Constant thresholds of wind erosion are widely used in climate models for 684 simplicity. Here, high-resolution MODIS Deep Blue dust optical depth (DOD) and 685 Atmos. Chem. Phys. Discuss., https://doi.org/10.5194/acp-2019-223 Manuscript under review for journal Atmos. Chem. Phys.  Table 1 Major dusty regions shown in Figure 1. Note that region names such as India and 1189 northern China are not exactly the same as their geographical definitions but also cover 1190 some areas from nearby countries. 1191 1192 Table 2 Simulation design 1193   Table 1 Major dusty regions shown in Figure 1. Note that region names such as India and 1303 northern China are not exactly the same as their geographical definitions but also cover 1304 some areas from nearby countries.