WRF-Chem v3.9 simulations of the East Asian dust storm in May 2017: modeling sensitivities to dust emission and dry deposition schemes

Dust aerosol plays an important role in the radiative budget and hydrological cycle, but large uncertainties remain for simulating dust emission and dry deposition processes in models. In this study, we investigated dust simulation sensitivity to two dust emission schemes and three dry deposition schemes for a severe dust storm during May 2017 over East Asia using the Weather Research and Forecasting model coupled with chemistry (WRF-Chem). Results showed that simulated dust loading is very sensitive to different dry deposition schemes, with the relative difference in dust loading using different dry deposition schemes ranging from 20 %–116 %. Two dust emission schemes are found to produce significantly different spatial distributions of dust loading. The difference in dry deposition velocity in different dry deposition schemes comes from the parameterization of collection efficiency from impaction and rebound effect. An optimal combination of dry deposition scheme and dust emission scheme has been identified to best simulate the dust storm in comparison with observation. The optimal dry deposition scheme accounts for the rebound effect and its collection efficiency from impaction changes with the land use categories and therefore has a better physical treatment of dry deposition velocity. Our results highlight the importance of dry deposition schemes for dust simulation.

1 Introduction

Dust aerosol is an important component in the atmosphere, and it can impact many processes of the Earth system. Through absorbing and scattering shortwave and longwave radiative fluxes, dust can alter the radiative budgets, which is called the direct effect (Chen et al., 2013; Kok et al., 2017; Zhao et al., 2010, 2011, 2012). Acting as cloud condensation nuclei (CCN) and ice nuclei (IN), dust can change cloud properties and precipitation, which are called the indirect effects (Creamean et al., 2013; Demott et al., 2010). Besides, dust aerosol can absorb solar radiation and change the atmospheric stability and therefore cloud formation, which is known as the semidirect effect (Hansen et al., 1997). Furthermore, natural dust is important for air quality assessments and has significant impacts on human health (Abuduwaili et al., 2010; Chen et al., 2019; Hofer et al., 2017; Jiménez-Guerrero et al., 2008; Ozer et al., 2007). Although great progress has been made in dust models and dust simulations in recent decades, large uncertainties remain in dust simulations (Huneeus et al., 2011; Prospero et al., 2010; Todd et al., 2008; Uno et al., 2006; Zender et al., 2004; Zhao et al., 2013).

A complete description of dust events includes dust emission, deposition, and transport processes. The differences in dust simulation mainly result from the uncertainties in dust emission, deposition, and transport processes in models. One uncertainty is from dry deposition processes. Dry deposition refers to the transport of particles from the atmosphere to the Earth's surface in the absence of precipitation (Seinfeld and Pandis, 2006). In most aerosol modeling, dry deposition velocity V_d is used to calculate the dry deposition flux, and V_d is usually modeled using the resistance-based approach (Pryor et al., 2008). In the resistance-based approach, V_d is determined by gravitational settling, aerodynamic resistance, and surface resistance. Surface resistance is determined by collection efficiency from Brownian diffusion, impaction, and interception and is corrected for particle rebound. Slinn (1982) proposed a semi-analytical description of particle collection efficiencies based on the wind tunnel studies, and many dry deposition schemes since then are variants of this model (Binkowski and Shankar, 1995; Giorgi, 1986; Peters and Eiden, 1992; Zhang et al., 2001). As the formulations for collection efficiencies from different dry deposition schemes are derived from measurements that have been obtained under different meteorological conditions and land surface types, there remains a large discrepancy of these formulations between different dry deposition models (Petroff et al., 2008).

At present, the comparisons of different dry deposition schemes with reliable field measurements are mainly focused on 1D dry deposition models (Hicks et al., 2016; Khan and Perlinger, 2017; Petroff et al., 2008; Ruijrok et al., 1995). For example, Hicks et al. (2016) compared five deposition models with observations and found that V_d predicted for particles less than 0.2 µm is consistent with the measurements, but predicted V_d can vary greatly in the size range of 0.3 to about 5 µm. However, few studies have been conducted to study how different dry deposition schemes affect aerosol concentrations and their spatial distribution in the 3D numerical models. Wu et al. (2018) compared the effects of different dry deposition schemes on black carbon simulation in a global climate model (CESM-CAM5) but did not examine how different dry deposition schemes affect aerosol concentrations for large-size aerosol particles (e.g., diameters>2.5 µm), such as dust.

Another uncertainty in dust simulation is the treatment of the dust emission process in models. Natural dust is typically emitted from dry, erodible surfaces when the wind speed is high. Dust emission process is closely related to soil texture, soil moisture content, surface conditions, atmospheric stability, and the wind velocity (Marticorena and Bergametti, 1995). Dust emission schemes are used to predict the dust emission flux and to describe the dust size distribution. Many studies have compared and evaluated the performance of different dust emission schemes (Kang et al., 2011; LeGrand et al., 2019; Su and Fung, 2015; Wu and Lin, 2013, 2014; Yuan et al., 2019; Zhao et al., 2010, 2006). These studies show a large diversity of simulated dust emission flux among different dust emission schemes. Zhao et al. (2006) implemented two dust emission schemes in the NARCM (Northern Aerosol Regional Climate Model) regional model and found that both schemes captured the dust mobilization episodes and produced the similar spatial distributions of dust loading over East Asia, but significant differences exist in the dust emission fluxes and surface concentrations. Kang et al. (2011) compared three dust emission schemes in WRF-Chem and found that the difference between the vertical dust fluxes derived from the three emission schemes can reach to several orders of magnitude. Yuan et al. (2019) found that one scheme strongly underestimated the dust emission, while another two schemes can better show the spatial and temporal variation in dust aerosol optical depth (AOD) based on WRF-Chem simulation of a storm outbreak in Central Asia. In another WRF-Chem study, Chen et al. (2017a) concluded that the dust emission differences mainly come from the dust emission flux parameterizations and differences in soil and surface input parameters in different dust emission schemes.

While dust emission schemes have been studied quite extensively, few studies have examined dust emission and dry deposition schemes simultaneously. As both dust emission schemes and dry deposition schemes contribute significantly to the uncertainties in dust simulations, evaluating dust emission schemes based on a single dry deposition scheme may be problematic, especially if the dry deposition schemes employed have deficiencies. For example, as a widely used regional model that has been coupled with a variety of dust emission schemes, the WRF-Chem model has been used in many studies to evaluate the performance of the dust emission schemes (LeGrand et al., 2019; Su and Fung, 2015; Wu and Lin, 2013, 2014; Yuan et al., 2019). But most of these studies use the GOCART aerosol scheme, and only one dry deposition scheme (Wesely et al., 1985) is coupled within the GOCART aerosol scheme. Zhang et al. (2019) compared the modeled dust deposition using the GOCART aerosol scheme in WRF-Chem with observed dust deposition and found that modeled dust deposition is highly underestimated by more than 1 order of magnitude compared to the observed deposition. This indicates that the dry deposition scheme (Wesely et al., 1985) in the GOCART aerosol scheme may not be suitable for dust simulation and needs to be further improved.

In this study, we adopted the MOSAIC aerosol scheme coupled within the WRF-Chem model to study how dry deposition schemes and dust emission schemes affect dust simulations by evaluating model results against observations. As the MOSAIC aerosol scheme includes several different dry deposition schemes, this allows us to choose more advanced dry deposition schemes. As the default MOSAIC aerosol scheme only includes the GOCART dust emission scheme, we further implemented the dust emission scheme Shao2011 (Shao et al., 2011) in the MOSAIC aerosol scheme, which allows us to compare these two widely used dust schemes along with multiple dry deposition schemes. The goals of this study are (1) to study dust simulation sensitivity to different dust emission schemes and dry deposition schemes and (2) to explore which combination of dust emission scheme and dry deposition scheme can better simulate dust storms in East Asia. The paper is organized as follows. Section 2 introduces the WRF-Chem model, dust emission schemes and dry deposition schemes used, experiments design and measurements. Section 3 analyzes the dust simulation sensitivity to dust emission schemes and dry deposition schemes and the comparisons with observations. Section 4 is the summary and discussion.

2 Methodology and measurements

2.1 Model description

In this study, WRF-Chem version 3.9 is used. WRF-Chem is built based on the regional mesoscale model WRF and fully coupled with gas and aerosol chemistry module (Grell et al., 2005). A summary of the settings used to configure the model is listed in Table 1. The Noah land surface model (Chen and Dudhia, 2001) and the Yonsei University (YSU) planetary boundary scheme (Hong et al., 2006) are used in this study. The global soil categorization data set from the United States Geological Survey (USGS) with 24 land categories are used. The Rapid Radiative Transfer Model for General Circulation (RRTMG) radiation scheme (Iacono et al., 2008) is used to calculate the longwave and shortwave radiation. The Grell–Freitas convective scheme (Grell and Freitas, 2014) and the Morrison two-moment microphysics scheme (Morrison et al., 2008) are used. The gas-phase chemistry module used is the Carbon Bond Mechanism version Z (CBMZ; Zaveri and Peters, 1999). The aerosol module used here is the Model for Simulating Aerosol Interactions and Chemistry with four bins (MOSAIC 4-bin) (Zaveri et al., 2008). The MOSAIC 4-bin aerosol scheme divides airborne particles into four size bins by their effective diameter (0.039–0.156, 0.156–0.625, 0.625–2.5, and 2.5–10.0 µm) to represent aerosol size distribution. The first three bins represent the Aitken mode and accumulation mode of aerosol. The last bin represents the coarse mode of aerosol. The MOSAIC aerosol scheme includes sulfate, methane sulfonate, nitrate, chloride, carbonate, ammonium, sodium, calcium, black carbon (BC), primary organic mass (OC), liquid water, and other inorganic mass (OIN). The OIN species include silica, other inert minerals, and trace metals. The emitted dust is assigned to the OIN class of MOSAIC to simulate the major aerosol processes. To study the sensitivity of dust simulation to different dust emission schemes and dry deposition schemes, we test two different dust emission schemes (see Sect. 2.2) and three dry deposition schemes (see Sect. 2.3) within MOSAIC.

Table 1 WRF-Chem configuration.

2.2 Dust emission schemes

Dust emission schemes include empirical schemes and schemes based on dust physical processes. Because of differences in input parameters and formulas to calculate dust flux, dust emission varies among different dust emission schemes. The Goddard Chemistry Aerosol Radiation and Transport (GOCART) dust emission scheme (Ginoux et al., 2001) is an empirical scheme and was implemented in MOSAIC by Zhao et al. (2010). The GOCART dust emission scheme within the MOSAIC aerosol scheme is called by setting dust_opt=13. The University of Cologne (UoC) dust emission schemes (Shao, 2001, 2004; Shao et al., 2011) (Shao schemes) are size-resolved dust emission scheme based on the wind erosion physical theory. The UoC dust emission scheme within the GOCART aerosol scheme is called by setting dust_opt=4. When the UoC dust emission scheme is selected, the user should also choose one of the UoC sub-options by setting dust_scheme=1 for Shao2001, dust_schme=2 for Shao2004, or dust_schme=3 for Shao2011. The Shao dust emission schemes are widely used for dust simulations in East Asia and have been found to perform well in simulating dust emission fluxes (Shao et al., 2011; Su and Fung, 2015; Wu and Lin, 2014), and Shao2011 (Shao et al., 2011) is a simplified version of Shao2004 (Shao, 2004). To test the sensitivity of dust simulation to different dust emission schemes, we implemented the Shao2011 dust emission scheme in the MOSAIC aerosol scheme. Each dust emission scheme is described in detail below.

2.2.1 GOCART

The formula of vertical dust flux in GOCART is approximated as

$$\begin{array}{}\text{(1)}& F_{\mathrm{p}}=\left\{\begin{array}{l}CS{s}_{\mathrm{p}}{u}_{\mathrm{10}}^{\mathrm{2}}\left(u_{\mathrm{10}}-u_{\mathrm{t}}\right)\mathrm{if}{u}_{\mathrm{10}}>u_{\mathrm{t}}\\ \mathrm{0}\mathrm{otherwise},\end{array}\right.\end{array}$$

where C is an empirical proportionality constant, and S is the source function that is determined by the erodibility factor (see Fig. 1). s_p is the fraction of each size class of the emitted dust. u₁₀ is the horizontal wind speed at 10 m. u_t is the threshold velocity below which the dust emission does not occur. u_t is calculated as

$$\begin{array}{}\text{(2)}& u_{\mathrm{t}}=u_{\mathrm{t}\mathrm{0}}\times (\mathrm{1}+\mathrm{1.2}{\log }_{\mathrm{10}}w),\end{array}$$

where u_t0 is the threshold velocity for dry soil, and w is the soil surface wetness. The formula of u_t0 is not from the original GOCART paper (Ginoux et al., 2001) but rather from Marticorena and Bergametti (1995).

$$\begin{array}{}\text{(3)}& u_{\mathrm{t}\mathrm{0}}=\mathrm{0.129}{\displaystyle \frac{{\left(\frac{{\mathit{\rho }}_{\mathrm{p}}g{d}_{\mathrm{p}}}{{\mathit{\rho }}_{\mathrm{a}}}\right)}^{\mathrm{0.5}}{\left(\mathrm{1}+\frac{\mathrm{0.006}}{{\mathit{\rho }}_{\mathrm{p}}g{d}_{\mathrm{p}}^{\mathrm{2.5}}}\right)}^{\mathrm{0.5}}}{{\left[\mathrm{1.928}{\left(a{\left(d_{\mathrm{p}}\right)}^x+b\right)}^{\mathrm{0.092}}-\mathrm{1}\right]}^{\mathrm{0.5}}}},\end{array}$$

where ρ_p is the density of particles, ρ_a is the density of air, d_p is particle diameter, a equals 1331, x equals 1.56, and b equals 0.38. The original GOCART dust emission scheme in the GOCART aerosol scheme (dust_opt=1) calculates the dust emission flux from 0.2 to 20 µm. For the GOCART dust scheme in the MOSAIC aerosol scheme (dust_opt=13), the total dust emissions from 0.2 to 20 µm are redistributed to the size bins of MOSAIC (0.039–0.156, 0.156–0.625, 0.625–2.5, and 2.5–10.0 µm) with mass fractions of 0 %, 0.38 %, 8.8 %, and 68.0 % (Kok, 2011; Zhao et al., 2013). We note that in addition to the size distribution, the values of empirical proportionality constant C are also different for the two GOCART dust emission scheme options. For dust_opt=13, C value is set to $\mathrm{1.0}\times {\mathrm{10}}^{-\mathrm{9}}$ kg s² m⁻⁵, which is consistent with the original GOCART dust emission scheme paper (Ginoux et al., 2001). For dust_opt=1, C value is set to $\mathrm{0.8}\times {\mathrm{10}}^{-\mathrm{9}}$ kg s² m⁻⁵.

Figure 1 Domain map for the WRF-Chem simulations. The color shading shows the erodibility factor, which is the fraction of erodible surface in each grid cell.

2.2.2 Shao2011

The Shao2011 dust emission scheme is a size-resolved dust emission scheme based on the wind erosion physical theory. The dust flux is determined by

$$\begin{array}{}\text{(4)}& F\left(d_i\right)=c_y{\mathit{\eta }}_{\mathrm{mi}}\left(\mathrm{1}+{\mathit{\sigma }}_{\mathrm{m}}\right){\displaystyle \frac{gQ}{u_{\ast }^{\mathrm{2}}}},\end{array}$$

where F(d_i) is the dust emission rate of particle size d_i; c_y is the dimensionless coefficient; η_mi is the mass fraction of free dust for a unit soil mass; σ_m is bombardment efficiency; Q is the saltation flux averaged over the range of sand particle sizes. In Shao2011, the erodibility factor is only used to constrain the potential emission regions. Dust emission is permitted in Shao2011 where the erodibility factor is greater than zero. As the Shao2011 scheme is a size-resolved dust emission scheme, it first calculates the emitted dust from 0.98 um to 20 µm with 40 size bins. Dust emissions from these 40 size bins are then grouped into the four size bins of the MOSAIC aerosol scheme (0.039–0.156, 0.156–0.625, 0.625–2.5, and 2.5–10 µm). The details of the Shao2011 dust emission scheme are described in Appendix A. There is a bug in calculating dust emission flux in Shao2011 scheme reported after WRF-Chem v3.9, and we have already corrected it in our simulation (see Appendix A). We should mention that the Shao2011 dust emission scheme used in this study is based on WRF-Chem v3.9 with some modifications from WRF-Chem v3.7.1. The differences of Shao2011 among different WRF-Chem versions are documented in Appendix B.

2.3 Dry deposition schemes

For dry deposition schemes, dry deposition velocity (V_d) is used to calculate dry deposition flux. V_d is determined by gravitational settling velocity (V_g), aerodynamic resistance (R_a), and surface resistance (R_s). There are three dry deposition schemes available in WRF-Chem coupled with the MOSAIC module and used in this study, and they are referred to as BS95 (Binkowski and Shankar, 1995), PE92 (Peters and Eiden, 1992), and Z01 (Zhang et al., 2001). Each dry deposition scheme will be described in detail below.

2.3.1 BS95

In the BS95 scheme (Binkowski and Shankar, 1995), V_d is expressed as

$$\begin{array}{}\text{(5)}& V_{\mathrm{d}}=V_{\mathrm{g}}+{\displaystyle \frac{\mathrm{1}}{R_{\mathrm{a}}+R_{\mathrm{s}}+R_{\mathrm{a}}{R}_{\mathrm{s}}{V}_{\mathrm{g}}}},\end{array}$$

where R_a and R_s are aerodynamic and surface resistance; V_g is the gravitational settling velocity and is given as

$$\begin{array}{}\text{(6)}& V_{\mathrm{g}}={\displaystyle \frac{{\mathit{\rho }}_{\mathrm{p}}{d}_{\mathrm{p}}^{\mathrm{2}}g{C}_{\mathrm{c}}}{\mathrm{18}\mathrm{\mu }}},\end{array}$$

where C_c is the Cunningham correction factor as a function of d_p and mean free path of air (λ), and μ is the viscosity dynamic of air. The surface resistance is calculated as

$$\begin{array}{}\text{(7)}& R_{\mathrm{s}}={\displaystyle \frac{\mathrm{1}}{u_{\ast }\left(E_{\mathrm{B}}+E_{\mathrm{IM}}\right)}},\end{array}$$

where E_B is collection efficiency from Brownian diffusion. E_B is calculated as follows:

$$\begin{array}{}\text{(8)}& E_{\mathrm{B}}=S{c}^{-\frac{\mathrm{2}}{\mathrm{3}}},\end{array}$$

where Sc is the Schmidt number, given by $Sc=\mathit{\nu }/D$. ν is the kinematic viscosity of air, and D is the particle Brownian diffusivity. E_IM is the collection efficiency due to impaction of the particle with the collecting surface (Gallagher, 2002). Impaction occurs when there are changes in the direction of airflow, and particles that cannot follow the flow will collide with the obstacle and stay on the surface due to the inertia (Giardina and Buffa, 2018). E_IM is given by

$$\begin{array}{}\text{(9)}& E_{\mathrm{IM}}={\mathrm{10}}^{-\frac{\mathrm{3}}{St}},\end{array}$$

where St is the Stokes number, given by

$$\begin{array}{}\text{(10)}& St={\displaystyle \frac{u_{\ast }^{\mathrm{2}}{V}_{\mathrm{g}}}{g\mathit{\nu }}}.\end{array}$$

St is the ratio of the particle stop distance to the characteristic length of the flow and describes the ability of particles to adopt the fluid velocity (Pryor et al., 2008; Seinfeld and Pandis, 2006).

2.3.2 PE92

In PE92 scheme (Peters and Eiden, 1992), the dry deposition velocity (V_d) is expressed as

$$\begin{array}{}\text{(11)}& V_{\mathrm{d}}=V_{\mathrm{g}}+{\displaystyle \frac{\mathrm{1}}{R_{\mathrm{a}}+R_{\mathrm{s}}}}.\end{array}$$

The formula of V_g and R_a is the same as in BS95, but the way to calculate R_s is quite different. In PE 92, R_s is parameterized as

$$\begin{array}{}\text{(12)}& R_{\mathrm{s}}={\displaystyle \frac{\mathrm{1}}{u_{\ast }\left(E_{\mathrm{B}}+E_{\mathrm{IM}}+E_{\mathrm{IN}}\right)R}},\end{array}$$

where E_IN is collection efficiency from interception, and R is the factor for particle rebound. E_IM, E_IN, and R are expressed as

$$\begin{array}{}\text{(13)}& E_{\mathrm{IM}}={\left({\displaystyle \frac{St}{\mathrm{0.8}+St}}\right)}^{\mathrm{2}},\end{array}$$

$$\begin{array}{}\text{(14)}& E_{\mathrm{IN}}={\displaystyle \frac{\left(\mathrm{0.0016}+\mathrm{0.0061}{z}_{\mathrm{0}}\right)d_{\mathrm{p}}}{\mathrm{1.414}\times {\mathrm{10}}^{-\mathrm{7}}}},\end{array}$$

$$\begin{array}{}\text{(15)}& R=e^{-\mathrm{2}\sqrt{St}}.\end{array}$$

z₀ is the roughness length, and d_p is particle diameter. The Stokes number is given by

$$\begin{array}{}\text{(16)}& St={\displaystyle \frac{{\mathit{\rho }}_{\mathrm{p}}{d}_{\mathrm{p}}^{\mathrm{2}}}{\mathrm{9}\mathit{\mu }{d}_{\mathrm{c}}}}u.\end{array}$$

u is the horizontal wind velocity and d_c is the diameter of the obstacle.

2.3.3 Z01

In the Z01 scheme (Zhang et al., 2001), the formula of V_d is the same as in the BS95 scheme (Eq. 5). Surface resistance R_s is calculated as

$$\begin{array}{}\text{(17)}& R_{\mathrm{s}}={\displaystyle \frac{\mathrm{1}}{{\mathit{ϵ}}_{\mathrm{0}}{u}_{\ast }\left(E_{\mathrm{B}}+E_{\mathrm{IM}}+E_{\mathrm{IN}}\right)R}},\end{array}$$

$$\begin{array}{}\text{(18)}& E_{\mathrm{B}}=S{c}^{-\mathit{\gamma }},\end{array}$$

where γ depends on land use categories (LUCs) and lies between 0.50 and 0.58.

E_IM is expressed as

$$\begin{array}{}\text{(19)}& E_{\mathrm{IM}}={\left({\displaystyle \frac{St}{\mathit{\alpha }+St}}\right)}^{\mathit{\beta }},\end{array}$$

where β equals to 2. α depends on LUC and lies between 0.6 and 100.0. The Stokes number is given by

$$\begin{array}{}\text{(20)}& St=V_{\mathrm{g}}{u}_{\ast }/gA\end{array}$$

over vegetated surfaces (Slinn, 1982) and

$$\begin{array}{}\text{(21)}& St=V_{\mathrm{g}}{u}_{\ast }^{\mathrm{2}}/g\mathit{\nu }\end{array}$$

over smooth surfaces or surfaces with bluff roughness elements (Giorgi, 1988). E_IN is the collection efficiency based on the relative dimensions of the particle to the collector diameter (Gallagher, 2002). Interception occurs when particles moving with the mean flow and the distance between an obstacle and particle center is less than half of the diameter. Then the particles will collide with and be collected by the obstacle. E_IN is expressed as

$$\begin{array}{}\text{(22)}& E_{\mathrm{IN}}={\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}{\left({\displaystyle \frac{d_{\mathrm{p}}}{A}}\right)}^{\mathrm{2}}\end{array}$$

over vegetated surfaces and E_IN=0 for nonvegetated surfaces, where A is the characteristic radius of collectors. A depends on LUC and lies between 2.0 and 10.0 mm. R is expressed as

$$\begin{array}{}\text{(23)}& R=e^{-\mathrm{1.0}\sqrt{St}}.\end{array}$$

The main differences in formulas used to calculate dry deposition velocity for three different dry deposition schemes are listed in Table 2. For R_s, PE92 and Z01 include the collection efficiency from interception (E_IN) and the rebound effect (R), while these two are neglected in BS95. For the E_IM parameterization, all three schemes use St to parameterize E_IM, but the formulas are quite different. BS95 has a different formula from PE92 and Z01, while the PE92 and Z01 have the same formula but with different coefficients. For PE92, the coefficient for E_IM is constant for all of the surface types. For Z01, the coefficients α and β for E_IM change with different surface types. For the E_IN parameterization, BS95 ignores this effect; PE92 and Z01 use different formulas and variables to calculate E_IN. When large particles (usually >5 µm) hit the nonsticky surface, they are liable to rebound from the surface if they have sufficient kinetic energy. The rebound factor R represents the fraction of particles that stick to the surface (Seinfeld and Pandis, 2006). For rebound effect, BS95 does not consider it; PE92 and Z01 use the same e-exponential form $e^{-b\sqrt{St}}$ to calculate the rebound effect with a different coefficient b. For PE92, b is 2.0; for Z01, b is 1.0. In addition, the parameterization of St is quite different for different dry deposition schemes. For BS95, the formulation of St tends to emphasize the nature of the flow field (Binkowski and Shankar, 1995; Pryor et al., 2008). For Z01, the formulation of St is from Slinn (1982) over vegetated surfaces and from Binkowski and Shankar (1995) over smooth surfaces. The formulation of St from Slinn (1982) and Peters and Eiden (1992) are focus on the individual obstacles (Pryor et al., 2008).

Table 2 Three dry deposition schemes.

2.4 Experiments design

We use WRF-Chem v3.9 with 20 km×20 km horizontal resolution and 35 vertical levels with model top pressure at 50 hPa. The domain covers most of East Asia (14–60^∘ N, 74–130^∘ E) as shown in Fig. 1. The simulation period is from 26 April to 7 May 2017 with time step of 60 s and frequency of output every hour. The time step between radiation physics calls is 20 min. During this period, a severe dust storm event originated from northwestern China and Outer Mongolia, and air quality deteriorated dramatically in a very short time in downwind areas (Guo et al., 2019; Zhang et al., 2018). Meteorological conditions are initialized and forced at the lateral boundaries using the 6-hourly National Center for Environmental Prediction Final (NCEP/FNL) Operational Global Analysis data at a resolution of $\mathrm{1}{}^{\circ }\times \mathrm{1}{}^{\circ }$. For meteorological conditions (such as wind speed and temperature), we reinitialized every 24 h using NCEP/FNL reanalysis data. For chemistry, the output of the aerosol field (such as the concentration of different aerosol species) from the previous 1 d run was used as the initial chemical conditions for the next 1 d run. Our simulation period is from 26 April to 7 May 2017, and one experiment consists of 12 runs, each of 1 d. In this way, the chemical fields are continuous, and we can also get more reliable meteorological conditions. The MOSAIC aerosol scheme was used for all of the simulations. Simulation results prior to 28 April are treated as model spinup for chemical initial condition and are not included in results presented in Sect. 3. The model results from 1 to 7 May are used for the dust loading and concentration analysis. And the model results from 28 April to 7 May are used for the dust emission analysis as the dust emissions before 1 May also have an influence on the dust concentration during 1–7 May. To study the dust simulation sensitivity to dust emission and dry deposition schemes, we run six experiments with two different dust emission schemes and three dry deposition schemes (See Table 3). The corresponding model configuration for dry deposition processes of the six experiments also listed in Table 3. We note here that the USGS LUCs should be selected for the Z01 dry deposition scheme.

Table 3 Model experiments and the corresponding model configuration in WRF-Chem.

2.5 Measurements

2.5.1 PM₁₀

Hourly surface observed PM₁₀ is used to compare with the simulated PM₁₀ from WRF-Chem. In China, hourly surface PM10 concentrations were collected from more than 1000 environmental monitoring stations (locations shown in results section) maintained by the Ministry of Environmental Protection (MEP). The hourly PM₁₀ data from 1 to 7 May 2017 were downloaded from https://quotsoft.net/air/ (last access: 3 April 2020). We colocated the PM₁₀ data with WRF-Chem simulation grids to evaluate model performance with different configurations.

2.5.2 MODIS AOD

Daily aerosol optical depth (AOD) from the Moderate Resolution Imaging Spectroradiometer (MODIS) is used to compare with our simulated AOD from WRF-Chem. The MODIS on board the Aqua satellite was launched by the NASA in 2002, and Aqua is a part of the A-Train satellite constellation. To compare modeled AOD with observations, we use AOD retrievals at 550 nm from MODIS AOD products on Aqua with daily gridded data at a resolution of $\mathrm{1}{}^{\circ }\times \mathrm{1}{}^{\circ }$ (MYD08_D3, Collection 6, combined dark target and deep blue AOD). The MODIS Aqua collection daily MYD08_D3 files were obtained from https://ladsweb.nascom.nasa.gov (last access: 23 April 2020). As Aqua passes through every region of Earth at around 13:30 local time, we extract the model simulation results at 13:00 to compare with the daily MODIS AOD. For the model results, first we divided the domain into different time zones according to the longitude. Then the model results at corresponding UTC when the local time is 13:00 are extracted. The colocated model AOD results for each day are used in comparison with daily MODIS AOD.

2.5.3 CALIPSO data

The vertical profile of aerosol extinction coefficient at a wavelength of 532 nm from the Cloud-Aerosol Lidar and Infrared Pathfinder Satellite observation (CALIPSO) satellite is used to evaluate model results. The CALIPSO launched on 28 April 2006 equipped with CALIOP (Cloud-Aerosol Lidar with Orthogonal Polarization). The CALIOP lidar provides an along-track observation of aerosol and cloud vertical profile. The vertical and horizontal resolutions for CALIOP from the surface to 8.2 km are 30 and 333 m, respectively. Above 8.2 km, the vertical and horizontal resolutions are 60 m and 1 km, respectively. We use the CALIPSO level 2 APro product (V4.20) to obtain the aerosol extinction coefficient (CAL_LID_L2_05kmAPro-Standard-V4-20). The CALIPSO data are available at https://www-calipso.larc.nasa.gov/ (last access: 23 April 2020).

3 Results

3.1 Meteorological conditions

Dust emission and transport processes are closely related to the meteorological conditions. So we first evaluated the model performance in simulating the synoptic conditions. Figure 2 shows the surface meteorological conditions during the dust event. Figure 2a, d, g, and j show the daily mean wind field at 10 m and daily mean temperature at 2 m from NCEP/FNL reanalysis data. The meteorological conditions at 700 hPa are shown in the Supplement (Fig. S1). This dust storm was triggered by the development of a Mongolian cyclone (Figs. S1c and 2d). With the strong northwest and southwest wind near the dust source region, emitted dust was transported to the southeast and northeast of China (Figs. S1c, e and 2d, g). Figure 2b, e, h, and k show the WRF-Chem-simulated daily mean wind field at 10 m and daily mean temperature field at 2 m. Figure. 2c, f, i, and l show the difference in daily mean wind speed at 10 m between WRF-Chem simulation and NCEP/FNL reanalysis data. The WRF-Chem model was able to simulate the wind speed well over the dust source regions (the Taklimakan Desert and the Gobi Desert) and eastern and southern China, where the differences were mostly in the range of −2.0–2.0 m s⁻¹. The wind speed is slightly underestimated near the center of the cyclone (Fig. 2c, f, i, l) and as this is away from dust source regions, we do not expect this underestimation causes large bias in dust emissions. The correlation coefficient (R) and root mean square error (RMSE) between WRF-Chem simulation and FNL reanalysis data for temperature at 2 m, U component of wind, V component of wind, and wind speed at 10 m during simulation period are shown in Table 4. The R for time-averaged temperature at 2 m, U component of wind, V component of wind, and wind speed at 10 m from 1 to 7 May are 1.0, 0.90, 0.86, and 0.82, respectively. The RMSE for time-averaged temperature at 2 m, U component of wind, V component of wind, and wind speed at 10 m from 1 to 7 May are 1.03, 1.08, 0.98, and 1.11, respectively. The R for temperature, U component of wind, V component of wind, and wind speed at 700 hPa from 1 to 7 May are 1.0, 0.94, 0.91, and 0.95, respectively (Table S1 in the Supplement). The RMSE for temperature, U component of wind, V component of wind, and wind speed at 700 hPa from 1 to 7 May are 0.67, 2.34, 2.70, and 1.76, respectively (Table S1). Overall, the correlation coefficients are generally large, and the RMSEs are generally small. This indicates that WRF-Chem performed well in simulating the meteorological conditions. We also compared the difference in the meteorological conditions in our six experiments and found that the difference is negligible (Fig. S2 and Table S2).

Table 4 Correlation coefficient (R) and root mean square error (RMSE) between WRF-Chem simulation and FNL reanalysis data for daily mean temperature at 2 m, U component of wind, V component of wind, and wind speed at 10 m during the dust event time period over the whole domain. The last two rows show the R and RMSE for the time-averaged temperature at 2 m, U component of wind, V component of wind, and wind speed at 10 m from 1 to 7 May.

Figure 2 The left two columns show the surface meteorological conditions during the dust event. The color contours show the daily mean temperature field at 2 m. Vectors represent the daily mean wind field at 10 m. Panels (a), (d), (g), and (j) show the NCEP/FNL reanalysis data. Panels (b), (e), (h), and (k) show the WRF-Chem simulation. The rightmost column shows the difference in daily mean wind speed at 10 m between WRF-Chem simulation and NCEP/FNL reanalysis data. The rectangles show the dust source regions. “TD” is the Taklimakan Desert. “GD” is the Gobi Desert.

3.2 Dust simulation sensitivity to dust emission schemes

In this section, we examine the changes in the simulated dust loading using different dust emission schemes. Figure 3 shows simulated mean dust loading for six experiments over the 7 d simulation period 1–7 May 2017. When using the same dry deposition scheme (BS95, PE92, or Z01), different dust emission schemes give very different dust spatial distributions. Compared with the Shao2011 scheme, GOCART has higher dust loading over the Taklimakan desert (TD) but has relatively lower dust loading over the Gobi Desert (GD), the south of Outer Mongolia, and most parts of northern China. The difference in the spatial distribution of dust loading is mainly caused by the different spatial distributions of dust emission flux from dust emission schemes, as shown in Fig. 4. As the dust emissions before 1 May also have an influence on the dust loading during 1–7 May, the total dust emissions from 28 April to 7 May are analyzed. The total dust emissions from 00:00 UTC on 28 April to 23:00 UTC on 7 May over the GD from GOCART and Shao2011 are 4.90 and 13.88 Tg, respectively. The total dust emissions from 00:00 UTC on 28 April to 23:00 UTC on 7 May over the TD from GOCART and Shao2011 are 7.16 and 2.75 Tg, respectively. Over the GD, the Shao2011 scheme has higher dust emission than GOCART, while over the TD, the GOCART scheme has higher dust emission than Shao2011 (Fig. 4c).

Figure 3 Spatial distribution of simulated mean dust loading for six experiments (a) GOBS95, (b) GOPE92, (c) GOZ01, (d) S11BS95, (e) S11PE92, and (f) S11Z01 over the 7 d simulation period from 00:00 UTC on 1 May to 23:00 UTC on 7 May 2017 (unit: mg m⁻²).

Figure 4 The simulated total dust emission (10⁻³ Tg) from two dust emission schemes: (a) Shao2011 and (b) GOCART from 00:00 UTC on 28 April to 23:00 UTC on 7 May 2017. (c) The total dust emission flux difference between Shao2011 and GOCART. The diameter of the emitted dust is less than 10 µm in both GOCART and Shao2011 dust emission schemes.

Figure 5a, c, e, and g show the spatial distribution of friction velocity, threshold friction velocity, the difference between friction velocity and threshold friction velocity, and the dust emission flux from Shao2011 at 06:00 UTC on 3 May. The areas where the friction velocity is greater than the threshold friction velocity is mainly located in the west of Inner Mongolia and the south of Outer Mongolia (Fig. 5e). This is consistent with Fig. 5g. When the friction velocity is larger than threshold friction velocity, dust can be emitted from the surface. Figure 5b, d, f, and h show the spatial distribution of wind speed at 10 m, threshold velocity, the difference between wind speed at 10 m and threshold velocity, and the dust emission flux from the GOCART dust emission scheme. Different from Shao2011, the dust emission regions from GOCART are not only determined by wind speed but also constrained by the erodibility factor (Eq. 1). From Fig. 5f, the threshold velocity is much smaller than the wind speed at 10 m in most areas. In these areas, GOCART uses Eq. (1) to calculate the dust emission flux, and the source function S depends on the erodibility factor. The dust emission flux in GOCART is directly scaled by erodibility factor. Figure 1 shows the erodibility factor, which describes the fraction of erodible surface in each grid cell. As shown in Fig. 5h, dust emission occurs where the wind speed is high and the erodibility factor is larger than 0.

Figure 5 Spatial distributions of (a) friction velocity (u_∗), (c) threshold friction velocity (u_∗t), and (e) the difference between u_∗ and u_∗t ($u_{\ast }-u_{\ast \mathrm{t}})$ from Shao2011 dust emission scheme at 06:00 UTC on 3 May 2017; (b) wind speed at 10 m (u₁₀), (d) threshold velocity (u_t), and the difference between u₁₀ and u_t (u₁₀−u_t) from the GOCART dust emission scheme at 06:00 UTC on 3 May 2017. Spatial distribution of (g) dust emission flux from Shao2011 and (h) dust emission flux from GOCART at 06:00 UTC on 3 May 2017.

Over the TD, Shao2011 produces lower dust emission flux than GOCART. One reason may be the formula used to calculate the threshold velocity (Eq. 3). The formula used to calculate threshold velocity is from Marticorena and Bergametti (1995), which was originally designed to calculate threshold friction velocity (see LeGrand et al., 2019 for details). This inconsistency leads to a very small threshold velocity in GOCART, which may result in dust emission at low wind speed. Another reason may be the incorrect soil particle size distribution over the TD (Wu and Lin, 2014). The incorrect soil particle size distribution can lead to the unreasonable dust emission flux in Shao2011 over the TD. Over the GD, the GOCART scheme has lower dust emission than the Shao2011 scheme. As mentioned by Su and Fung (2015), the erodibility factor over the GD is highly underestimated and needs to be improved for the GOCART dust emission scheme.

As we mentioned in Sect. 2.2.2, Shao2011 used in this study is based on WRF-Chem v3.9 with some modifications from WRF-Chem v3.7.1. The modified Shao2011 simulates better dust loading than the original Shao2011 scheme in WRF-Chem v3.9 (not shown). Simulated dust emission fluxes are quite different when using two versions of the Shao2011 scheme in WRF-Chem v3.9 and WRF-Chem v3.7.1, which is mainly caused by different soil particle size distributions in two versions. The differences in Shao2011 among different WRF-Chem versions are documented in Appendix B.

3.3 Dust simulation sensitivity to dry deposition schemes

In this section, we analyze dust simulation sensitivity to different dry deposition schemes using the six experiments. For simulated dust loading using the GOCART dust emission scheme (Fig. 3a–c), compared to the BS95 dry deposition scheme, PE92 and Z01 produce higher dust loading over the dust source regions and remote regions. The relative differences in mean dust loading from PE92 and Z01 relative to BS95 are 20 % and 59 %, respectively. As for the simulated dust loading using the Shao2011 dust emission scheme (Fig. 3d–f), PE92 and Z01 schemes also produce higher dust loading than the BS95 scheme, and the differences relative to BS95 are 72 % and 116 %, respectively. This indicates that dust simulation is very sensitive to dry deposition schemes.

Figure 6a shows the modeled dry deposition velocity over desert surface. As desert dust mass is mainly concentrated in the large particle size range, our dry deposition analysis focuses on the coarse mode (2.5–10 µm) (Kok, 2011; Zhao et al., 2013). The reference diameter of the coarse mode is defined at 5 µm (Fig. 6). BS95 produces larger V_d than PE92 and Z01 in the coarse aerosol mode. Larger V_d leads to larger dry deposition and thus lower dust loading, consistent with the lower simulated dust loading from the BS95 scheme discussed above (Fig. 3). In Eq. (5), the dry deposition velocity is comprised of gravitational velocity, aerodynamic resistance, and surface resistance. The diversity of different dry deposition schemes mainly comes from the way to parameterize surface resistance, and differences from gravitational settling and aerodynamics resistance are small (not shown), consistent with previous studies (e.g., Bergametti et al., 2018). Figure 6b shows the surface resistance from different schemes as a function of particle diameter (d_p). In the coarse aerosol mode, Z01 produces the largest surface resistance, followed by PE92 and BS95. Larger surface resistance causes smaller dry deposition velocity in Z01, thus larger dust concentration as shown in Fig. 3.

Figure 6 (a) Dry deposition velocity (V_d), (b) surface resistance (R_s), (c) surface collection efficiency from impaction (E_IM), (d) surface collection efficiency from interception (E_IN), (e) rebound (R), and (f) Stokes number (St) as a function of particle diameter (d_p) over desert surface computed using different dry deposition schemes (BS95, PE92, and Z01). The colored dots indicate values at the reference diameter of 5 µm.

The surface collection efficiency is comprised of Brownian diffusion, impaction, and interception and is corrected for particle rebound (see Eq. 12). Collection from Brownian diffusion is most important for the smaller particles, while collection from impaction and interception play a more important role for large particles in surface collection processes. Figure 6c shows the surface collection efficiency from impaction (E_IM) from different schemes as a function of particle diameter. BS95 gives the largest E_IM, and Z01 gives the smallest. Based on field observation data, Slinn (1982) used a semiempirical fit for smooth surface (Eq. 9), and Binkowski and Shankar (1995) adopted this formula for E_IM and used it for all land surface types. Peters and Eiden (1992) uses Eq. (19) to describe E_IM, with α equals 0.8, and β equals 2 to get the best fit for the data collected over a spruce forest (Eq. 13). In Zhang et al. (2001) scheme, α varies with LUC, and β is chosen as 2 (Eq. 19). For BS95 and PE92, the formula of E_IM is derived from a specific land surface type, but they have been applied to all land surface types in WRF-Chem. This may lead to large uncertainties for dry deposition over the whole domain with different surface types. As the E_IM of Z01 varies with LUC, Z01 may have a better physical treatment of E_IM than the other two dry deposition schemes.

Figure 6d shows the surface collection efficiency from interception (E_IN). E_IN depends on the particle diameter and the characteristic radius of the collectors (Seinfeld and Pandis, 2006). E_IN is important for large particles on hairs at the leaf surface and is negligible over nonvegetated surfaces such as the desert surface we analyzed here (Chamberlain, 1967; Slinn, 1982; Zhang et al., 2001). In BS95, the effect of interception is not considered. In the original PE92 scheme as described in Peters and Eiden (1992), E_IN is also not considered. But in the PE92 scheme used in WRF-Chem, E_IN increases with particle diameter as in Eq. (14). In Z01, the effect of interception is considered as Eq. (22) over vegetated surface and is not considered for nonvegetated surface (as shown in Fig. 6d over desert surface type). The parameterization of E_IN partially results in the difference in surface resistance between PE92 and the other two dry deposition schemes.

Figure 6e shows the rebound factor from different dry deposition schemes. Rebound and resuspension have long been recognized as a mechanism by which the surface can act as sources of particles (Pryor et al., 2008). Due to limited knowledge of particle rebound and resuspension processes, most dry deposition models adopted the form of the rebound effect as $R=e^{-b\sqrt{\mathrm{St}}}$ suggested by Slinn (1982) (Zhang and Shao, 2014; Zhang et al., 2001), while some dry deposition schemes do not include the rebound effect with R=1.0 (Binkowski and Shankar, 1995; Petroff and Zhang, 2010; Zhang and He, 2014). BS95 does not consider the rebound effect. b is equal to 2.0 for the PE92 scheme and 1.0 for the Z01 scheme. Another difference between PE92 and Z01 is the threshold particle diameter for including the rebound effect. Rebound effect is included for PE92 when particles are larger than 0.625 µm and for Z01 when particles are larger than 2.5 µm. In summary, the smaller E_IM and rebound factor lead to larger R_s in Z01, while the larger E_IM leads to smaller R_s in BS95, and the moderate E_IM and rebound effect give a moderate R_s for PE92.

Figure 6f shows the Stokes number from different dry deposition schemes. Over smooth surfaces, the formula of St for BS95 and Z01 is the same, as shown in Eq. (10). In PE92, St is calculated using Eq. (16), which is similar to the formula used in Slinn (1982). BS95 and Z01 schemes give a larger St than PE92. The Stokes number is used to calculate both R and E_IM. The difference in Stokes numbers and the different formulas of R and E_IM lead to the different R and E_IM among different dry deposition schemes (Fig. 6c and e).

Our discussion indicates that Z01 has a better physical treatment of dry deposition velocity, as Z01 considers the rebound effect, and E_IM changes with LUC. The Z01 scheme has also been documented to agree better with measured dry deposition fluxes and dry deposition velocity (e.g., Zhang et al., 2012; Connan et al., 2018). Zhang et al. (2012) compared the dry deposition fluxes measured at five sites in Taiwan with the modeled dry deposition fluxes and found that the measured dry deposition fluxes can be reproduced reasonably well using the Z01 scheme. Connan et al. (2018) conducted experimental campaigns on-site to determine dry deposition velocity of aerosols and found that the Z01 scheme is most suitable for operational use in the size range 0.2–10 µm. All of these indicate that the Z01 dry deposition scheme is more physically meaningful than other two dry deposition schemes.

3.4 Comparisons with observations

To better evaluate the performance of different experiments, we compared the model results with observations. Figure 7 shows hourly observed PM₁₀ concentrations over observational sites at 02:00 UTC on 4 May 2017 (10:00 Beijing Time, BJT, on 4 May 2017). Very high PM₁₀ values (>1000 µg m⁻³) are observed in northern China. Figure 8 compares simulated PM₁₀ in six experiments with observed PM₁₀. During the comparison, the observational sites closest to the model grids are paired up. The correlation coefficients (R), root mean square errors (RMSE) between model and observations, and the mean simulated and observed PM₁₀ for all of the sites over the five regions during the 7 d period 1–7 May are marked in Fig. 8. The simulated PM₁₀ of all of the six experiments have obviously underestimated the observations. Among all of these experiments, GOBS95 has the lowest average PM₁₀ concentration, with a value of 26.45 µg m⁻³, and S11Z01 has the largest one, with a value of 105.17 µg m⁻³, the closest one to the observed mean value of 172.70 µg m⁻³. S11Z01 gives a large R of 0.77 and the smallest RMSE of 96.14 compared to other experiments. Table 5 shows the R and RMSE between the model and observations for PM₁₀ for six experiments over five subregions and over the whole of China. Over the TD, GOBS95 gives the largest R and smallest RMSE. Over the GD, GOZ01 and S11Z01 give a better performance compared with other experiments. For other regions (North China Plain, NCP, Northeast Plain, NEP, and the middle and lower reaches of the Yangtze River plain, YR), S11Z01 gives a relatively larger R and smallest RMSE. For all of the stations in total, S11Z01 gives a larger R of 0.83 and the smallest RMSE of 82.98. Overall, the S11Z01 experiment has the best performance for simulating this dust storm.

Figure 7 Five subregions and observed PM₁₀ concentrations. Key: “1” represents the Taklimakan Desert (TD), “2” represents the Gobi Desert (GD), “3” represents the Northeast Plain (NEP), “4” represents the North China Plain (NCP), and “5” represents the middle and lower reaches of Yangtze River plain (YR). The colored dots represent observed PM₁₀ concentrations over observational sites at 02:00 UTC on 4 May 2017 (10:00 Beijing Time, BJT, on 4 May 2017).

Figure 8 Simulated PM₁₀ versus observed PM₁₀ for six experiments (a) GOBS95, (b) GOPE92, (c) GOZ01, (d) S11BS95, (e) S11PE92, and (f) S11Z01 over the 7 d simulation period 1–7 May 2017. “Obs mean” is mean PM₁₀ over 1–7 May from observation; “Model mean” is mean PM₁₀ over 1–7 May from simulation; “R” is the correlation coefficient between model and observations; “RMSE” is the root mean square error. Different color dots represent different regions as shown in Fig. 6.

Figure 9 Simulated and observed mean AOD over the simulation period 1–5 May. Panels (a)–(e) show the distribution of daily mean aerosol optical depth (AOD) at 550 nm derived from MODIS-Aqua. Panels (f)–(j) show the WRF-Chem-simulated AOD for the GOBS95 experiment. Panels (k)–(n) show the WRF-Chem-simulated AOD for the S11Z01 experiment. The model results are extracted from the simulation results at 13:00 local time for each region to match the MODIS observation time (for details see Sect. 2.5). All of the other experiments and for the period 1–6 May are shown in Fig. S1. Grid points without valid MODIS AOD retrieval are masked for both observational and model results.

Table 5 Correlation coefficient (R) and root mean square error (RMSE) between the model and observations for PM₁₀ over five subregions and for all of the stations over whole China in Fig. 7 for six experiments listed in Table 3.

Figure 9 shows the MODIS-observed daily mean AOD and WRF-Chem-simulated AOD over the simulation period 1–5 May. For strong dust storms like the one we examined here, dust particles contribute the most to AOD, and AOD therefore can represent the dust loading in the atmosphere. To match the MODIS AOD observation time, simulated AOD at 13:00 local time is used for comparison (see Sect. 2.5 for details). For each $\mathrm{1}{}^{\circ }\times \mathrm{1}{}^{\circ }$ grid with observed AOD from MODIS, the average value of simulated AOD from WRF-Chem in this grid is calculated. Grid points without valid MODIS AOD retrieval are masked for both observational and model results in Fig. 9. A major dust emission event occurred over the GD on 3 May (Fig. 9c). Shao2011 simulated well the dust emission event over the GD on 3 May (Fig. 9m), while GOCART obviously underestimated dust emission over the GD (Fig. 9h). On 4 May, emitted dust from the GD was transported to northeast China, and the highest AOD values for this case study were observed in northern China (Fig. 9d). As the GD is the main dust source region of this dust storm, Shao2011 correctly captured the emission phase of this dust event. Simulated AOD values from the S11Z01 configuration produced the closest match to the observed daily MODIS AOD with respect to the magnitude and spatial pattern (Figs. 9n and S3). For a more quantitative comparison, Table 6 shows the correlation coefficient (R) and root mean square error (RMSE) between the model and observed AOD for six experiments during 1–7 May. Overall, the S11Z01 experiment gives a larger correlation coefficient, and the RMSE is almost the same among different experiments; the correlation coefficient is still lower than 0.5. The low correlation may partly come from the spatial and temporal limitation of satellites and the difficulties to retrieve aerosol in the vicinity of clouds for satellites.

Table 6 Correlation coefficient (R) and root mean square error (RMSE) between the model and MODIS observation for AOD for six experiments over the 7 d simulation period 1–7 May 2017.

To evaluate the model performance in simulating the vertical profile of dust aerosol, we compared the extinction coefficient from the model and from CALIPSO (Fig. 10). Figure 10 shows the simulated and observed aerosol extinction profiles at 532 nm at 18:00 UTC on 4 May. The trajectory of CALIPSO passes East Asia (Fig. 10d). All of the six experiments show the similar dust location in the atmosphere, which is consistent with the CALIPSO observation. However, the magnitude of dust concentration differs substantially. The simulated extinction coefficients using GOCART dust emission schemes are significantly underestimated compared to the CALIPSO observation (Fig. 10a, b, and c), while the modeled extinction coefficients using Shao2011 dust emission scheme agrees better with observation though they are still underestimated (Fig. 10e, f and g). Among all of the six experiments, results from S11Z01 agree the best with observation.

Figure 10 Simulated and observed aerosol extinction profiles at 532 nm at 18:00 UTC on 4 May. Panel (d) shows the CALIPSO trajectory. Panel (h) shows the CALIPSO observed extinction coefficient. Panels (a)–(g) show the WRF-Chem-simulated extinction coefficient from six experiments.

In summary, both ground and satellite observations indicate that the S11Z01 experiment yields the best performance in simulating this dust storm. As we discussed in Sect. 3.2, the Z01 dry deposition scheme indeed has a better physical treatment and performs better than some other dry deposition schemes.

4 Summary and discussion

In this study, we analyzed the dust simulation sensitivity to different dust emission schemes and dry deposition schemes. In order to compare different dust emission schemes, the Shao2011 dust emission scheme has been implemented into the MOSAIC aerosol scheme in WRF-Chem v3.9. Six model experiments were conducted to simulate the dust storm in May 2017 over East Asia, with two dust emission schemes (GOCART and Shao2011) and three dry deposition schemes (BS95, PE92, and Z01). The simulation results of different experiments were evaluated against surface and satellite observations.

Our results show that dust loading is very sensitive to different dry deposition schemes. The relative difference in dust loading in different experiments range from 20 % to 116 % when using different dry deposition schemes. The difference in dry deposition velocity in different dry deposition schemes comes from the parameterization of surface resistance, and difference in surface resistance mainly comes from the parameterization of collection efficiency from impaction and rebound effect. In addition, different dust emission schemes result in different spatial distributions of dust loading, as dust emission fluxes in dust source regions differ substantially among different dust emission schemes, which is mainly attributed to differences in the threshold conditions for dust emission and in formulas and parameters for calculating dust emission flux. We noted that the Shao2011 dust emission scheme is different among different WRF-Chem versions, and a significant difference exists in the simulated dust emission fluxes between WRF-Chem v3.9 and WRF-Chem v3.7.1, which is mainly caused by differences in soil particle size distributions used in two versions (see Appendix B).

Compared with both surface PM₁₀ station observations and MODIS AOD, the Shao2011 dust emission scheme coupled with the Z01 dry deposition scheme produces the best simulation for the dust storm in East Asia. For PM₁₀, the S11Z01 experiment gives a larger R of 0.83 and the smallest RMSE of 82.98 of all of the stations (Table 5). The spatial distribution of AOD during the simulation period obtained by S11Z01 agrees the best with MODIS AOD (Fig. 9), with the largest R and a relatively small RMSE (Table 6). Our analysis indicates Z01 accounts for the rebound effect, and E_IM changes with LUC and therefore has a better physical treatment of dry deposition velocity than the two other dry deposition schemes. Previous studies have also shown that the Z01 scheme agrees better with measured dry deposition fluxes and dry deposition velocity (e.g., Zhang et al., 2012; Connan et al., 2018). The Shao2011 dust emission scheme has larger dust emission fluxes than the GOCART dust emission scheme over the Gobi Desert, and the transport of dust emitted from the Gobi Desert is the most important source of dust weather in northern China (Chen et al., 2017b). Compared with daily MODIS AOD (Fig. 9), our results indicate that dust emission from Shao2011 is better for this dust event, in terms of dust spatial and temporal distributions. We note that our results are obtained from simulations of a dust storm over a short period, and longer simulations are desirable in the future to test whether the optimal scheme here still produces best simulations.

This study highlights the importance of dry deposition process in dust simulation. Future studies on dust simulation should pay attention to improve dry deposition schemes as well as the dust emission schemes. Additional field measurements of dry deposition process and comparisons with model results are required to reduce the uncertainties in dust simulation.

Appendix A: Description of the Shao2011 dust emission scheme

Here we describe the Shao2011 dust emission scheme in more detail as a supplement to the Sect. 2.2.2 of this article. The total saltation flux Q in Eq. (4) is calculated as

$$\begin{array}{}\text{(A1)}& Q=\underset{d_{\mathrm{1}}}{\overset{d_{\mathrm{2}}}{\int }}Q\left(d\right)p_{\mathrm{m}}\left(d\right)\mathit{\delta }d,\end{array}$$

where d₁ and d₂ define the upper and lower limits of saltation particle size. p_m(d) is the minimally disturbed particle-size distribution. The saltation flux Q for each particle size d is calculated as

$$\begin{array}{}\text{(A2)}& Q\left(d\right)=\left(\mathrm{1}-c_{\mathrm{f}}\right)c_{\mathrm{0}}{\displaystyle \frac{{\mathit{\rho }}_{\mathrm{a}}}{g}}{u}_{\ast }^{\mathrm{3}}\left(\mathrm{1}-{\displaystyle \frac{u_{\ast \mathrm{t}}}{u_{\ast }}}\right){\left(\mathrm{1}+{\displaystyle \frac{u_{\ast \mathrm{t}}}{u_{\ast }}}\right)}^{\mathrm{2}},\end{array}$$

where c_f is the fraction of vegetation cover, (1−c_f) means the fraction of erodible surface area, and c₀ is a coefficient. u_∗t is the threshold friction velocity, and u_∗ is the friction velocity. When u_∗ is larger than u_∗t, it calculates the dust emission flux. Before WRF-Chem v4.0, there is a bug in calculating the saltation flux Q(d_s) in Shao2011. The last term is miscalculated as 1$+{\left(\frac{u_{\ast \mathrm{t}}}{u_{\ast }}\right)}^{\mathrm{2}}$ in WRF-Chem codes (LeGrand et al., 2019). In WRF-Chem v4.0 and later versions, this bug is fixed, and we also fixed this bug in our simulations.

The threshold friction velocity u_∗t is calculated as

$$\begin{array}{}\text{(A3)}& u_{\ast \mathrm{t}}=u_{\ast t\mathrm{0}}{f}_{\mathit{\lambda }}{f}_{\mathit{\theta }},\end{array}$$

where u_∗t0 is the ideal threshold friction velocity when soil is dry, bare, and free of crust and salt, f_λ is the correction functions for surface roughness, and f_θ is the correction functions for soil moisture. The ideal threshold friction u_∗t0 is calculated as

$$\begin{array}{}\text{(A4)}& u_{\ast \mathrm{t}\mathrm{0}}=\sqrt{a_{\mathrm{1}}{\displaystyle \frac{{\mathit{\rho }}_{\mathrm{p}}}{{\mathit{\rho }}_{\mathrm{a}}}}gd+{\displaystyle \frac{a_{\mathrm{2}}}{{\mathit{\rho }}_{\mathrm{a}}d}}},\end{array}$$

where a₁ and a₂ are constant. ρ_p and ρ_a are particle and air density. d is the particle diameter.

The correction functions for surface roughness f_λ are calculated as

$$\begin{array}{}\text{(A5)}& f_{\mathit{\lambda }}={\left[\left(\mathrm{1}-m\mathit{\sigma }\mathit{\lambda }\right)\left(\mathrm{1}+m{\mathit{\beta }}_{\mathrm{0}}\mathit{\lambda }\right)\right]}^{\frac{\mathrm{1}}{\mathrm{2}}},\end{array}$$

where m is a constant, σ is the ratio of roughness-element basal area to frontal area, λ is the frontal area index, and β₀ is the ratio of the drag coefficient of an isolated roughness element on the surface to the drag coefficient of the substrate surface itself.

The mass fraction of free dust η_mi is calculated as

$$\begin{array}{}\text{(A6)}& {\mathit{\eta }}_{\mathrm{mi}}=\underset{d-\frac{\mathrm{\Delta }{d}_i}{\mathrm{2}}}{\overset{d+\frac{\mathrm{\Delta }{d}_i}{\mathrm{2}}}{\int }}{p}_{\mathrm{m}}\left(d\right)\mathit{\delta }d,\end{array}$$

where p_m(d) is the minimally disturbed particle-size distribution, which is regarded as a composite of several lognormal distributions, and p_m(d) is expressed as

$$\begin{array}{}\text{(A7)}& p_{\mathrm{m}}\left(d\right)={\displaystyle \frac{\mathrm{1}}{d}}\sum _{i=\mathrm{1}}^J{\displaystyle \frac{w_j}{\sqrt{\mathrm{2}\mathit{\pi }}{\mathit{\sigma }}_j}}\exp (-{\displaystyle \frac{{\left(\ln d-\ln D_j\right)}^{\mathrm{2}}}{\mathrm{2}{\mathit{\sigma }}_j^{\mathrm{2}}}}).\end{array}$$

Soil samples collected from experiment sites are used to determine the particle-size distribution.

Appendix B: The Shao2011 dust emission scheme in different versions of WRF-Chem

As we noted in Sect. 3.2, the Shao2011 scheme in different versions of WRF-Chem can produce significantly different dust emission fluxes. Here we document differences in Shao2011 among different WRF-Chem versions.

Figure B1 Free dust fraction for 12 soil types as a function of particle diameter (d_p). The red lines represent the free dust fraction in WRF-Chem v3.8 and later versions. The blue lines represent the free dust fraction before WRF-Chem v3.8. The colors of the soil type font in the upper left corner of the plot are different. In WRF-Chem v3.8 and later versions, each soil type has a corresponding free dust fraction distribution. In versions before WRF-Chem v3.8, several soil types share a free dust fraction distribution. The same soil type font color indicates that a free dust fraction is shared among these soil types in versions before WRF-Chem v3.8.

The first difference is c₀, where c₀ is a coefficient used to calculate the saltation flux as in Eq. (A2). In versions before WRF-Chem v3.8, c₀ is equal to 0.5; in WRF-Chem v3.8 and later versions, c₀ is equal to 2.3 (Table B1).

The second difference is β₀. β₀ is a coefficient used to calculate the correction function for surface roughness f_λ in Eq. (A5). In versions before WRF-Chem v3.8, β₀ is 90; in WRF-Chem v3.8 and later versions, β₀ is 200 (Table B1).

The third difference is caused by the minimally disturbed particle-size distribution p_m(d) (see Eq. A6). p_m(d) is used to calculate the free dust fraction η_mi (see Eq. A6). Free dust fraction is the fraction of dust that has low enough binding energy so that it can be easily lifted from the surface by either aerodynamic forces or mechanical abrasion (Shao, 2001). The η_mi is used to calculate the dust emission rate in Eq. (4). Twelve soil types are included in all WRF-Chem versions. In WRF-Chem v3.8 and later versions, each soil type has a corresponding p_m(d) as listed in Table B1 from Shao et al. (2010); in versions before WRF-Chem v3.8, there are only four p_m(d) as listed in Table 1 from Shao (2004) for 12 soil types (Fig. B1). For example, (f) sand and (g) loamy sand soil types use the same free dust fraction distribution in versions before WRF-Chem v3.8. As shown in Table B2, the loam and clay loam are the two soil types with the largest percentage, while the other soil types account for a very small percentage. From Fig. B1c and e, for loam and clay loam soil types, the free dust fraction is very small in the particle size range 0–10 um in WRF-Chem v3.8 and later versions, almost all close to 0; however, in the versions before WRF-Chem v3.8, the free dust fraction is relatively high. In different WRF-Chem versions, the total saltation flux Q is the same, but dust emission flux F(d_i) is different due to different free dust fractions (see Eq. 4). With a smaller free dust fraction, the dust emission flux is smaller in WRF-Chem v3.8 and later versions.

To examine the importance of these changes, we run four experiments to quantify the contribution of each factor (Table B3). For the control run, c₀ is 2.3, β₀ is 200, and p_m(d) has 12 distributions based on WRF-Chem v3.8 or later versions. For the case 1 experiment, β₀ is changed to 90, the one used in WRF-Chem v3.7.1 and all other parameters are kept the same as in the control run. The dust emission in case 1 is 1.35 times higher than the control run. For the case 2 experiment, c₀ is changed to 0.5, the one used in WRF-Chem v3.7.1, and all other parameters remain the same as in the control run. The dust emission in case 2 is 21 % of the dust emission of the control run. For the case 3 experiment, p_m(d) is adopted from WRF-Chem v3.7.1 and has four distributions, and all other parameters remain the same as in the control run. The dust emission in case 3 is 13 times higher than the control run. This indicates that the difference in dust emission between different versions of Shao2011 scheme is mainly caused by the change in p_m(d). As p_m(d) is determined by soil particle size distribution, this also highlights the need to improve the accuracy of soil texture.

Table B1 Differences in Shao2011 dust emission scheme between different WRF-Chem versions.

Table B2 Percentage of each soil type in the whole East Asia domain.

Table B3 Sensitivity of the simulated total dust emission from the Shao2011 to model parameters over the 7 d simulation period 1–7 May 2017. The multiple of the dust emission of different cases is calculated with respect to the control run.

We should mention that the Shao2011 dust emission scheme we used in this study is based on WRF-Chem v3.9 with the soil particle size distribution from WRF-Chem v3.7.1, which simulates better dust loading compared with observations. Compared with the original Shao2011 scheme in WRF-Chem v3.9, the total dust emission simulated in our experiments during 1–7 May is 13 times higher.