Retrieving the global distribution of the threshold of wind erosion from satellite data and implementing it into the Geophysical Fluid Dynamics Laboratory land–atmosphere model (GFDL AM4.0/LM4.0)

Dust emission is initiated when surface wind velocities exceed the threshold of wind erosion. Many dust models used constant threshold values globally. Here we use satellite products to characterize the frequency of dust events and land 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 the wind erosion threshold (V_threshold) over dry and sparsely vegetated surfaces. 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 PM₁₀ dust concentrations, and the seasonal cycle of DOD are better captured over the “dust belt” (i.e., northern 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 V_threshold, 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 cycles as well as dust forecasting.

1 Introduction

Mineral dust is one of the most abundant aerosols by mass and plays an important role in the climate system. Dust particles absorb and scatter solar and terrestrial radiation, thus modifying the local energy budget and consequently atmospheric circulation patterns. Studies have shown that the radiative effect of dust can affect a wide range of environmental processes. Dust is shown to modulate western African (e.g., Miller and Tegen, 1998; Miller et al., 2004; Mahowald et al., 2010; Strong et al., 2015) and Indian (e.g., Jin et al., 2014, 2015, 2016; Vinoj et al., 2014; Solmon et al., 2015; Kim et al., 2016; Sharma and Miller, 2017) monsoonal precipitation. During severe droughts in North America, there is a positive feedback between dust and the hydrological cycle (Cook et al., 2008, 2009, 2013). African dust is also found to affect Atlantic tropical cyclone activity (e.g., Dunion and Velden, 2004; Wong and Dessler, 2005; Evan et al., 2006; Strong et al., 2018). When deposited on snow and ice, dust reduces the surface reflectivity, enhancing net radiation and accelerating snow and ice melting and consequently affecting runoff (e.g., Painter et al., 2010, 2018; Dumont et al., 2014). Dust can serve as ice nuclei and affect the formation, lifetime, and characteristic of clouds (e.g., Levin et al., 1996; Rosenfield et al., 1997; Wurzler et al., 2000; Nakajima et al., 2001; Bangert et al., 2012), perturbing the hydrological cycle. Iron- and phosphorus-enriched dust is also an important nutrient for the marine and terrestrial ecosystems and thus interacts with the ocean and land biogeochemical cycles (e.g., Fung et al., 2000; Jickells et al., 2005; Shao et al., 2011; Bristow et al., 2010; Yu et al., 2015).

Given the importance of mineral dust, many climate models incorporate dust emission schemes to simulate the life cycle of dust aerosols (e.g., Donner et al., 2011; Collins et al., 2011; Watanabe et al., 2011; Bentsen et al., 2013). Mineral dust particles are lifted from dry and bare soils into the atmosphere by saltation and sandblasting. This process is initiated when surface winds reach a threshold velocity of wind erosion. The value of this wind erosion threshold depends on soil and surface characteristics, including soil moisture, soil texture, and particle size, and the presence of pebbles, rocks, and vegetation residue (e.g., Gillette et al., 1980; Gillette and Passi, 1988; Raupach et al., 1993; Fécan et al., 1999; Zender et al., 2003; Mahowald et al., 2005), and this thus varies spatially and temporally (Helgren and Prospero, 1987). Due to a lack of in situ data at a global scale and uncertainties on these dependencies, most dust and climate models prescribe a spatially and temporally constant threshold of wind erosion for surface 10 m wind (e.g., around 6 to 6.5 m s⁻¹) over dry surfaces for simplicity (e.g., Tegen and Fung, 1994; Takemura et al., 2000; Uno et al., 2001; Donner et al., 2011). For instance, in the Geophysical Fluid Dynamics Laboratory coupled land–atmosphere model AM4.0/LM4.0 (Zhao et al., 2018a, b), a constant threshold of 6 m s⁻¹ is used. On the other hand, some models, such as the Hamburg version of the European Centre for Medium-Range Weather Forecasts (ECMWF) model Hamburg Aerosol Module (ECHAM-HAM), Hadley Centre Global Environmental Model, version 2, Earth System (HadGEM2-ES), and ICOsahedral Nonhydrostatic – Aerosol and Reactive Trace gases (ICON-ART), parameterize the constant dry threshold friction velocity (usually a function of soil particle size, soil, and air density) or threshold wind velocity with dependencies on soil moisture, surface roughness length, and vegetation coverage (e.g., Takemura et al., 2000; Ginoux et al., 2001; Zender et al., 2003; Cheng et al., 2008; Jones et al., 2011; Rieger et al., 2017).

The threshold of wind erosion may be approximately inferred using observations. For instance, Chomette et al. (1999) used the Infrared Difference Dust Index (IDDI) and 10 m winds reanalysis from the ECMWF between 1990 and 1992 to calculate the threshold of wind erosion over seven sites over the Sahel and Sahara. The IDDI was used to determine whether there was a dust event for subsequently calculating an emission index defined as the number of dust events to the total number of potential events. The distribution of surface wind speed was matched with the emission index, and the threshold of wind erosion was determined when the emission index was around 0.9. The resulting average threshold of wind erosion ranged from 6.63 m s⁻¹ at a Sahelian site to about 9.08 m s⁻¹ at a Niger site, consistent with the model results by Marticorena et al. (1997).

Later, Kurosaki and Mikami (2007) used World Meteorological Organization (WMO) station data from March 1998 to June 2005 to examine the threshold wind speed in eastern Asia. Using the distribution of surface wind speed and associated weather conditions (i.e., with or without dust emission events), they approximated a dust emission frequency by dividing the number of dust events by the total number of observations for each wind bin, and then they determined threshold wind speeds at the 5 % and 50 % levels, corresponding to the most favorable and normal land surface conditions for dust emission, respectively. They found that the derived threshold wind speed varied in space and time, with a larger seasonal cycle in grassland regions, such as northern Mongolia, and smaller seasonal variations in desert regions, such as the Taklimakan and Gobi deserts and the Loess Plateau. Cowie et al. (2014) applied a similar method over northern Africa, using wind data observed between 1984 and 2012, and they focused on threshold winds at the 25 %, 50 %, and 75 % levels.

Draxler et al. (2010) derived the distribution of the threshold of friction velocity over the US by matching the frequency of occurrence (FoO) of Moderate Resolution Imaging Spectroradiometer (MODIS) Deep Blue (Hsu et al., 2004) aerosol optical depth (AOD) above 0.75 with the FoO of friction velocities extracted from the North American Mesoscale (NAM) forecast model at each grid point. This new threshold and a soil characteristics factor was then incorporated into the Hybrid Single-Particle Lagrangian Integrated Trajectory (HYSPLIT) model (Draxier and Hess, 1998) to forecast dust surface concentrations. It was found that major observed dust plume events in June and July 2007 were successfully captured by the model. Later, Ginoux and Deroubaix (2017) used FoO derived from the MODIS Deep Blue dust optical depth (DOD) record to retrieve the wind erosion threshold of surface 10 m winds over eastern Asia.

For individual dust events, the threshold of friction velocity can also be determined by fitting a 2nd-order Taylor series to dust saltation flux measurements (Barchyn and Hugenholtz, 2011; Kok et al., 2014b).

Nonetheless, a global distribution of the threshold of wind erosion with observational constraints that may be implemented in climate models is still lacking. In this study, we propose a method to retrieve the monthly global threshold of wind erosion (hereafter, V_threshold) for dry and sparsely vegetated surface (i.e., under favorable conditions for dust emission) using high-resolution satellite products and reanalysis datasets. This two-dimensional threshold of surface 10 m winds is then implemented into the Geophysical Fluid Dynamics Laboratory (GFDL) coupled land–atmosphere model, AM4.0/LM4.0 (Zhao et al., 2018a, b). The benefits of using this spatially and temporally varying threshold in simulating present-day climatology and seasonal cycles of dust are analyzed by comparing the model results with observations.

The data and method used to retrieve the threshold of wind erosion are detailed in Sect. 2. The distribution of the derived V_threshold and its implication in the climate model is presented in Sect. 3. Section 4 discusses the uncertainties associated with this method, and major conclusions are summarized in Sect. 5.

2 Data and methodology

In this section we first introduce the satellite products, observational data, and reanalyses used to retrieve the threshold of wind erosion and validate model output (Sect. 2.1). The processes to retrieve the threshold of wind erosion are detailed in Sect. 2.2. The uncertainties of V_threshold associated with the retrieval criteria and selection of surface wind datasets are discussed in Sect. 2.3. Section 2.4 introduces the GFDL AM4.0/LM4.0 model, its dust emission scheme, and simulation designs.

2.1 Data

2.1.1 Satellite products

1. MODIS Aqua and Terra dust optical depth. DOD is column-integrated extinction by mineral particles. Here daily DOD is retrieved from MODIS Deep Blue aerosol products (collection 6, level 2; Hsu et al., 2013; Sayer et al., 2013): aerosol optical depth (AOD), single-scattering albedo (ω), and the Ångström exponent (α). All the daily variables are first interpolated to a 0.1 × 0.1^∘ grid using the algorithm described by Ginoux et al. (2010). We require the single-scattering albedo at 470 nm to be less than 0.99 for dust due to its absorption of solar radiation. This separates dust from scattering aerosols, such as sea salt. Then a continuous function relating the Ångström exponent, which is highly sensitive to particle size (Eck et al., 1999), to fine-mode AOD established by Anderson et al. (2005; their Eq. 5) is used to separate dust from fine particles. In short, DOD is retrieved using the following equation:

   $$\begin{array}{}\text{(1)}& \text{DOD}=\text{AOD}\times (\mathrm{0.98}-\mathrm{0.5089}\mathit{\alpha }+\mathrm{0.0512}{\mathit{\alpha }}^{\mathrm{2}}).\end{array}$$

   Details about the retrieval process and estimated errors are summarized by Pu and Ginoux (2018b). High-resolution MODIS DOD products (0.1 × 0.1^∘) have been used to identify and characterize dust sources (Ginoux et al., 2012; Baddock et al., 2016) and examine the variations in dustiness in different regions (e.g., Pu and Ginoux, 2016, 2017, 2018b).

   Following the recommendation from Baddock et al. (2016), who found the dust sources are better detected using DOD with a low-quality flag (i.e., quality assurance flag, QA, equals 1, following the category of retrieval quality flags in MODIS Deep Blue products; Hsu et al., 2013) than with a high-quality flag (i.e., QA = 3); as retrieved aerosol products were poorly flagged over dust source regions, we also use DOD with the flag of QA = 1. Both daily DOD retrieved from Aqua and Terra platforms are used by averaging the two when both products are available or using either one when only one product is available. Since Terra passes the Equator from north to south around 10:30 local time (LT) and Aqua passes the Equator from south to north around 13:30 LT, an average of the two combines the information from both morning and afternoon hours. This process also largely reduces missing data (Pu and Ginoux, 2018b). This combined daily DOD, hereafter MODIS DOD, is available from January 2003 to December 2015 at a resolution of 0.1 × 0.1^∘ grid. Note that due to the temporal coverage of MODIS products, the diurnal variations in dust (e.g., O'rgill and Sehmel, 1976; Mbourou et al., 1997; Knippertz, 2008; Schepanski et al., 2009) are not included in current study.

2. Soil moisture. Soil moisture is an important factor that affects dust emission (Fécan et al., 1999). Daily surface volumetric soil moisture (VSM) retrievals derived from similar calibrated microwave (10.7 GHz) brightness temperature observations from the Advanced Microwave Scanning Radiometer for the Earth Observing System (AMSR-E) onboard the NASA Aqua satellite (from June 2002 to October 2011) and the Advanced Microwave Scanning Radiometer 2 (AMSR2) sensor onboard the JAXA GCOM-W1 satellite (from July 2012 to June 2017) from the University of Montana (Du et al., 2017a, b) was used to retrieve the wind erosion threshold. Both AMSR-E and AMSR2 sensors provide global measurements of polarized microwave emissions at six channels, with ascending and descending orbits crossing the Equator at around 13:30 and 01:30 LT respectively. The VSM retrievals are derived from an iterative retrieval algorithm that exploits the variable sensitivity of different microwave frequencies and polarizations and minimizes the potential influence of atmosphere, vegetation, and surface water cover on the soil signal. The VSM record represents surface (top ∼ 2 cm) soil conditions and shows favorable global accuracy and consistent performance (Du et al., 2017b), particularly over areas with low-to-moderate vegetation cover that are also more susceptible to wind erosion, although cautions are needed when examining long-term trends due to the small biases between AMSR-E and AMSR2. The horizontal resolution of the product is about 25 km×25 km, and the daily product from January 2003 to December 2015 is used. The ascending and descending obit VSM retrievals are averaged to get the mean VSM for each day.

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

4. Leaf area index (LAI). Vegetation protects soil from the effects of wind and thus modulates dust emission (e.g., Marticorena and Bergametti, 1995; Zender et al., 2003). While dense vegetation coverage can increase surface roughness and reduce near-surface wind speed, the roots of vegetation can increase soil cohesion and further reduce wind erosion. LAI describes the coverage of vegetation with a unit of m² m⁻², i.e., leaf area per ground area. Here monthly LAI retrieved by Boston University from MODIS onboard Aqua (Yan et al., 2016a, b; Ranga Myneni and Taejin Park, Boston University, personal communication, 2016) with a resolution of 0.1 × 0.1^∘ from 2003 to 2015 is used. The root mean square error of the product is 0.66, with some overestimation of LAI in sparsely vegetated regions (Yan et al. 2016b; Garrigues et al., 2008).

2.1.2 Reanalysis

Surface wind speed is a critical factor that affects wind erosion. Here 6-hourly 10 m wind speed from the NCEP/NCAR (National Centers for Environmental Prediction/National Center for Atmospheric Research) reanalysis (Kalnay et al., 1996, hereafter NCEP1) on a T62 Gaussian grid (i.e., 192 longitude grids equally spaced and 94 latitude grids unequally spaced) is used. The NCEP1 is a global reanalysis with relatively long temporal coverage, from 1948 to the present. We chose to use the NCEP1 reanalysis mainly because surface winds in the GFDL AM4.0 model are nudged toward the NCEP1, and we preferred to use the reanalysis surface wind that is closest to the model climatology.

ERA-Interim (ECMWF Reanalysis Interim; Dee et al., 2011) is a global reanalysis produced from ECMWF. It provides high spatial resolution (about 0.75^∘ or 80 km) 6-hourly, daily, and monthly reanalysis from 1979 to present day. Soil temperature from the ERA-Interim is used to determine the regions where wind erosion may be prohibited by the frozen surface. Monthly temperature of the first soil layer (0 to 0.07 m) from 2003 to 2015 is used.

In order to quantify the uncertainties of the retrieved threshold wind erosion in association with the selection of reanalysis products, surface 10 m winds from 6-hourly ERA-Interim and hourly ERA5 (Hersbach and Dee, 2016) are both examined. The ERA5 is the latest reanalysis product from the ECMWF, with a horizontal resolution of about 31 km and hourly temporal resolution.

2.1.3 Station data

Multiple ground-based datasets are used to validate AM4.0/LM4.0 simulated aerosol and dust optical depth and surface dust concentrations.

1. AERONET. The AErosol RObotic NETwork (AERONET; Holben et al., 1998) provides quality-assured cloud-screened (level 2) aerosol measurements from sun photometer records. In this paper we used the data products of the version 3.0 AERONET processing routine. To examine model-simulated DOD, we used coarse-mode AOD (COD; i.e., radius > 0.6 µm) at 500 nm processed by the Spectral Deconvolution Algorithm (O'Neill et al., 2003; hereafter SDA). SDA COD monthly data are first screened to remove those months with less than 5 d of records. To get the annual means, years with less than 5 months of records were removed. Only stations with records of at least 3 years during the period were used to calculate the 2003–2015 climatology (the same time period when MODIS DOD is available). Overall, records from 313 stations were obtained.

   AERONET monthly aerosol optical thickness (AOT) data around 550 nm (e.g., 500, 551, 531, 440, 675, 490, 870 nm, etc.) and the Ångström exponents across the dual wavelength of 440–675, 440–870, and 500–870 nm are used to calculate AOD at 550 nm (τ₅₅₀). If AOT for 551, 555, 531 or 532 nm exist, then these values are directly used as AOD 550 nm. Otherwise, the AOT at wavelength λ_A (less than 550 nm; i.e., τ_A), AOT at wavelength λ_B (larger than 550 nm; i.e., τ_B), and Ångström exponent between wavelengths λ_A and λ_B (α) are used to derive AOD 550 nm using the following equations:

   $$\begin{array}{}\text{(2)}& {\displaystyle }& {\displaystyle }{\mathit{\tau }}_{\mathrm{550}}={\mathit{\tau }}_A{\left({\displaystyle \frac{\mathrm{550}}{{\mathit{\lambda }}_A}}\right)}^{-\mathit{\alpha }}\text{if}{\mathit{\tau }}_A\text{i s available},\text{(3)}& {\displaystyle }& {\displaystyle }{\mathit{\tau }}_{\mathrm{550}}={\mathit{\tau }}_B{\left({\displaystyle \frac{\mathrm{550}}{{\mathit{\lambda }}_B}}\right)}^{-\mathit{\alpha }}\text{if}{\mathit{\tau }}_B\text{is available}.\end{array}$$

   In a manner similar to the process of screening SDA COD data, monthly AOD 550 nm data with less than 3 d of records in a given month are removed. When calculating the annual means we excluded years having less than 5 months of records. Finally, to calculate the climatology of 2003–2015, only stations with at least 3 years of records during this period are used, totaling 351.

   We also developed a method to derive DOD at 550 nm from AOD at 550 nm based on the relationship between the Ångström exponent and fine-mode AOD established by Anderson et al. (2005; their Eq. 5). This adds a few more sites over the Sahel than the SDA COD stations. DOD is calculated by subtracting fine-mode AOD from total AOD. Due to the large uncertainties of single-scattering albedo in AERONET records over regions where AOD is lower than 0.4 (e.g., Dubovik and King, 2000; Holben et al., 2006; Andrews et al., 2017), we did not use single-scattering albedo to screen AOD to further separate dust from scattering aerosols. Therefore, the derived AERONET DOD over coastal stations may be contaminated by sea salt.

2. RSMAS surface dust concentration. The Rosenstiel School of Marine and Atmospheric Science (hereafter the RSMAS dataset) at the University of Miami collected mass concentrations of dust, sea salt, and sulfate over stations globally, with most of the stations on islands (Savoie and Prospero, 1989). The dataset has been widely used for model evaluation (e.g., Ginoux et al., 2001; Huneeus et al., 2011).

   Only stations with records longer than 4 years were used, and of those stations only those years with at least 8 months of data are used for calculating climatological annual means. So, in total 16 stations are used. Station names, locations, and record lengths are listed in Table S1 of the Supplement. We compare the climatology of annual mean surface dust concentration with model output during 2000–2015. Note that since most station records end earlier than 1998, the dataset largely represents the climatology during the 1980s and 1990s. Thus the discrepancies between model output and the RSMAS data include both model biases and the difference in surface dust concentration from the 1980s to the 2000s.

3. IMPROVE surface fine dust concentration. The Interagency Monitoring of Protected Visual Environments (IMPROVE) network has collected near-surface particulate matter 2.5 (PM_2.5) samples in the US since 1988 (Malm et al., 1994; Hand et al., 2011). IMPROVE stations are located in national parks and wilderness areas, and PM_2.5 sampling is performed twice weekly (Wednesday and Saturday; Malm et al., 1994) prior to 2000 and every third day afterwards. Fine dust (with an aerodynamic diameter less than 2.5 µm) concentration is calculated using the concentrations of aluminum (Al), silicon (Si), calcium (Ca), iron (Fe), and titanium (Ti) by assuming oxide norms associated with predominant soil species (Malm et al., 1994; their Eq. 5). This dataset has been widely used to study variations in surface fine dust in the US (e.g., Hand et al., 2016, 2017, Tong et al., 2017; Pu and Ginoux, 2018a). Here only monthly data with at least 50 % of the daily data available in a month (i.e., at least 5 records) are used. Since station coverage over the central US increases after 2002 (e.g., Pu and Ginoux, 2018a), monthly station data from 2002 to 2015 are used and interpolated to a 0.5 × 0.5^∘ grid using inverse distance weighting interpolation. The gridded data are used to evaluate modeled surface fine dust concentrations.

4. LISA PM₁₀ surface concentration. Surface PM₁₀ concentrations from stations from the Sahelian Dust Transect, which was deployed in 2006 under the framework of the African Monsoon Multidisciplinary Analysis International Program (Marticorena et al., 2010), were used to examine the surface dust concentration over the Sahelian region. The data are maintained by the Laboratoire Interuniversitaire des Systèmes Atmosphériques (LISA) in the framework of the International Network to study Deposition and Atmospheric composition in Africa (INDAAF; Service National d'Observation de l'Institut National des Sciences de l'Univers, France). Three stations are located within the pathway of Saharan and Sahelian dust plumes moving towards the Atlantic Ocean. Here hourly PM₁₀ concentrations from the stations in Banizoumbou (Niger, 13.54^∘ N, 2.66^∘ E), Cinzana (Mali, 13.28^∘ N, 5.93^∘ W), and M'Bour (Senegal, 14.39^∘ N, 16.96^∘ W) from 2006 to 2014 are used. The hourly station data are averaged to obtain daily and monthly mean records to compare with model output.

2.1.4 Other data

Soil depth from the Food and Agriculture Organization of the United Nations (FAO/IIASA/ISRIC/ISS-CAS/JRC, 2009) on a 0.08 × 0.08^∘ resolution is used to examine whether the soil depth is too shallow (i.e., less than 15 cm) for wind erosion.

Figure 1 (a–e) Frequencies of occurrence (FoO; unit: days per season) in each season and the annual mean. (f–j) Threshold of wind erosion (V_threshold; unit: m s⁻¹) derived from satellite products and reanalyses for each season and the annual mean using DOD_thresh=0.2 (or 0.02). Black boxes in (f) denote nine dust source regions as listed in Table 1.

2.2 Retrieving the threshold of wind erosion

The monthly climatological threshold of wind erosion is retrieved by matching the frequency distribution of the MODIS DOD at a certain level, namely, DOD_thresh, with the frequency distribution of surface 10 m winds from the NCEP1 reanalysis over the period from 2003 to 2015. The process can be summarized by the following steps:

• Step 1. Since dust is emitted from the dry and sparsely vegetated surface, the daily DOD data are first masked out to remove the influences of non-erodible factors and unfavorable environmental conditions that are known to prevent dust emission using criteria as follows: daily VSM less than 0.1 cm³ cm⁻³; monthly LAI less than 0.3; monthly snow cover less than 0.2 % (since the snow cover percentage is rounded up to an integer in the MODIS product, this criterion actually requires no snow cover); monthly top-layer soil temperature higher than 273.15 K, i.e., over unfrozen surface; and soil depth thicker than 15 cm. These criteria approximate the most favorable land surface conditions for wind erosion.

  Similar criteria have been used in previous studies to detect or confine dust source regions. For instance, Kim et al. (2013) used the normalized difference vegetation index (NDVI) for less than 0.15, soil depth greater than 10 cm, surface temperature greater than 260 K, and without snow cover to mask topography-based dust source function. LAI less than 0.3 has been used as a threshold for dust emission in the Community Land Model (Mahowald et al., 2010; Kok et al., 2014a), while gravimetric soil moisture ranging from 1.01 % to 11.2 % depending on soil clay content is recommended to constrain dust emission (Fécan et al., 1999). The uncertainties associated with small variations in the retrieval criteria are further quantified and discussed in Sect. 2.3.

• Step 2. Masked daily DOD from Step 1 is then interpolated to a 0.5 × 0.5^∘ grid using bilinear interpolation. This is close to the horizontal resolution of the GFDL AM4.0/LM4.0 model used in this study. Then the cumulative frequency distribution of daily DOD from 2003 to 2015 is derived at each grid point for each month.

• Step 3. Daily maximum surface wind speed is first derived from 6-hourly NCEP1 surface winds and then interpolated to a 0.5 × 0.5^∘ grid. Following Ginoux and Deroubaix (2017), we use maximum daily wind speed instead of daily mean wind speed, largely because dust emission only occurs when wind speed is strong enough, and the emission magnitude is roughly proportional to the third power of surface wind speed in empirical estimations. The cumulative frequency distribution of daily maximum surface wind from 2003 to 2015 is then calculated at each grid point for each month.

• Step 4. A minimum value of DOD (i.e., DOD_thresh) is used to separate dust events from background dust. The cumulative frequency (as a percentage) of dust events passing this threshold is compared to the cumulative frequency of surface winds. The minimum surface winds with the same frequency correspond to the threshold of wind erosion, V_threshold (see a schematic diagram in Fig. S1 in the Supplement). This operation is performed for all grid points for each month. Ginoux et al. (2012) used DOD_thresh=0.2 to separate dust events from background dust and quantify the FoO of local dust events. Similarly, DOD_thresh=0.2 is used here in major dusty regions (northern Africa, the Middle East, India, northern China), while for less dusty regions, such as the US, South America, South Africa, and Australia, DOD_thresh=0.02 is used. The reason to use a lower DOD_thresh for less dusty regions is because (i) the overall dust emission in these regions are at least 10 times smaller than major dusty regions, such as northern Africa (e.g., Huneeus et al., 2011) and (ii) the frequency distribution of DOD in these regions also peaks at a much lower DOD band (see discussion in Sect. 3.3). We also tested the DOD_thresh=0.5 for dusty regions and DOD_thresh=0.05 for less dusty regions, and results are discussed in Sects. 2.3 and 3.1.

  Figure 1a–e show the seasonal and annual mean FoO (days when DOD is greater than DOD_thresh) using DOD_thresh=0.2 or 0.02. The shaded area covers major dust sources, and the pattern is very similar to that obtained by Ginoux et al. (2012; their Fig. 5), although there are some differences, largely due to the masked DOD (i.e., from Step 1) used in this study and a lower threshold in less dusty regions. The higher FoO in northern Africa during summer in comparison with other seasons is consistent with the summer peak of the frequency of the dust source activation derived from the Meteosat Second Generation (MSG) images (Schepanski et al., 2007; their Fig. 1). The relatively high value of FoO over the northern Sahel to southern Sahara is also consistent with dust emission frequency derived from the Meteosat Second Generation Spinning Enhanced Visible and Infrared Imager (Evan et al., 2015; their Fig. 1).

  Note that the selections of masking criteria in Step 1 and DOD_thresh in Step 4 are empirical and can add uncertainties to this method. Also, we approximate dust emission using the cumulative frequency of DOD, which may overestimate dust emission in regions where the contribution of transported dust is significant, and thus we may underestimate the V_threshold in those regions. These uncertainties are further discussed in the following section.

Table 1 Major dust source regions shown in Fig. 1. Note that region names such as India and northern China are not exactly the same as their geographical definitions, as these regions also cover some areas from nearby countries.

Table 2 Sensitivity of the annual mean wind erosion threshold (m s⁻¹) to the selection of different retrieval criteria. Note the setting of the last column is the same as DOD_thresh=0.2 or 0.02, except surface DOD (sDOD) from Aqua is used over northern Africa. Here DOD_thresh=0.2 or 0.5 is applied to dusty regions, i.e., the Sahel, Sahara, Arabian Peninsula, northern China, and India, while DOD_thresh=0.02 or 0.05 is applied to less dusty regions, i.e., the US, South Africa, South America, and Australia.

2.3 Sensitivities of V_threshold to retrieval criteria and the selection of reanalysis surface winds

Table 2 shows variations in derived annual mean V_threshold averaged in nine dust source regions (see Table 1 for locations) following slight changes of retrieval criteria: soil moisture, LAI, snow coverage, and DOD_thresh. When the soil moisture threshold is changed from 0.1 to 0.15 cm³ cm⁻³ or without the soil moisture constraint, the variations in V_threshold are quite small, ranging from 0.01 to about 0.73 m s⁻¹ (Table 2). Similarly, changes of LAI criteria from 0.15 to 0.5 m² m⁻² or snow coverage from 0.2 % to 10 % slightly change V_threshold – within 1 m s⁻¹ over most regions. On the other hand, V_threshold is quite sensitive to the selection of DOD_thresh. V_threshold would increase about 1 to 3 m s⁻¹ if using DOD_thresh=0.5 for dusty regions (0.05 for less dusty regions) instead of DOD_thresh=0.2 (or 0.02). For instance, using DOD_thresh=0.5 increases the averaged annual mean V_threshold over the Sahara from 4.6 m s⁻¹ (using DOD_thresh=0.2) to about 7.6 m s⁻¹.

As mentioned earlier, dust event frequency can be overestimated in regions with a high ratio of transported dust, and consequently V_threshold would be underestimated. Here we provide a rough estimation about the influence of transported dust on V_threshold over northern Africa. It is hard to separate local dust emission and transported dust in the column-integrated DOD, so we use surface DOD data (sDOD, Juliette Paireau, unpublished data), i.e., DOD from the surface to about 400 m, to approximate the component of DOD due to local emission. sDOD is derived by using the DOD vertical profile from the Cloud-Aerosol Lidar with Orthogonal Polarization (CALIOP; Winker et al., 2004., 2007) to first calculate a ratio of near-surface DOD (0–400 m) to total DOD (0–12 km) and then multiplying the ratio with daily MODIS Aqua DOD over northern Africa from 2003 to 2014. Using sDOD, V_threshold over the Sahel would increase from 3.2 to 6.0 m s⁻¹, while over the Sahara, V_threshold would increase from 4.6 to 7.7 m s⁻¹ (Table  2, last column).

Table 3 Sensitivity of the annual mean wind erosion threshold (m s⁻¹) to surface wind speeds from different reanalyses (DOD_thresh=0.2 or 0.02).

How V_threshold would change when using surface winds from different reanalyses is examined in Table 3. Surface winds from the ERA-Interim product produce higher V_threshold values than the NCEP1 product by 0.2 to 2.2 m s⁻¹. Using surface winds from ERA5 also would increase V_threshold by 1 to 1.6 m s⁻¹ over northern Africa and about 1.5 m s⁻¹ over Australia but create smaller differences in other regions.

In short, V_threshold is less sensitive to small changes in the criteria to define a favorable, dry, and sparsely vegetated land surface condition for wind erosion than the choices of DOD_thresh or surface wind speeds from different reanalysis products. Over northern Africa, not separating transported dust from total DOD may lead to an underestimation of V_threshold up to 3 m s⁻¹ based on a rough estimation. However, due to the large uncertainties in quantifying transported dust and the regional converge of the sDOD dataset, we chose not to incorporate the results from sDOD to the global V_threshold.

2.4 Simulation design

We will examine if the observation-constrained, spatially and temporally varying V_threshold would improve dust simulation in GFDL AM4.0/LM4.0. AM4.0/LM4.0 is a coupled land–atmosphere model newly developed at GFDL (Zhao et al., 2018a, b). It uses the recent version of the GFDL Finite-Volume Cubed-Sphere dynamical core (FV3; Putman and Lin, 2007), which was developed for weather and climate applications with both hydrostatic and non-hydrostatic options. Some substantial updates have been incorporated into AM4.0, such as an updated version of the model radiation transfer code, an alternate topographic gravity wave drag formulation, a double-plume model representing shallow and deep convection, a “light” chemistry mechanism, and modulation on aerosol wet removal by convection and frozen precipitation (Zhao et al., 2018a, b). Here we used a model version with 33 vertical levels (with a model top at 1 hPa) and cube sphere with 192×192 grid boxes per cube face (approximately 50 km grid size).

The aerosol physics is based in large part on that of GFDL AM3.0 (Donner et al., 2011), but it has a simplified chemistry where ozone climatology from AM3.0 simulation (Naik et al., 2013) is prescribed. AM4.0 simulates the mass distribution of five aerosols: sulfate, black carbon, organic carbon, dust, and sea salt. Dust is partitioned into 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 (bin 4), and 6–10 µm (bin 5). The dust emission scheme follows the parameterization of Ginoux et al. (2001), as shown in the following equation:

$$\begin{array}{}\text{(4)}& F_p=C\times S\times s_p\times V_{\text{10\hspace{0.17em}m}}^{\mathrm{2}}(V_{\text{10\hspace{0.17em}m}}-V_{\text{t}})(\text{if}{V}_{\text{10\hspace{0.17em}m}}>V_{\text{t}}),\end{array}$$

where F_p is the flux of dust of particle size class p, C is a scaling factor with a unit of µg s² m⁻⁵, here C is set to $\mathrm{0.75}\times {\mathrm{10}}^{-\mathrm{9}}$, S is the source function based on topographic depressions (Ginoux et al., 2001), s_p is fraction of each size class, and V_10 m is the surface 10 m wind speed, and V_t=6 m s⁻¹ is the threshold of wind erosion.

Table 4 Simulation design.

Three simulations with prescribed sea surface temperature (SST) and sea ice (Table 4) were conducted from 1999 to 2015, with the first year discarded for spin up. The Atmospheric Model Intercomparison Project-style (AMIP-style) SST and sea ice data (Taylor et al., 2000) are from the Program for Climate Model Diagnosis and Intercomparison (PCMDI), which combined HadISST (Hadley Centre Global Sea Ice and Sea Surface Temperature; Rayner et al., 2003) from the UK Met Office before 1981 and NCEP Optimum Interpolation (OI) v2 SST (Reynolds et al., 2002) afterwards. The surface winds in the simulations are nudged toward the NCEP1 reanalysis with a relaxation timescale of 6 h (Moorthi and Suarez, 1992). Note that the nudged surface winds are actually weaker than the surface wind speed simulated by the standard version of AM4.0/LM4.0 without nudging, so the overall magnitude of dust emission is lower than the standard version. Here we choose not to retune the dust emission scheme but instead test the usage of V_threshold, which theoretically provides a more physics-based way to improve dust simulation. We also choose to keep the tuning factor C (Eq. 4) the same in all simulations to better examine the effects of implementing the newly developed V_threshold.

In the Control run, the default model setting is used for dust emission, with a prescribed 6 m s⁻¹ threshold of wind erosion (cf. Ginoux et al., 2019). In the V_thresh12 mn simulation, the observation-based climatological monthly V_threshold is used to replace the constant wind erosion threshold. The default source function S in Eq. (4) only allows dust emission over bare ground by masking out regions with vegetation cover. Since LAI masking is already applied in the retrieval of V_threshold (i.e., LAI < 0.3), we choose to use a source function that is the same as the default source function S but without vegetation masking, i.e., S^′ (Fig. S2 in the Supplement). This allows the influence of the spatial and temporal variations in V_threshold to be fully examined. The combination of source function S^′ and V_threshold also extends dust source from bare ground to sparsely vegetated areas as outlined by V_threshold, e.g., over central North America, central India, and part of Australia, and that can increase dust emission in these regions. The pattern of extended dust source area largely resembles the vegetated dust source identified by Ginoux et al. (2012; their Fig. 15b) and Kim et al. (2013; their Fig. 9). All the other settings are the same as the Control run. The V_threshAnn simulation is the same as the V_thresh12 mn but uses the annual mean V_threshold for each month. Since the same SST and sea ice are prescribed for all simulations and land use dose not change much during the short duration of simulation, the differences in simulated dynamic vegetation by LM4.0 among the three simulations are actually very small and can be ignored (see Figs. S3–S4 in the Supplement).

Figure 2 (a–i) Frequency distribution of annual mean V_threshold (black bars) in each region (black boxes in Fig. 1) and V_threshold for dusty seasons, i.e., MAM (green) and JJA (blue) for regions in the Northern Hemisphere and SON (orange) and DJF (grey) for regions in the Southern Hemisphere. The mean (averaged over all grid points in the region, without area weight) and ± 1 standard deviation of V_threshold in each region are shown on the top right of each plot.

3 Results

3.1 Thresholds of wind erosion with DOD_thresh=0.2 (or 0.02) and DOD_thresh=0.5 (or 0.05)

Figure 1f–j show the derived threshold of wind erosion for each season and annual mean using DOD_thresh=0.2 (or 0.02). The seasonal variations in the wind erosion threshold are largely due to the variations in DOD and surface wind frequency distributions that are in turn associated with variations in land surface features, such as soil moisture, soil temperature, snow cover, and vegetation coverage in each month. V_threshold is generally lower in MAM (March–April–May) and JJA (June–July–August) (SON (September–October–November) and DJF (December–January–February)) for Northern (Southern) Hemisphere dusty regions than in other seasons, consistent with higher FoO in these seasons. V_threshold values are also lower in regions with a high FoO (Fig. 1a–e).

The distributions of V_threshold for the annual mean (black bars) and dusty seasons (color lines; MAM and JJA for the Northern Hemisphere and SON and DJF for the Southern Hemisphere) for each dust source region (see Fig. 1f and Table 1 for locations) are shown in Fig. 2a–i. In the Sahel and Sahara, the annual mean V_threshold peaks around 4 and 4.5–5.5 m s⁻¹, respectively (Fig. 2a–b). This magnitude is lower than indicated from previous studies based on station observations in the region, e.g., Helgren and Prospero (1987) found the threshold velocity over eight stations in northwestern Africa ranged from 6.5 to 13 m s⁻¹ during summer in 1974. Chomette et al. (1999) and Marsham et al. (2013) also reported higher wind erosion thresholds around 6–9 m s⁻¹ at individual stations. On the other hand, Cowie et al. (2014) found that the annual threshold of wind erosion at the 25 % level, i.e., when surface condition is favorable for dust emission, can be lower than 6 m s⁻¹ at some sites in the Sahel (their Fig. 5). Several factors may contribute to the discrepancies. Firstly, studies suggest that reanalysis datasets may underestimate surface wind speed in spring and for monsoon days in Africa (e.g., Largeron et al., 2015), and therefore they could lead to a lower value of V_threshold than that derived from station observations. In fact, Bergametti et al. (2017) found even 3-hourly wind speed records at stations may miss short events with a high wind speed. As shown in Table 3, among the reanalysis wind products tested here, NCEP1 actually produced a lower V_threshold in northern Africa than the other two reanalyses. Secondly, using the DOD frequency to approximate dust emission may lead to an overestimation of dust emission over regions such as the southern Sahel where transported dust is a large component and consequently an underestimation of V_threshold. Based on our rough estimation, V_threshold in northern Africa can be underestimated by up to 3 m s⁻¹ (Sect. 2.3). In addition, different analysis time periods or methods to retrieve the wind erosion threshold may also contribute to the differences.

Figure 3 (a–e) Threshold of wind erosion (V_threshold; unit: m s⁻¹) derived from satellite products and reanalyses for each season and the annual mean using DOD_thresh=0.5 (or 0.05). Black boxes in (a) denote nine dust source regions as listed in Table 1.

The annual mean V_threshold in the Arabian Peninsula is a bit higher, with mean values at 5.2 m s⁻¹ (Fig. 2c). The V_threshold over northern China is even higher, with an annual mean of 7.8 m s⁻¹. This is consistent with the results of Kurosaki and Mikami (2007), who found that under favorable land surface conditions the threshold wind speed ranges from 4.4±0.6 m s⁻¹ in the Taklimakan Desert to 6.9±1.2 m s⁻¹ over the Loess Plateau and around 9.8±1.6 m s⁻¹ in the Gobi Desert. These values are also consistent with Ginoux and Deroubaix (2017) who found that the regional mean wind erosion threshold over northern China ranges from 6.5 to 9.1 m s⁻¹. In India, the V_threshold peaks at about 4.5 m s⁻¹ and 6.5 m s⁻¹, respectively (Fig. 2e). The second peak is probably related to anthropogenic dust sources over the central Indian subcontinent (Ginoux et al., 2012). We also note that in the Northern Hemisphere, V_threshold in dusty seasons is shifted towards lower values than the annual mean (blue and green lines in Fig. 2a–f), but it is similar to the annual mean in the Southern Hemisphere (especially South America and Australia), indicating stronger influences of surface variability in the Northern Hemisphere.

Figure 3 shows the seasonal mean and annual mean global V_threshold using DOD_thresh=0.5 (or 0.05). The corresponding distribution of annual mean V_threshold in each region is shown in Fig. S5 in the Supplement. The derived V_threshold is generally higher than using DOD_thresh=0.2 (or 0.02), especially over northern Africa, the Arabian Peninsula, India, and Asia (Fig. 3 and Table 2). The results are thus closer to previous station-based studies over northern Africa. On the other hand, over northern China, V_threshold is around or greater than 8 m s⁻¹ (Fig. 3e), slighter higher than previous estimates (e.g., Kurosaki and Mikami, 2007; Ginoux and Deroubaix, 2017).

In the following section, we will exam if the spatially and temporally varying V_threshold would improve model simulation of the DOD spatial pattern, seasonal variations, frequency distribution, and surface dust concentrations in the GFDL AM4.0/LM4.0. Results using V_threshold with DOD_thresh=0.2 (or 0.02) are shown in Sect. 3.2–3.3, and results using V_threshold with DOD_thresh=0.5 (or 0.05) are briefly discussed in Sect. 4.

3.2 V_threshold in the GFDL AM4.0/LM4.0 model

In this section we analyze the model output using the default setting (Control; Table 4), 12-month (V_thresh12 mn), and annual mean V_threshold (V_threshAnn) by comparing model results with multiple observational datasets and MODIS DOD.

Figure 4 Climatology of annual mean AERONET (a) AOD (550 nm) and (b) SDA COD (500 nm) averaged over 2003–2015.

Figure 5 Scatter plot of simulated annual mean (a) AOD and (b) COD in the Control run versus AERONET AOD and COD (a, b). The relative difference (as a percentage) (c) between modeled AOD and AERONET AOD and (d) between modeled COD and AERONET COD (c–e). (e) The relative contribution of DOD to COD in the model.

3.2.1 Climatology of AOD and DOD

In order to compare the model results with observations, we first show the climatology of AERONET AOD and COD from 2003 to 2015. The length of records for each station is shown in Fig. S6 in the Supplement. As shown in Fig. 4, annual mean global AOD is highest over Africa, the Arabian Peninsula, the Indian subcontinent, and southeastern Asia. In the latter two regions, high sulfate concentrations (e.g., Ginoux et al., 2006) and organic carbon from biomass burning in southeastern Asia (e.g., Lin et al., 2014) contribute substantially to the total AOD. The SDA COD shows the optical depth due to coarse aerosols, which includes both dust and sea salt, and sea salt over coastal regions or islands can be a major contributor. Here, high values (> 0.2) are largely located over dusty regions such as northern Africa, the Arabian Peninsula, and northern India (Fig. 4b).

Figure 5a–b show the scatter plots of modeled AOD and COD in the Control run versus AERONET AOD and COD, respectively. Here column-integrated extinction from both dust and sea salt is used to calculated COD in the model. The relative differences (as a percentage) between AM4.0 output and AERONET station data are also shown (Fig. 5c–d). The percentage of DOD to total COD in the model is displayed at the bottom (Fig. 5e). The simulated AOD is lower than that from AERONET over northern Africa, the Middle East, and western India, largely due to low values of COD simulated in these regions (Fig. 5d). Besides these regions, the COD over North America, South America, South Africa, and northern Eurasia is also, for the most part, underestimated by the model. Dust is the dominant contributor to the COD value over most of these low COD regions, except over central-to-eastern North America and central South America (Fig. 5e).

Figure 6 Same as Fig. 5 but for the V_thresh12 mn simulation.

COD (and effectively DOD given its dominance in most regions) was better simulated in the subsequent model run using a prescribed 12-month V_threshold in terms of both magnitude and spatial pattern. Figure 6 shows the results from the V_thresh12 mn simulation. COD is better captured, while the AOD effectively moves from a negative to a slightly positive bias (Fig. 6a–d). Most sites over northern Africa and the Middle East show a relatively small difference with AERONET COD (Fig. 6d). Over the Indian subcontinent, COD is overestimated, while over North America, excluding the east coast, northern Eurasia, and part of South America, COD is also better captured than in the Control run.

These improvements are largely associated with a better simulation of DOD in the “dust belt” (i.e., northern Africa and the Middle East). Figure 7 shows the DOD at 550 nm derived from AERONET AOD (see methodology for details) versus that from the V_thresh12 mn simulation. Over most stations in the Sahel, Mediterranean coasts, and central Middle East, the relative differences between modeled and observed DOD is within ±25 %.

Figure 7 (a) Climatology (2003–2015) of AERONET DOD (550 nm) over major dusty regions and (b) scatter plot of modeled DOD in the V_thresh12 mn simulation versus AERONET DOD. (c) The relative difference (as a percentage) between modeled DOD and AERONET DOD in the V_thresh12 mn simulation.

Figure 8 shows the regionally averaged annual mean DOD over nine dusty regions from MODIS and three simulations. The Control run largely underestimates DOD in all regions, while the magnitude of DOD is better captured in the V_thresh12 mn and V_threshAnn simulations, although slightly overestimated in the Sahel and greatly overestimated over Australia. In general, DOD simulated by the V_threshAnn run using a constant annual mean V_threshold is higher than that simulated by the V_thresh12 mn run, consistent with the higher dust emission in the V_threshAnn run (Table S2 in the Supplement). Lack of the soil moisture constraint in the model, which is a very important element in capturing the variation of DOD in Australia (Evans et al., 2016), may contribute to the large overestimation of DOD in Australia.

3.2.2 Climatology of surface dust concentration

While DOD is a key parameter associated with the climate impact of dust, surface dust concentration is an important factor affecting local air quality. Here we compare the modeled surface dust concentration with RSMAS station observations. Model output is averaged from 2000 to 2015 to form the annual climatology. Consistent with the DOD output, the Control run largely underestimates surface dust concentrations at almost all of the sites (except sites 9 and 15; Fig. 9, top panel). The underestimation is reduced in the V_threshAnn simulation (Fig. 9, middle panel), with seven stations having model/observation ratios between 0.5 and 2 (white triangles). Over the coastal US (e.g., sites 16 and 13), dust concentrations are overestimated, consistent with the overestimation of DOD over the US and the Sahel (Fig. 8). Dust concentrations in Australia and the east coast of China are also overestimated by more than fivefold. Surface dust concentration is further improved in the V_thresh12 mn simulation (Fig. 9, bottom), with eight stations showing a model/observation ratio between 0.5 and 2 and only four stations overestimating or underestimating dust concentrations by more than 5 times.

Figure 8 Regionally averaged annual mean DOD (2003–2015) over nine regions from the Control (grey), V_thresh12 mn (orange), and V_threshAnn (yellow) simulations and MODIS (black).

Figure 9 Scatter plots (a, c, e) of model simulated (from top to bottom are the Control, V_threshAnn, and V_thresh12 mn simulations) surface dust concentration (µg m⁻³) versus the climatology of observed surface dust concentration from RSMAS stations (Savoie and Prospero, 1989). Spatial pattern of surface dust concentration from model output (shading; b, d, f) and the ratio between modeled and RSMAS station observed surface dust concentration (color triangles, with upward triangles indicating overestimation and downward triangles indicating underestimation). A total of 16 stations were used, and the numbers in each triangle (b, d, f) and grey dots (a, c, e) indicate the stations. The one–one, one–two, and one–five lines are plotted in solid, dashed, and dash-dotted lines in the scatter plots. Statistics in the scatter plots are calculated in logarithmic space.

Simulated surface fine dust concentration (calculated as dust bin $\mathrm{1}+\mathrm{0.25}\times$ dust bin 2) in the US is compared with gridded IMPROVE data (Fig. 10). While the Control run largely underestimates surface fine dust concentration, the simulated concentration is overall too high in the V_threshAnn run. The spatial pattern of fine dust concentration is better captured in the V_thresh12 mn run, with higher values over the southwestern US, but the magnitude is still overestimated, and additional dust hot spots are simulated over the northern Great Plains and the Midwest, which are not shown in the IMPROVE data. Such an overall overestimation may be attributed to a lack of soil moisture modulation in the dust emission scheme. The way in which dust bins are partitioned in the model can add uncertainties to model's representation of surface fine dust concentrations as well. On the other hand, the relatively low spatial coverage of IMPROVE sites over the northern Great Plains and Midwest (e.g., Pu and Ginoux, 2018a) may also add uncertainties to the data itself.

3.2.3 Seasonal cycles

Figure 11 compares the seasonal cycle of DOD from three simulations with MODIS DOD in nine dusty regions. The seasonal cycle of gridded AERONET COD (as an approximation of DOD; on a 0.5 × 0.5^∘ grid) is also shown. Since the gridded COD may have large uncertainties over regions with only a few stations, such as the Sahel, Sahara, northern China, and South Africa, MODIS DOD is used as the main reference in the comparison. Seasonal cycles are better captured by the V_thresh12 mn simulation in the Sahel, the Sahara, and the Arabian Peninsula (Fig. 11a–c), although the spring and summer peak in the Sahel is overestimated, and the winter minimum in the Sahara is underestimated. The MAM peak of MODIS DOD in northern China is missed by both V_thresh12 mn and V_threshAnn simulations (Fig. 11d), while the JJA peak over India is largely overestimated (Fig. 11e). Over the US dusty region, the seasonal cycle in the V_thresh12 mn simulation is slightly underestimated compared to MODIS DOD but overestimated from May to August in the V_threshAnn simulation (Fig. 11f). DOD is underestimated in South Africa in all three simulations (Fig. 11g). Over South America, the peak from October to February is roughly captured by the V_thresh12 mn run but is overestimated by the V_threshAnn run (Fig. 11h). The seasonal cycles of DOD in Australia are very similar in all three simulations and largely resemble that in MODIS, although both the V_thresh12 mn and V_threshAnn simulations overestimate the DOD by about an order of magnitude.

Figure 12 shows the seasonal cycle of COD from 12 AERONET SDA sites over northern Africa and nearby islands (see Fig. S7 in the Supplement for site locations) along with MODIS DOD and DOD simulated in three runs. The magnitude of AERONET COD and MODIS DOD in these sites are very similar, despite missing values at sites 1, 4, 5, 8, and 11 and a smaller value at site 2 in MODIS. Over most of the sites, the seasonal cycle is better captured in the V_thresh12 mn and V_threshAnn simulations than the Control run, although the peak over Cairo_EMA_2 (site 12) is slightly underestimated, which is consistent with the underestimation of annual mean DOD in the area (Fig. 7).

We also examined the seasonal cycle of PM₁₀ surface concentration at three Sahelian INDAAF stations (see Fig. S7 in the Supplement for site locations) from the LISA project. Figure 13a–c show PM₁₀ surface dust concentration (here dust dominates total PM₁₀ concentration) from the Control, V_thresh12 mn, and V_threshAnn simulations versus observed PM₁₀ concentrations from three LISA sites. PM₁₀ concentrations in these sites peak during boreal winter and spring and reach minima from July to September. These seasonal variations are associated with the dry northerly Harmattan wind in boreal winter and spring that transports Saharan dust southward to the Guinean coast and the scavenging effect of monsoonal rainfall in boreal summer that removes surface dust (Marticorena et al., 2010; Fiedler et al., 2015). While the Control run does not capture the seasonal cycles in these sites, the V_thresh12 mn run largely captures the spring peak and summer minimum, although the magnitude is overestimated. In all three sites, the simulated concentration in the V_threshAnn run is larger than that in the V_thresh12 mn run, especially from boreal fall to early spring. Such an overestimation is probably due to the prescribed constant annual mean V_threshold, which is lower than it would be during the less dusty season (i.e., boreal fall to winter) and thus increases dust emission and surface concentrations.

Figure 13d–f show the seasonal cycle of DOD from three AERONET sites co-located with LISA INDAAF stations and from three simulations. The V_thresh12 mn and V_threshAnn simulations largely captured the seasonal cycle of DOD at these sites. The overestimation of near-surface PM₁₀ dust concentrations (Fig. 13a–c) and the generally well-captured column-integrated DOD (Fig. 13d–f) indicate that the model likely underestimates dust concentration in the atmospheric column above the surface, which needs further investigation in future studies.

3.2.4 A dust storm over the US northern Great Plains on 18 October 2012

Can the AM4.0/LM4.0 with the prescribed V_threshold better represent individual dust events? Here we examine a major dust storm captured by a MODIS Aqua true color image on 18 October 2012 (https://earthobservatory.nasa.gov/images/79459/dust-storm-in-the-great-plains, last access: November 2018) over the US northern Great Plains. There was a severe drought in 2012 with anomalously low precipitation centered over the central US (e.g., Hoerling et al., 2014). The dry conditions favored dust storm development when there were intensified surface winds. However, this storm was not predicted by the forecast models, such as the Goddard Earth Observing System version 5 (GEOS-5; Rienecker et al., 2008) and Navy Aerosol Analysis Prediction System (NAAPS; Witek et al., 2007; Reid et al., 2009; Westphal et al., 2009).

Figure 10 Annual mean surface fine dust concentration (µg m⁻³) from IMPROVE stations (left column), three simulations (middle column), and the differences between model and observation (right column) for 2002–2015.

Figure 11 Seasonal cycle of DOD from MODIS (black), the Control (grey), V_thresh12 mn (orange), and V_threshAnn (yellow) runs, with gridded AERONET SDA COD (blue) averaged over nine regions. The annual mean of each dataset in each region is listed on the top of the plot. Note that in (i) DOD from the V_thresh12 mn and V_threshAnn simulations is multiplied by 0.1 for clarity.

Figure 12 Seasonal cycle of DOD over 12 AERONET SDA sites (see Fig. S7 in the Supplement for locations) from the Control (grey), V_thresh12 mn (orange), and V_threshAnn (yellow) simulations, along with DOD from MODIS (blue) and COD from AERONET (black dotted line). All values are averaged over 2003–2015. The location (coordinates) and the name (due to space, only the first seven characters are shown) of the sites are listed at the top of each plot.

Figure 13 (a–c) Seasonal cycle of PM₁₀ surface concentration (black) over three sites from the LISA project, along with PM₁₀ surface dust concentration from the Control (grey), V_thresh12 mn (orange), and V_threshAnn (yellow) simulations. Error bars are ± 1 standard deviation of the daily mean in each month averaged over 2006–2014 (units: µg m⁻³). (d–f) Seasonal cycle of DOD (550 nm) from three AERONET sites co-located with LISA sites (blue) versus that modeled by the Control (grey), V_thresh12 mn (orange), and V_threshAnn (yellow) simulations.

Figure 14 (a–c) Daily DOD from MODIS. (d–f) Daily DOD simulated by the V_thresh12 mn run along with anomalies (with reference to the 2000–2015 mean) of surface wind vectors (m s⁻¹) from 17 October to 19 October 2012. Only DOD over land is shown. Missing values in MODIS DOD (a–c) are plotted in grey shading.

As shown in Fig. 14, MODIS DOD captures this event, with a peak value above 0.5 over southwest Nebraska and northern Kansas on 18 October 2012. The V_thresh12 mn run also largely captures this event (Fig. 14, bottom panel), although the Control run totally misses it (not shown). In the model, the dust storm appears in South Dakota and Nebraska on 17 October 2012, along with the anomalous southwesterly winds. It reaches a maximum on 18 October in association with intensified anomalous southwesterly winds at the surface and an anomalous low-pressure system at 850 hPa (Fig. S8 in the Supplement). Note that the modeled dust storm center is located a bit northeastward compared to the MODIS DOD pattern and it also has greater magnitude and covers a larger area. On 19 October both the anomalous low-pressure system and surface wind speeds weaken, and the dust storm dissipates, with slightly elevated DOD levels over a region extending over the lower Mississippi River basin and the Midwest. This is somewhat consistent with MODIS records, which also shows slightly higher DOD levels over Tennessee and northern Alabama on 19 October regardless of the large area of missing values.

Figure 15 Frequency distribution (as a percentage) of regionally averaged daily DOD from MODIS (black) versus that from the Control (light blue) and V_thresh12 mn (orange) simulations for the Sahara, the Sahel, the Arabian Peninsula, northern China, India, the western-to-central US, South America, South Africa, and Australia from 2003 to 2015. The x axis denotes the ranges of DOD (the bin spacing for dusty regions is 0.05 and for less dusty regions is 0.01), and the y axis is the percentage of occurrence. The light-green boxes denote the averaging areas. For regions in the Northern Hemisphere frequency in MAM is shown, while for regions in the Southern Hemisphere frequency in SON is shown.

3.3 Frequency distribution of DOD in the model versus that from MODIS

Figure 15 shows the frequency distribution of regional mean DOD during one dusty season (MAM in the Northern Hemisphere and SON in the Southern Hemisphere) for nine regions. Results from MODIS, the Control, and V_thresh12 mn runs are shown in black, blue, and orange lines, respectively. In most dusty regions, such as the Sahara, the Sahel, the Arabian Peninsula, India, and northern China, MODIS DOD frequency largely peaks between 0.2 to 0.4, while DOD frequency peaks at a much lower level between 0.02 to 0.08 in less dusty regions, such as the US, South America, South Africa, and Australia. This also justifies our selection of DOD_thresh of 0.02 (instead of 0.2) in the less dusty regions. The DOD distribution in the Control run is biased low and peaks around 0.05 in those dusty regions and between 0 and 0.01 in less dusty regions. The frequency is much better captured in the V_thresh12 mn run over the Arabian Peninsula and the Sahel, slightly improved but still biased low over the Sahara, northern China, India, and the US. The modeled frequency in the V_thresh12 mn run is biased high in Australia (peaks outside the maximum of the x axis, not shown) and shows little improvement over South Africa and South America. The overall improvement of DOD frequency using the time-varying 2-D V_threshold occurs mostly over major dusty regions, which is consistent with the improvements in DOD climatology and seasonal cycle in the model simulations.

4 Discussion

A global distribution of the threshold of wind erosion is retrieved using high-resolution MODIS DOD and land surface constraints from relatively high-resolution satellite products and reanalyses. While this climatological monthly V_threshold provides useful information about the spatial and temporal variations in wind erosion threshold, there are some uncertainties associated with it. Here DOD frequency is derived using MODIS and other satellite products; thus the uncertainties in the satellite products are inherited in the derived DOD frequency distribution. Due to the cloud screening processes of MODIS products, dust activities over cloud-covered regions may be underestimated. Also, DOD frequency is derived based on daily observations over a 13-year record so that some variability of dust emission associated with alluvial sediments deposited by seasonal flooding may be not captured. Diurnal variability of dust emission and short-duration events such as haboobs are also not included. Since DOD is a column-integrated variable, it includes both locally emitted and remotely transported dust. When using DOD frequency distribution to approximate dust emission, it may overestimate dust emission in regions where transported dust is dominated and lead to an underestimation of V_threshold. Future studies to better quantify the influences of transported dust are needed to further improve the quantitative retrieval of V_threshold.

Previous studies found that over regions such as northern Africa, reanalysis products may underestimate surface wind speed in spring and monsoon seasons but overestimate it during dry nights (e.g., Largeron et al., 2015). This is largely because mechanisms such as density current that can enhance surface wind speed are not parameterized in the atmospheric models to produce the reanalysis products, while coarse spatial and temporal sampling may also contribute to the underestimation of reanalysis wind speeds. The selection of surface winds from different reanalysis products also affects the derived V_threshold. Among the three reanalyses examined here, V_threshold, derived from the NCEP1 reanalysis, shows slightly lower values than others.

In addition, V_threshold is derived by matching the frequency distribution of DOD at certain levels (i.e., DOD_thresh) with the frequency distribution of daily maximum wind. An issue is that selecting a value of DOD_thresh is quite empirical. The influences of soil properties such as soil cohesion, particle size, and particle compositions on the threshold of wind erosion (e.g., Fécan et al., 1999; Alfaro and Gomes, 2001; Shao, 2001; Kok et al., 2014b) are not explicitly examined here and will need further investigation.

The influences of V_threshold on AM4.0/LM4.0 results are twofold. On the one hand, it modifies the default constant threshold of wind erosion (V_t in Eq. 4) by allowing spatial and temporal variations of wind erosion threshold over bare ground, i.e., within the domain of default dust source function S (Fig. S9a–e in the Supplement). On the other hand, it slightly extends the potential emission area to sparsely vegetated regions as outlined by V_threshold (Fig. S9f–j in the Supplement). Which effect dominates? Taking the V_thresh12 mn simulation as an example, Fig. S10 shows the differences of dust emission with the Control run. The increase of dust emission in the V_thresh12 mn simulation (also summarized in Table S2 in the Supplement) is largely associated with the enhanced emission over the bare ground (Fig. S10a–e in the Supplement), mainly over the regions with a reduced wind erosion threshold (Fig. S9a–e in the Supplement). The increased emission over sparsely vegetated areas over regions such as the southern Sahel, India, and Australia plays a minor role. This is consistent with Kim et al. (2013), who found that global dust emission in the Georgia Institute of Technology–Goddard Ozone Chemistry Aerosol Radiation and Transport (GOCART) model is dominated by emission from the bare ground.

The major benefit of using the spatially and temporally varying V_threshold is that it improves the simulation of the DOD spatial pattern (Figs. 6–7), seasonal cycle (Figs. 11–13), and frequency distribution (Fig. 15) as well as the spatial pattern of surface dust concentrations (Figs. 9–10), which cannot be achieved by simply modifying the global tuning factor (i.e., C in Eq. 4) to fit the observations such as surface concentrations or optical depth.

The default setting in the Control run produced relatively low global dust emissions (978 Tg yr⁻¹) in comparison with the Aerosol Comparisons between Observations and Models (AeroCom) multi-model median (1123 Tg yr⁻¹; Huneeus et al., 2011) or a previous estimation based on MODIS DOD (1223 Tg yr⁻¹; Ginoux et al., 2012). So we also conducted a test run (Control II) to increase global dust emission in the Control run to about 1232 Tg yr⁻¹ by enlarging C in Eq. (4). The magnitude of DOD slightly increases, e.g., over the Sahel annual mean increases from 0.07 to 0.09; however, there is no improvement in terms of the seasonal cycle or spatial pattern, as expected.

We also examined the performance of V_threshold using DOD_thresh=0.5 (or 0.05) in the AM4.0/LM4.0. Similarly, we conducted simulations with 12-month V_threshold (V_thresh12 mn II) and annual mean V_threshold (V_threshAnn II), all using the same tuning factor as in Control II. We found similar improvement in the DOD seasonal cycle and weaker improvement in the DOD spatial pattern and frequency distribution and surface dust concentrations (except with the IMPROVE data over the US and surface concentrations over the Sahel, where dust concentrations are previously overestimated). This is largely because higher V_threshold leads to lower global dust emissions in the V_threshAnn II (1961 Tg yr⁻¹) and V_thresh12 mn II simulations (1705 Tg yr⁻¹) and overall lower DOD. Over the Mediterranean coast, Europe, and northern Asia, the DOD spatial pattern is not as well captured in the V_thresh12 mn II run as in the V_threshold12 mn run, likely due to relatively high V_threshold in these regions.

5 Conclusion

While dust aerosols play an important role in the Earth's climate system, large uncertainties exist in modeling its life cycle (e.g., Huneeus et al., 2011; Pu and Ginoux, 2018b). Constant thresholds of wind erosion are widely used in climate models for simplicity. Here, high-resolution MODIS Deep Blue dust optical depth and surface wind speeds from the NCEP1 reanalysis, along with other land surface factors that affect wind erosion, such as soil moisture, vegetation cover, snow cover, soil temperature, and soil depth, were used to develop a time-varying two-dimensional climatological threshold of wind erosion, V_threshold, based on the seasonal variations of DOD and surface wind distribution frequencies. V_threshold is generally lower in dusty seasons, i.e., MAM and JJA (SON and DJF) in the Northern (Southern) Hemisphere.

The climatological monthly V_threshold was then incorporated into the GFDL AM4.0/LM4.0 model to examine the potential benefits relative to the use of a constant threshold. In comparison with the simulation using the default setting of a globally constant threshold of wind erosion (6 m s⁻¹), both the magnitude of DOD and surface dust concentrations are increased and closer to observations. However, different from modifying the global tuning factor (i.e., C in Eq. 4) to increase the overall magnitudes of DOD or surface dust concentrations, we found the spatially and temporally varying V_threshold largely improves the simulation of the spatial pattern, seasonal cycle, and frequency distribution of DOD over Northern Hemisphere dusty regions, such as northern Africa and the Arabian Peninsula, and slightly improves the simulation over India, the western-to-central US, and northern China. The seasonal cycle of DOD is also slightly improved in South America, although it changes little in South Africa. The incorporation of V_threshold leads to an overestimation of DOD in Australia, likely in association with the absence of soil moisture constraints on dust emission in the model.

The spatial pattern of surface dust concentrations is also improved when spatially and temporally varying V_threshold is incorporated. The fine dust concentration in the US is also better captured, with the maximum of the annual mean largely located over the southwestern US, although the magnitude is overestimated.

A constant annual mean V_threshold is also tested in the model, and it is found to overestimate DOD over dusty seasons in the Arabian Peninsula, US, India, Australia, and South America. Surface PM₁₀ concentrations in the Sahel during boreal fall and winter seasons are also largely overestimated with this setting. The results indicate the importance of including the seasonal cycle of V_threshold in the model. Using time-varying 2-D V_threshold, the model was also able to capture a strong dust storm in the US Great Plains in October 2012, which created deadly accidents, while some dust forecasting models failed to reproduce it.

Finally, this method to retrieve the global threshold of wind erosion can be conducted under different resolutions or surface wind reanalyses or applied to surface friction velocity datasets to match the resolution/scheme of dust models and may help improve their simulations and forecasting of dust distribution. As discussed in Sect. 4, there are uncertainties associated with this method, and future studies to better quantify the influence of transported dust to the overall DOD frequency distribution and incorporate station-based surface wind records into the retrieval process will further improve the dataset.