A revised mineral dust emission scheme in GEOS-Chem: improvements in dust simulations over China

Mineral dust plays a significant role in climate change and air quality, but large uncertainties remain in terms of dust emission prediction. In this study, we improved treatment of the dust emission process in a global 3-D chemical transport model (GEOS-Chem v12.6.0), by incorporating the geographical variation of aerodynamic roughness length (Z₀), smooth roughness length (Z_0s) and soil texture and by introducing the Owen effect and the formulation of the sandblasting efficiency α by Lu and Shao (1999). To investigate the impact of the modifications incorporated in the model, several sensitivity simulations were performed for a severe dust storm during 27 March to 2 April 2015 over northern China. Results show that simulated threshold friction velocity is very sensitive to the updated Z₀ and Z_0s field, with the relative difference ranging from 10 % to 60 % compared to the original model with a uniform value. The inclusion of the Owen effect leads to an increase in surface friction velocity, which mainly occurs in the arid and semi-arid regions of northwest China. The substitution of a fixed value of α assumed in the original scheme with one varying with friction velocity and soil texture based on observations reduces α by 50 % on average, especially over regions with sand texture. Comparisons of sensitivity simulations and measurements show that the revised scheme with the implementation of updates provides more realistic threshold friction velocities and PM₁₀ mass concentrations. The performance of the improved model has been evaluated against surface PM₁₀ observations as well as MODIS aerosol optical depth (AOD) values, showing that the spatial and temporal variation of mineral dust are better captured by the revised scheme. Due to the inclusion of the improvement, average PM₁₀ concentrations at observational sites are more comparable to the observations, and the average mean bias (MB) and normalized mean bias (NMB) values are reduced from −196.29 µg m⁻³ and −52.79 % to −47.72 µg m⁻³ and −22.46 % respectively. Our study suggests that the erodibility factor, sandblasting efficiency and soil-related properties which are simply assumed in the empirical scheme may lack a physical mechanism and spatial–temporal representativeness. Further study and measurements should be conducted to obtain a more realistic and detailed map of these parameters in order to improve dust representation in the model.

1 Introduction

Mineral dust is typically produced by wind erosion from regions with arid and semi-arid surfaces in the world and exerts significant impacts on the atmospheric radiation balance (Tegen et al., 1996; DeMott et al., 2010; Kumar et al., 2014; Saidou Chaibou et al., 2020a), climate (DeMott et al., 2003; Mahowald and Kiehl, 2003; Zhao et al., 2012; Chen et al., 2014; Chin et al., 2014), air quality (Giannadaki et al., 2014; Tian et al., 2019) and human health (Goudie, 2014; Tong et al., 2017). The dust emission process has been recognized as a leading contributor to dust aerosol loading. Global mineral dust particles are mainly emitted from North Africa, the Arabian Peninsula, Central Asia, East Asia, Australia and North America, with East Asia (including the Gobi and Taklimakan deserts) accounting for ∼ 20 % of the global dust emission (Ginoux et al., 2004; Nagashima et al., 2016).

In order to properly reproduce the dust emission process, many dust emission schemes have been developed and implemented in both global and regional chemical transport models (CTMs) (e.g., Marticorena and Bergametti, 1995; Lu and Shao, 1999; Alfaro and Gomes, 2001; Shao, 2001, 2004; Shao et al., 2011; Zender et al., 2003; Kok, 2011a, b). Nevertheless, some intercomparison studies demonstrated that there are large discrepancies among different dust emission models (Uno et al., 2006; Todd et al., 2008; Huneeus et al., 2011; Su and Fung, 2015; Ridley et al., 2016; Chen et al., 2017; Ma et al., 2019; Wu et al., 2019; Saidou Chaibou et al., 2020b; Zhao et al., 2020). Ma et al. (2019) quantitatively evaluated the performance of three dust schemes in the Weather Research and Forecasting model coupled with Chemistry (WRF-Chem), two schemes in both CHIMERE and CMAQ and one scheme in CAMx during a dust episode over northern China. Large differences between observed surface PM₁₀ concentrations and the modeling results of each model were found. Among schemes in WRF-Chem, AFWA and UOC_Shao2004 are better correlated with observations compared to GOCART but tend to overestimate dust concentrations. Kang et al. (2011) compared the performance of three dust emission schemes in WRF-Chem over East Asia, showing that the difference of dust emission fluxes between the three schemes ranges from an order of 10¹ to 10². Ridley et al. (2016) showed that the estimated global dust aerosol optical depth (AOD) varies by over a factor of 5 among four global models (including GEOS-Chem, WRF-Chem, CESM and MERRAero), and dust emissions across North Africa are overestimated, while emissions from Asia and the Middle East are underestimated overall. An intercomparison of 14 CTMs as part of the Model Inter-Comparison Study for Asia (MICS-Asia) phase III project (Chen et al., 2019) showed that nearly all participant models underestimate PM₁₀ levels and that current CTMs have difficulty producing similar dust emissions when adopting different dust schemes.

The uncertainties in dust emission models can be attributed to a number of issues, such as threshold friction velocity, surface wind speed, soil texture, particle size distribution, other soil/surface parameters (e.g., soil moisture, vegetation cover and aerodynamic roughness length) and different physical mechanisms (Tegen, 2003; Zhao et al., 2013; Liu et al., 2018; Chen et al., 2019). Darmenova et al. (2009) conducted a detailed comparison between two schemes developed by Marticorena and Bergametti (1995) and Shao et al. (1996), indicating that wind friction velocity is a significant factor in simulating dust emission, while the aerodynamic roughness length as well as vegetation cover may play an important role at higher wind speed. Many sensitivity experiments have been conducted and have shown that the modeled threshold friction velocity can be modified by soil moisture (Cheng et al., 2008; Mokhtari et al., 2012; Gherboudj et al., 2015; Ju et al., 2018), soil texture (Menut et al., 2013; Gherboudj et al., 2015; Perlwitz et al., 2015a, b; Kontos et al., 2018) and surface roughness (Cheng et al., 2008; Astitha et al., 2012; Menut et al., 2013), which in turn affects the predicted dust emission. In addition, a more accurate value of the sandblasting mass efficiency (α) has been reported to be a crucial factor for a better performance of the dust emission flux (Mokhtari et al., 2012; Klingmüller et al., 2018; Kontos et al., 2018; Ma et al., 2019).

Based on the above studies, it is necessary to take these key parameters, including soil-related properties and empirical input parameters, into full consideration in a dust emission parameterization. Unfortunately, due to limited observations, many of these parameters are not included well in the model. For example, most dust models simply assume constant values of aerodynamic roughness length and soil clay fraction (Ginoux et al., 2001; Tegen et al., 2002; Zender et al., 2003), ignoring their temporal and spatial variability, which may cause uncertainties in the estimated surface friction velocity and threshold friction velocity. During recent decades, with the development of observation technology, detailed information on the surface characteristics appropriate for global and regional models has been provided (Laurent et al., 2005, 2008; Prigent et al., 2005, 2012; Shangguan et al., 2014; Perlwitz et al., 2015a, b). Therefore, adopting more accurate and detailed soil datasets is expected to improve the dust model performance.

In this study, we present an improvement of the dust emission scheme in the GEOS-Chem model by incorporating the updated soil texture and aerodynamic roughness length with spatial variability, the Owen effect and the drag partition correction factor, as well as the updated formulation of the sandblasting efficiency, which together significantly improve the prediction of the dust emission flux and concentrations over China. The objective is to obtain more realistic surface friction velocity (u_∗) and threshold friction velocity (u_∗t) by considering the effect of soil moisture, surface roughness and soil texture, thus improving the representation of dust emission in the model.

Section 2 gives a detailed description of the GEOS-Chem model and the modifications of the improved scheme, as well as numerical experiments and data description. Sensitivity results are compared in Sect. 3.1 to examine the impacts of the modifications. Section 3.2 presents the comparisons of the improved scheme and original version with observations, to evaluate the performance of the improved scheme. Uncertainties, limitations and future improvements of the emission scheme are discussed in Sect. 3.3, followed by a summary in Sect. 4.

2 Model and measurements

2.1 Model description

The GEOS-Chem model is a global three-dimensional chemical transport model driven by assimilated meteorology. In this work, we use a nested version of GEOS-Chem (v12.6.0) driven by the GEOS-FP-assimilated meteorological field, with a spatial resolution of 0.25^∘ × 0.3125^∘ and 72 vertical levels for the China region (15–55^∘ N, 70–140^∘ E) during the period of 27 March to 2 April 2015. The nested-grid GEOS-Chem is developed by Wang et al. (2004) with lateral boundary conditions provided by a global simulation (at 2^∘ × 2.5^∘ resolution in this study). Many nested-grid GEOS-Chem simulations have been evaluated and applied to the analysis of gaseous and aerosol species over China (Chen et al., 2009; Lin et al., 2014; Zhang et al., 2015; Wang et al., 2013; Li et al., 2019).

GEOS-Chem includes detailed atmospheric chemical mechanism and online aerosol calculations. In this work, we simulate the dust emission with a combination of the dust entrainment and deposition (DEAD) mobilization scheme (Zender et al., 2003) and Global Ozone Chemistry Aerosol Radiation and Transport (GOCART) source function. Mineral dust aerosols are distributed across four size bins (with radius bins of 0.1–1.0, 1.0–1.8, 1.8–3.0 and 3.0–6.0 µm). The mass fractions of each size bin are parameterized by the optimized dust particle size distribution proposed by Zhang et al. (2013). Dry deposition velocities for dust aerosols are computed with the gravitational settling scheme of Fairlie et al. (2007) and aerosol deposition scheme from Zhang et al. (2001). The wet deposition scheme, which includes scavenging in convective updrafts, as well as the rainout and washout of soluble tracers, is described in Liu et al. (2001). Aerosol optical depth is derived online from aerosol concentrations with an externally mixed assumption using relative humidity (RH)-dependent aerosol optical properties from Latimer and Martin (2019). Dust optics are from Ridley et al. (2012).

2.2 Improvement on the dust emission scheme in GEOS-Chem

The standard dust emission scheme in GEOS-Chem is based on a semi-empirical formulation developed by Zender et al. (2003) and is combined with the GOCART source function (Ginoux et al., 2001). In this scheme, the vertical dust flux (F) is proportional to the horizontal saltation flux (Q_s), which is the function of surface friction velocity (u_∗) and threshold friction velocity (u_∗t):

$$\begin{array}{}\text{(1)}& {\displaystyle }F=\left(\mathrm{1}-A_{\mathrm{s}}\right)S\mathit{\alpha }{Q}_{\mathrm{s}}\text{(2)}& {\displaystyle }{Q}_{\mathrm{s}}=C_z{\displaystyle \frac{{\mathit{\rho }}_{\mathrm{air}}}{g}}{u}_{\ast }^{\mathrm{3}}\left(\mathrm{1}-{\displaystyle \frac{u_{\ast t}}{u_{\ast }}}\right){\left(\mathrm{1}+{\displaystyle \frac{u_{\ast t}}{u_{\ast }}}\right)}^{\mathrm{2}}{u}_{\ast }>u_{\ast t},\end{array}$$

where α is the vertical-to-horizontal flux ratio or sandblasting efficiency, based on the soil clay content (Marticorena and Bergametti, 1995). S is based on the GOCART source function (see Supplement Fig. S1), also named the soil erodibility factor, representing the grid cell fraction of the bare land suitable for mobilization. A_s is the fraction of snow-covered surface. C_z is the saltation constant (C_z=2.61).

According to the equation, u_∗t, u_∗ and α are the key input parameters in the accurate prediction of dust emission flux. u_∗t is used to describe the characteristics of soil and land surface conditions, representing the resistance of surface to wind erosion. In the standard dust scheme, u_∗t is calculated using a semi-empirical formulation as a function of air density and soil particle density (Iversen and White, 1982). Furthermore, two correction terms, including soil moisture correction (Fécan et al., 1999) and drag partition correction (Marticorena and Bergametti, 1995), are also applied to modify u_∗t. It should be noted that in the original scheme, the drag partition correction term is eliminated.

u_∗ is a description of surface wind speed, which mainly depends on 10 m wind speed taken from meteorological field assuming neutral stability (Bonan, 1996). The Owen effect, which represents a positive feedback between saltation process and friction speed (Owen, 1964), is often adopted in models to modify u_∗. However, the Owen effect is eliminated in the original scheme.

Sandblasting efficiency α is parameterized according to the empirical relation described by Marticorena and Bergametti (1995) (MB95), which depends on the soil clay content (M_clay) and is restricted to M_clay < 20 %:

$$\begin{array}{}\text{(3)}& \mathit{\alpha }=\mathrm{100}{e}^{\left(\mathrm{134}{M}_{\mathrm{clay}}-\mathrm{6}\right)\ln \mathrm{10}}.\end{array}$$

However, in the global model, α tends to be overly sensitive to M_clay. Due to this reason, a globally fixed value of M_clay= 20 % is assumed in the current model (Zender et al., 2003).

It should be noted that some input parameters, data or formulations are quite simplified and need to be improved based on the original dust scheme described above. For example, the aerodynamic roughness length (Z₀), the smooth roughness length (Z_0s) and the mass fraction of clay in the soil (M_clay) are assumed to be uniform constants, despite the fact that they may vary with time and location. As a result, the simulation of related processes, such as the drag partition effect or soil moisture effect, may lack spatial representation. Therefore further modifications to these variables should be made in order to obtain more realistic dust emissions. Figure 1 presents the schematic diagram of the dust emission schemes in the standard model and the modifications incorporated into this study. The details of the parameterization options and required input parameters are presented in the following sections.

Figure 1 Schematic diagram and required input parameters in the dust emission scheme in GEOS-Chem.

2.2.1 Soil type and soil texture data

In the model, M_clay can have an impact on u_∗t through modifying the soil moisture correction term, thus influencing the modeled dust emission flux. The soil moisture correction term, defined as f_w, is parameterized according to Fécan et al. (1999), which accounts for the increase of u_∗t with soil moisture content.

$$\begin{array}{}\text{(4)}& {\displaystyle }{f}_w=\left\{\begin{array}{cc}\mathrm{1}& w<w^{\prime }\\ \left[\mathrm{1}+\mathrm{1.21}{\left(w-w^{\prime }\right)}^{\mathrm{0.68}}\right]& w\ge w^{\prime }\end{array}\right.\text{(5)}& {\displaystyle }{w}^{\prime }\left(\mathit{\%}\right)=a\left(\mathrm{0.0014}{M}_{\mathrm{clay}}^{\mathrm{2}}+\mathrm{0.17}{M}_{\mathrm{clay}}\right),\end{array}$$

where w is the gravimetric soil moisture and w^′ is the soil residual moisture.

With the increase of soil moisture, soil cohesion can be enhanced, particularly over regions with high clay content, thus inhibiting the sandblasting process to some extent. However, as stated above, M_clay is assumed to be a constant equal to 20 % in the original scheme, which can cause uncertainty in dust prediction. In the improved scheme, we employ the gridded data of clay content from the Global Soil Dataset for use in Earth System Models (GSDE) (Shangguan et al., 2014), which is based on the Digital Soil Map of the World (FAO, 1995, 2003) and various regional and national soil databases. Figure 2 shows the updated M_clay from Shangguan et al. (2014) with a horizontal resolution of 2^∘ × 2.5^∘ at the global scale. Compared to the original fixed value of 20 %, the updated M_clay is generally lower in most of the dust source areas over East Asia.

Figure 2 Updated input data of global M_clay (%). Data are derived from Shangguan et al. (2014) and are re-gridded to 2^∘ × 2.5^∘ horizontal resolution in the model.

2.2.2 Surface roughness length

The drag partition is used to describe the impact of roughness elements (such as rocks, pebbles and vegetation) on u_∗t. According to Marticorena and Bergametti (1995), the roughness correction term, f_d, is a function of the aerodynamic roughness length (Z₀) and the smooth roughness length (Z_0s):

$$\begin{array}{}\text{(6)}& f_d=\mathrm{1}-{\displaystyle \frac{\ln \left(\frac{z_{\mathrm{0}}}{z_{\mathrm{0}\mathrm{s}}}\right)}{\ln \left[\mathrm{0.7}{\left(\frac{\mathrm{12}\mathrm{255}\mathrm{cm}}{z_{\mathrm{0}\mathrm{s}}}\right)}^{\mathrm{0.8}}\right]}},\end{array}$$

where the required roughness lengths are set as the constant values of Z₀=0.01 and Z_0s=0.0033 cm globally.

Z₀ represents the roughness length of the overlying non-erodible elements (solid obstacles, such as rocks), which transfers part of the wind momentum from the atmosphere to the surface, dissipating the shear force for particle saltation. Prigent et al. (2005) derived global aerodynamic roughness length in arid and semi-arid areas, retrieved from the ERS-1 satellite measurements with a horizontal resolution of 0.25^∘ × 0.25^∘ . Here we apply the global monthly mean Z₀ fields provided by Prigent et al. (2005) and then re-grid the map to 2^∘ × 2.5^∘ horizontal resolution for the incorporation into GEOS-Chem. As Fig. 3 shows, compared to the fixed constant assumed in the original version, the updated global Z₀ is generally higher. Over northern China, the Z₀ value ranges from approximately 0.01 cm over desert regions to 0.07 cm.

Figure 3 Monthly updated input data of global aerodynamic roughness length (Z₀) (cm) in March (a) and April (b). Data are derived from Prigent et al. (2005) and are re-gridded to 2^∘ × 2.5^∘ horizontal resolution in the model.

Z_0s characterizes the roughness length of the uncovered, bare erodible surface. Instead of setting Z_0s to a fixed value, some studies suggested that Z_0s can be estimated as $\mathrm{1}/\mathrm{30}$ of the coarse-mode mass median diameter (MMD) of soil particles, which will provide a more realistic representation of soil texture distribution (Marticorena and Bergametti, 1995; Laurent et al., 2006; Mokhtari et al., 2012). In the improved version, we adopt this empirical formula, based on updated soil texture classification (Mokhtari et al., 2012; Xi and Sokolik, 2015), to estimate Z_0s:

$$\begin{array}{}\text{(7)}& z_{\mathrm{0}\mathrm{s}}={\displaystyle \frac{\mathrm{MMD}}{\mathrm{30}}},\end{array}$$

where MMD is the median diameter of the coarsest mode for various soil textures shown in Table 1. The corresponding Z_0s values for different soil types are listed in Table 1, and their global distribution are shown in Fig. 4. The soil texture map is obtained based on the Harmonized World Soil Database (HWSD; http://www.iiasa.ac.at/Research/LUC/External-World-soil-database/HTML, last access: 24 April 2020), which provides global sand, silt and clay contents at 30 arcsec resolution. The soil texture dataset is re-gridded to 2^∘ × 2.5^∘ resolution and then is applied to identify the global soil texture by using the United States Department of Agriculture (USDA) soil texture triangle (based on the amount of sand, clay and silt contents; http://soils.usda.gov/technical/aids/investigations/texture/, last access: 24 April 2020). There are 12 classes of soil defined by USDA. It can be seen from Fig. 5 that loam, sandy loam and clay loam are the dominant soil types over China. Among them, sandy loam and loam occupy the major part of northwest China.

Figure 4 Updated global map of smooth roughness length (Z_0s) (cm) estimated from the empirical relationship with soil texture.

Figure 5 Global soil texture map based on the USDA classification.

2.2.3 Sandblasting efficiency α

Sandblasting efficiency α is important in the dust emission calculation as it is used to convert the horizontal saltation mass flux to a vertical dust mass flux. In the original scheme, α is simply expressed as a function of M_clay, which is a fixed value of 20 %. The assumption in the original scheme might cause uncertainty in the modeled flux and make the spatial variation less representative (Mokhtari et al., 2012).

In order to reduce this uncertainty, a more physically based function from Lu and Shao (1999) (LS99) is adopted in our study. Based on wind tunnel experiments carried out by Rice et al. (1996a, b), Lu and Shao (1999) derived the expression of α through theoretical calculation and some simplifications:

$$\begin{array}{}\text{(8)}& \mathit{\alpha }={\displaystyle \frac{C_{\mathit{\alpha }}gf{\mathit{\rho }}_{\mathrm{b}}}{\mathrm{2}p}}\left(\mathrm{0.24}+C_{\mathit{\beta }}{u}_{\ast }\sqrt{{\displaystyle \frac{{\mathit{\rho }}_{\mathrm{p}}}{p}}}\right),\end{array}$$

where f is the fine particle content in the soil volume; p is the soil plastic pressure, which represents the magnitude of the surface resistance (N m⁻²); ρ_b and ρ_p are the bulk soil density and particle density, respectively; g is the gravitational acceleration (m s⁻²); u_∗ is the friction velocity (m s⁻¹); and C_α and C_β are empirical constants. Among these parameters, the values of ρ_b and p depend upon different soil textures. Some studies (Kang et al., 2011; Foroutan et al., 2017; Ma et al., 2019) have implemented this formulation in the model and proposed the proper range of these parameters over different soil types.

Table 1 Input soil aggregate size distribution parameters dependent on soil texture classification following USDA.

Including three-mode log-normal parameters (mass fraction n (%), mass median diameter MMD (µm) and geometric standard deviation σ) and smooth eolian roughness length z_0s (cm).

Many measurements from laboratory experiments and field observations have demonstrated the close relationship between α and u_∗ (Gillette et al., 1997; Gomes et al., 2003; Rajot et al., 2003; Roney and White, 2006; Macpherson et al., 2008; Panebianco et al., 2016; Zhang et al., 2016). To improve the original scheme, we extract α from these measurements over different soil types, based on the expression of LS99, as depicted in Fig. 6 and Table 2.

Figure 6 Updated sandblasting efficiency α as a function of surface friction velocity u_∗ following Lu and Shao (1999) for sand, loamy loam, loam and clay and observations from the literature.

Table 2 Input soil-related parameters for different soil texture used in calculation of the sandblasting efficiency α.

2.3 Experiment design

Several sensitivity experiments (Table 3) are conducted to assess the model performance of the modifications in the improved scheme. “Control” is the control run with the dust emission scheme originally implemented by Fairlie et al. (2007). Sen_mclay, Sen_owen, Sen_ratio, Sen_drag and Sen_Z₀Z_0s are the same as the control run but include the modification of M_clay, the Owen effect, the sandblasting efficiency, the drag partition effect and the updated surface roughness length (Z₀ and Z_0s) respectively. Sen_all represents the simulation with the improved scheme which accounts for all the modification described above.

Table 3 Sensitivity experiment design and description.

2.4 Measurements

The data used in this study include the Moderate Resolution Imaging Spectrometer (MODIS) Level 3 AOD data, hourly observational data of surface PM₁₀ concentration and meteorological field data taken from the Meteorological Information Comprehensive Analysis and Process System (MICAPS). The data used in this study are for the period of 27 March to 2 April 2015.

MODIS aerosol products are used to evaluate model results of AOD. MODIS AOD at 550 nm is obtained from the daily level-3 product from Aqua satellites (MYD08_D3, 1^∘ × 1^∘ gridded data) and is combined with Deep Blue retrievals which can provide AOD over bright surfaces (i.e., desert regions).

Hourly surface observed PM₁₀ concentration data, collected from about 1000 environmental monitoring stations maintained by the Chinese Ministry of Environmental Protection (MEP; http://www.mee.gov.cn, last access: 24 January 2020), are used to validate the model performance of surface dust concentrations.

Meteorological fields of wind speed taken from the Meteorological Information Combine Analysis and Process system (MICAPS) developed by the Chinese National Meteorological Center (NMC) are used for the evaluation of the wind field in the model. Figures S2–S3 show that the 10 m wind field used in the model scheme generally agrees well with the MICAPS observations over most sites. However, comparisons of the averaged surface wind field between the model input and observations (Fig. S4) show that although the circulation patterns in the model are identical with the observations, surface wind speed in the model tends to be larger than observations (which was also found by Wang et al., 2014), especially over western and northeastern Inner Mongolia.

3 Results and discussion

3.1 Sensitivity study

In order to assess the sensitivity of the dust emission to the modified input parameters or physical processes, several numerical experiments are conducted and compared. Figure 7a presents the relative difference (%) of averaged u_∗t during the study period between these sensitivity simulations and the control run. The u_∗t values simulated by the control run are generally small, with values less than 0.3 m s⁻¹ (not shown). Wu (2013) indicated that u_∗t values over source regions in northern China calculated by Zender et al. (2003) are generally lower (with values ranging from 0.2 to 0.25 m s⁻¹) than the measurement (with values ranging from 0.34 to 0.69 m s⁻¹) and the values calculated by Shao (2004), which is closer to the observations. The sensitivity simulations show that the update of M_clay in Sen_mclay can lead to higher u_∗t over northern China and lower u_∗t over southern China than the control simulation, which overestimates M_clay over northern China and underestimates it over southern China by setting M_clay to 20 %. In northern China, particularly in arid and semi-arid regions, the updated M_clay will decrease the soil moisture threshold w^′ and increase soil moisture term f_w, thus leading to a slight increase in u_∗t (with magnitude < 10 %). The inhibition of dust emission by surface roughness elements is not taken into account in the original scheme, i.e., f_d=1. Some studies (Darmenova et al., 2009; Menut et al., 2013) have demonstrated f_d as a function of Z₀ and Z_0s, implying that f_d values increase with Z₀ and decrease with Z_0s. Compared to the fixed values used in the original scheme, updated Z₀ fields used in Sen_Z₀Z_0s are generally larger, and updated Z_0s fields are smaller. Therefore, f_d values are increased significantly, particularly over the regions with non-erodible elements (larger Z₀). Results show that u_∗t is increased when considering the drag partition effect (increased by 10 % in Sen_drag with constant Z₀ field), particularly with the updated Z₀ and Z_0s field (increased by 10 %–60 % in Sen_Z₀Z_0s). In general, due to the inclusion of Z₀, Z_0s and M_clay, f_d and f_w are modified, which results in significant alteration in u_∗t (ranging from −8 %–72 % in Sen_all) over China. It can be found that the modification of f_d due to updated Z₀ and Z_0s makes more contribution to the increase in u_∗t.

Figure 7 Relative difference (%) in simulated averaged threshold friction velocity u_∗t (a), surface friction velocity u_∗ (b), sandblasting efficiency α (c) and emission flux (d) between sensitivity simulations and the control run during the study period.

The relative difference of u_∗ with respect to the control run is also compared in Fig. 7b. Considering the Owen effect in Sen_owen leads to an increase in u_∗ by 0 %–39 %, especially over northwest China where surface wind is strong. Modeled u_∗ is obtained from u_10 m and Z₀ under neutral conditions (Bonan, 1996). It can be seen that updated Z₀ in Sen_Z₀Z_0s can modify u_∗ by influencing the boundary-layer exchange properties. u_∗ over northern China is generally increased by 10 %–22 %, with higher values of Z₀ in Sen_Z₀Z_0s, while it is slightly decreased over the Taklimakan and Gobi deserts. In Sen_all, modeled u_∗ is increased by 5 %–50 % over most parts of China.

Figure 7c presents the percentage difference in terms of sandblasting efficiency α. In the original version, α is set as a uniformly constant value (around 0.04) due to the assumption of a fixed M_clay. In Sen_ratio and Sen_all, the u_∗-dependent ratio, following LS99, which varies with different soil texture according to observations, is adopted. On average, α is decreased by 50 % with the modification in Sen_ratio and Sen_all. The largest reduction occurs over regions with sand texture, such as over the Taklimakan and Gobi deserts.

As seen from Fig. 7d, the simulated dust emission flux (F) varies significantly among different experiments. Due to the inclusion of updated M_clay, the soil moisture term increases in Sen_mclay, which leads to higher u_∗t and lower F over most regions. Accounting for the Owen effect in Sen_owen results in an increase in F of 0 %–314 %, particularly over the northern part of Gansu Province and northwestern Inner Mongolia. A significant reduction in arid and semi-arid regions of northern China is caused by updated α (Sen_ratio). In Sen_drag and Sen_Z₀Z_0s, the alteration of F values ranges from −100 % to −4 % and −100 % to 50 %, respectively, due to the inclusion of f_d with constant Z₀ and updated Z₀ and Z_0s. Due to the combined effects of the modifications, F simulated by Sen_all is generally reduced over northern China, except in some regions of northwest China, where the Owen effect plays a dominant role.

Five sites closer to the dust source area or significantly influenced by dust storms (Xilinguole, Huhehaote, Jiuquan, Kuele and Akesu; locations shown in Fig. S1) are selected to evaluate the performance of the control run and sensitivity simulations. Comparisons of the modeled u_∗t (Fig. 8) show that in all sites, modeled u_∗t values are increased in Sen_mclay, Sen_drag, Sen_Z₀Z_0s and Sen_all, compared with the original model, with the highest u_∗t simulated by Sen_all. Modeled u_∗t values increase from 0.23–0.25 m s⁻¹ in the control run to 0.29–0.37 m s⁻¹ in Sen_all. The reported u_∗t values over arid and semi-arid regions of China are around 0.3–0.5 m s⁻¹ (Wang et al., 2009). Wu (2013) summarized that u_∗t range from 0.34–0.69 m s⁻¹ over East Asia and indicated that u_∗t values calculated by Zender et al. (2003) are relatively lower, ranging from 0.2 to 0.25 m s⁻¹. It is apparent that modeled u_∗t values are greatly increased in the revised simulation, which is much closer to the observed values. This improvement is mainly attributed to the update of Z₀ and Z_0s. Comparisons between the modeled averaged PM₁₀ concentrations and the observational values in five sites show that PM₁₀ levels simulated by Sen_all are closer to observations than many other cases. In summary, Sen_all shows the better agreement with the observations in terms of u_∗t and PM₁₀ concentrations.

Figure 8 Comparisons of simulated averaged threshold friction velocity u_∗t (a) and PM₁₀ concentrations (b) at selected sites. Black dots in (b) indicate the observed averaged PM₁₀ concentrations.

3.2 Comparison between the improved scheme and the original scheme with observations

In order to validate the model performance of the improved scheme, time series of the observed surface PM₁₀ concentrations are compared against the modeled values from the control run (the original scheme) and Sen_all (the improved version) during a dust episode from 27 March to 2 April 2015 over northern China. The intensity and evolution of this dust event have been described by Wang et al. (2017), illustrating that dust particles were mainly emitted from Mongolia and the Inner Mongolia region of China and that a dust backflow event took place over northern China on 29 March. Figure 9 compares the hourly modeled PM₁₀ concentrations and observed values for nine selected sites (locations shown in Fig. S3), which are closer to the dust sources or severely affected by the dust event. It shows that the dust concentrations are generally underestimated in the control run, particularly when dust concentrations are quite high, indicating that the original scheme has difficulty in accurately reproducing the dust emission process. Sen_all generally reproduce the PM₁₀ levels better than the control run. Both the magnitude and the temporal evolution of PM₁₀ concentrations are captured in Sen_all quite well, with peak values much closer to the observations. Among these sites, Sen_all shows better performance over North China, e.g., Beijing, Tianjin and Huhehaote. But both the control run and Sen_all fail to capture the peak values from 29 to 30 March. During this period, dust particles, mixed with anthropogenic pollutants, flew back due to the south wind over North China (Wang et al., 2017). Uncertainties in the meteorological field and dust heterogeneous reactions in the model may cause the model bias.

Figure 9 Temporal variation of hourly PM₁₀ concentrations from observations (black dots) and simulations of the control run (blue line) and Sen_all (red line) during the study period at nine selected sites.

For specific periods, however, modeled peak values of some sites occur earlier (several hours) than the observations at some sites (e.g., Beijing and Tianjin in 28 March), which could be considered as being a result from the uncertainty in the wind field used in the model. It shows that the surface wind is stronger in the model than the observations (Fig. S2), which may lead to stronger transport of the dust from source regions to downwind areas such as Beijing, Tianjin and Kuele. Instead, modeled and observed peak values of some sites in the source regions (e.g., Huhehaote, Xilinguole and Hami) almost occur simultaneously.

In order to quantify the performance of the model result, some statistical parameters, including the mean values, correlation coefficient (R), mean bias (MB) and normalized mean bias (NMB), are calculated and listed in this paper. The statistical performance for the modeled surface PM₁₀ concentrations from the control run and Sen_all against observations is presented in Table 4. It shows that dust concentrations at all selected sites are significantly underestimated in the control run, especially over northwest regions, with the MB and NMB values ranging from −163.5 to −503.61 µg m⁻³ and −64.61 % to −68.48 % respectively. It is obvious that Sen_all with updated modification greatly improves the dust concentration prediction, with mean values more comparable to the observations, and the average MB and NMB values reduce from −196.29 µg m⁻³ and −52.79 % in the control run to −47.72 µg m⁻³ and −22.46 % respectively. The largest improvement occurs at northwest stations (e.g. Hami, Akesu and Kuele), which are located close to the Taklimakan desert. Over other regions, such as North China (e.g., Beijing, Tianjin, Huhehaote and Xilinguole), the model performance of Sen_all is slightly better than the control run.

Table 4 Statistics for observed and simulated (control and Sen_all) surface PM₁₀ concentrations at selected sites.

Although the MB and NMB values of most stations are generally lower and the mean values are much closer to observations for the Sen_all simulation, i.e., modifications included improve the underestimation in the control run to some extent, the dust concentrations are still generally underestimated. For stations closer to Gobi desert, such as Xilinguole, Jiuquan and Baiyin, dust concentrations are greatly underestimated, with NMB < −30 %, which is likely due to the uncertainty in the erodibility factor over Gobi desert used in our study (Ginoux et al., 2001). Similarly, Su and Fung (2015) evaluated the performance of dust emission schemes in WRF-Chem over East Asia, pointing out that the erodibility factor from Ginoux et al. (2001) over the Gobi desert is significantly underestimated, which may result in the underestimation of the dust emission over the Gobi desert. Given that the simulated dust emission flux is directly scaled by the erodibility factor, we suggest that the erodibility factor used in our model needs to be updated or improved.

As stated above, although the model can capture the overall temporal variations of surface dust concentrations, modeled peak values occur earlier (about 6 h) than the observations over several stations, which may be attributed to strong transport due to stronger surface wind used in the model. It should be noted here that this model bias contributes a lot to the simulation error, leading to smaller R and greater MB and NMB values. R values will be greatly improved if this bias is eliminated, implying that the input assimilated meteorological field is important for dust emission simulation and needs to be further evaluated and adjusted.

The averaged modeled surface PM₁₀ concentrations with and without modifications (Sen_all and control run respectively) and observational values at ∼ 1400 stations over China during the study period are compared in Fig. 10. It shows that dust concentrations are generally underestimated in the control run (NMB = −16 %, regression slope = 0.4), which could be attributed to the crude representation of soil properties, roughness length and other related elements. Incorporating improvements into the scheme makes the modeling result much closer to the observations, with R values increasing from 0.6 to 0.7, NMB values changing from −16 % to −11 % and the regression slope ranging from 0.4 to 0.6. However, the improved model still tends to underestimate dust concentrations. Unrealistic soil properties (e.g., soil texture or roughness length) and an insufficiently accurate potential source map (the erodibility factor) used to scale the dust emission flux could be the possible causes.

Figure 10 Comparison of modeled and measured surface PM₁₀ concentrations at observational sites. The dotted line is the 1:1 line. Model results are taken from the control run and Sen_all respectively.

To further investigate the model performance, spatial distributions of averaged simulated surface PM₁₀ concentrations from the control run and Sen_all and their comparisons against observations are presented in Fig. 11. Results show that both the control run and Sen_all can reproduce the pattern of dust concentrations in the study region, with high values located over northwest China, North China and some areas of northeast China, indicating that GEOS-Chem can represent the main features of dust emission and transport during the dust storm. It is found that for most sites in the control run, the simulated magnitudes are close to the observational values but are underestimated over northwest China (where Gobi and Taklimakan deserts are located) and the North China Plain. The simulated values from Sen_all are generally larger than the control run and are more consistent with measurements both in magnitude and in area extent, especially over the desert region of northwest China. However, dust concentrations are still underestimated over the North China Plain, possibly due to the outdated source map or some potential dust source regions over Inner Mongolia not being included well. In addition, the missing mechanism of secondary aerosol sources in the model such as heterogeneous reactions could also cause the model bias (Zheng et al., 2015; Cheng et al., 2016).

Figure 11 Comparison of simulated averaged PM₁₀ surface concentrations from the control run (a) and Sen_all (b) with the observed values.

Figure 12 shows the spatial distribution of simulated averaged AOD from the control run and Sen_all as well as MODIS AOD for the study period. For better comparison, simulated AOD at 13:00 local time (Aqua passage time) is extracted. Results show that the control run reproduces the major regions with high AOD values, e.g., eastern China, but with lower magnitude. The control run also fails to capture the high-AOD area over the Taklimakan desert, while Sen_all could capture it. Compared with the control run, Sen_all generally reproduces the spatial coverage and magnitude of the observed AOD.

Figure 12 Spatial distribution of MODIS-retrieved AOD at 550 nm (a) and simulated AOD at 550 nm from the control run (b) and from Sen_all (c). The simulation results are extracted at 13:00 local time to match the MODIS observation time.

3.3 Discussion

In our study, we point out that the erodibility factor (S) in the model may introduce uncertainty in modeled dust concentrations, especially over Gobi desert. Several studies indicated that S from Ginoux et al. (2001) over the Gobi Desert has been significantly underestimated and needs to be improved (Su and Fung, 2015; Zeng et al., 2020). Wu and Lin (2014) have demonstrated that the potential source regions in the southeast of Mongolia and the middle-east of Inner Mongolia are not well characterized by the S from the GOCART scheme, which results in the underestimation of dust concentration in this area and its downwind regions. In addition, the source function may not provide precise enough information about the recent expansion of dust source areas over northern China, with desertification and deforestation (Ku and Park, 2013). Studies have demonstrated that implementing a physically based parameterization instead of an empirical dust source function, which is usually time-invariant and lacks physical treatment (Kok et al., 2014a, b), or adopting the dynamic dust source function (Xi and Sokolik, 2015) could improve the representation of dust emission. Therefore, the dust source function should be precisely established with new updates and higher resolution using various measurements.

In terms of sandblasting efficiency α, many modeling studies as well as observational analysis have investigated its magnitude and expression, but the results may vary greatly (Kang et al., 2011; Ma et al., 2019). The formulation for α used in our improved scheme is based on LS99, which establishes the relationship between α and u_∗, along with other soil-related parameters dependent on soil textures. In this study, we derived α for different soil types based on the reported values, but uncertainties still remain due to limited available measurements. In addition to the expression from MB95 and LS99, there are other α formulations proposed by Shao et al. (1996) (Shao96) and Shao (2004) (Shao04). Different from the empirical function, expressions of Shao96 and Shao04 are more sophisticated, which are the function of u_∗t and u_∗ respectively, along with some size information of soil particles. Comparisons of different formulations of α for different soil types (Kang et al., 2011; Ma et al., 2019) have shown that the variation in α can reach up to several orders of magnitude, and few equations could reproduce the measured positive correlation between α and u_∗, suggesting that α for different soil texture should be further investigated and observed to improve the model accuracy.

In this study, surface conditions including erodibility factor, soil texture, clay content and surface roughness length play a significant role in improving the model performance of u_∗, u_∗t and F. We conclude that substituting globally fixed values of these properties with more realistic and physical-based ones could reduce the model uncertainty and improve the understanding of the dust emission mechanism. In the physically based scheme, the importance of accurate input surface properties, including soil particle size distribution (Darmenova et al., 2009; Kok, 2011a, b), soil texture (Shao et al., 2011; Foroutan et al., 2017), surface roughness length (Darmenova et al., 2009; Kontos et al., 2018) and soil plastic pressure (Lu and Shao, 1999; Kang et al., 2011), has also been highlighted by many studies. Therefore, accurate and abundant observation data of soil-related properties are urgently needed, particularly over dust source region. Moreover, various and comprehensive observation methods (e.g., experimental data and field and satellite observations) are recommended in order to correct and update the input data.

4 Summary and conclusion

In this study, we revised the treatments of dust emission processes by considering the effect of soil moisture, surface roughness and soil texture, as well as the Owen effect and a more physically based formulation of the sandblasting efficiency in GEOS-Chem version 12.6.0, in order to improve dust simulation over China. Several sensitivity simulations were conducted during a severe dust storm between 27 March to 2 April 2015 over northern China to analyze the effects of these modifications on u_∗, u_∗t and the emission flux.

In the improved scheme, we substituted assumed constant values of Z₀ and M_clay in the original version with geographical variation maps obtained from the measurement provided by Prigent et al. (2005) and Shangguan et al. (2014) respectively. Z_0s and sandblasting efficiency were calculated with formulations based on soil texture data from the FAO dataset, which is more physically based than the original version. In addition, the Owen effect and drag partition correction factor were considered in the improved version.

Sensitivity results showed that the modified f_d and f_w by inclusion of the updated Z₀, Z_0s and M_clay resulted in significant alteration in u_∗t (ranging from −8 %–72 %) over China. u_∗t was increased when including the drag partition effect, particularly with the updated Z₀ and Z_0s field (increased by 10 %–60 %), which brought the modeled u_∗t much closer to the measurements. Considering the Owen effect increased modeled u_∗ by 0 %–39 %, especially over northwest China where surface wind is strong. In general, modeled u_∗ was increased by 5 %–50 % over most parts of China due to the inclusion of the Owen effect and updated Z₀. In terms of sandblasting efficiency, it was decreased by 50 % on average with the updated u_∗-dependent ratio following LS99, with the largest reduction occurring over regions with sand texture. Due to the combined effect of updated treatments, the emission flux simulated by the improved scheme was generally decreased over northern China, except in some regions of northwest China, where the Owen effect played a dominant role. Better agreement between the improved model results and observational values was achieved in terms of the u_∗t and surface PM₁₀ concentrations in selected typical sites over northern China.

Compared with both surface PM₁₀ observations and MODIS AOD, the revised dust emission scheme produced a better performance in both temporal and spatial variation. Result showed that the dust concentrations were generally underestimated at selected sites in the original scheme, particularly when dust concentrations were high. For the improved scheme, both the magnitude and the temporal evolution of PM₁₀ concentrations were well captured, with peak values much closer to the observations. According to the statistics, with the implementation of the updates, averaged PM₁₀ values at selected sites were more comparable to the observations, and the average MB and NMB values were reduced from −196.29 µg m⁻³ and −52.79 % in the control run to −47.72 µg m⁻³ and −22.46 % respectively. However, for some sites closer to Gobi desert, dust concentrations were still underestimated, which was likely attributed to the uncertainty in the erodibility factor over Gobi desert. Comparison of the model results and observed averaged PM₁₀ concentrations at ∼ 1000 stations showed that the revised scheme improved the model performance, with R values increasing from 0.6 to 0.7 and NMB values changing from −16 % to −11 %. Moreover, the improved scheme demonstrated better performance in reproducing the spatial distribution of AOD than the original scheme, particularly over the desert region of northwest China.

In summary, this study indicated that compared to the original scheme, the revised dust emission scheme had an overall better agreement with the measurements. However, more physically based schemes and more detailed up-to-date input parameters should be further investigated and observed to improve the accuracy of model.