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.


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 Kok et al., 2017;Zhao et al., 2010Zhao et al., , 2011Zhao et al., , 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 Published by Copernicus Publications on behalf of the European Geosciences Union.
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 dif-ferent dust emission schemes (Kang et al., 2011;LeGrand et al., 2019;Su and Fung, 2015;Lin, 2013, 2014;Yuan et al., 2019;Zhao et al., 2010Zhao et al., , 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;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 Y. Zeng et al.: WRF-Chem simulations of the East Asian dust storm 2127  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

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 4bin 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.

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(Shao, , 2004 (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 Su and Fung, 2015;Wu and Lin, 2014), and Shao2011  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.

GOCART
The formula of vertical dust flux in GOCART is approximated as 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 10 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 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).   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 1.0 × 10 −9 kg s 2 m −5 , which is consistent with the original GOCART dust emission scheme paper (Ginoux et al., 2001). For dust_opt=1, C value is set to 0.8 × 10 −9 kg s 2 m −5 .

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 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.

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 MO-SAIC 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. In the BS95 scheme (Binkowski and Shankar, 1995), V d is expressed as where R a and R s are aerodynamic and surface resistance; V g is the gravitational settling velocity and is given as 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 where E B is collection efficiency from Brownian diffusion. E B is calculated as follows: where Sc is the Schmidt number, given by Sc = ν/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 where St is the Stokes number, given by 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).

PE92
In PE92 scheme (Peters and Eiden, 1992), the dry deposition velocity (V d ) is expressed as 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 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 z 0 is the roughness length, and d p is particle diameter. The Stokes number is given by u is the horizontal wind velocity and d c is the diameter of the obstacle.

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 where γ depends on land use categories (LUCs) and lies between 0.50 and 0.58. E IM is expressed as where β equals to 2. α depends on LUC and lies between 0.6 and 100.0. The Stokes number is given by over vegetated surfaces (Slinn, 1982) and 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.

Y. Zeng et al.: WRF-Chem simulations of the East Asian dust storm
Then the particles will collide with and be collected by the obstacle. E IN is expressed as 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 The main differences in formulas used to calculate dry deposition velocity for three different dry deposition schemes are listed in Table 2 (Seinfeld and Pandis, 2006). For rebound effect, BS95 does not consider it; PE92 and Z01 use the same e-exponential form e −b √ 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).

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 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 1 • × 1 • . 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.

PM 10
Hourly surface observed PM 10 is used to compare with the simulated PM 10 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 10 data from 1 to 7 May 2017 were downloaded from https://quotsoft.net/air/ (last access: 3 April 2020). We colocated the PM 10 data with WRF-Chem simulation grids to evaluate model performance with different configurations.

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 1 • ×1 • (MYD08_D3, Collection 6, combined dark target and deep  3 Results

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 −1 . The wind speed is slightly underestimated near the center of the cyclone (Fig. 2c, f, (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).

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 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.
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.

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.
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 √ 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.

Comparisons with observations
To better evaluate the performance of different experiments, we compared the model results with observations. Figure 7 shows hourly observed PM 10 concentrations over observational sites at 02:00 UTC on 4 May 2017 (10:00 Beijing Time, BJT, on 4 May 2017). Very high PM 10 values (> 1000 µg m −3 ) are observed in northern China. Figure 8 compares simulated PM 10 in six experiments with observed PM 10 . 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 10 for all of the sites over the five regions during the 7 d pe- 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 * − u * t ) from Shao2011 dust emission scheme at 06:00 UTC on 3 May 2017; (b) wind speed at 10 m (u 10 ), (d) threshold velocity (u t ), and the difference between u 10 and u t (u 10 − 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. riod 1-7 May are marked in Fig. 8. The simulated PM 10 of all of the six experiments have obviously underestimated the observations. Among all of these experiments, GOBS95 has the lowest average PM 10 concentration, with a value of 26.45 µg m −3 , and S11Z01 has the largest one, with a value of 105.17 µg m −3 , the closest one to the observed mean value of 172.70 µg m −3 . 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 10 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 ex-periments. 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 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   , (e) S11PE92, and (f) S11Z01 over the 7 d simulation period 1-7 May 2017. "Obs mean" is mean PM 10 over 1-7 May from observation; "Model mean" is mean PM 10 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. at 13:00 local time is used for comparison (see Sect. 2.5 for details). For each 1 • × 1 • 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 Table 5. Correlation coefficient (R) and root mean square error (RMSE) between the model and observations for PM 10 over five subregions and for all of the stations over whole China in Fig. 7 for six experiments listed in Table 3.
Region R/RMSE GOBS95 GOPE92 GOZ01 S11BS95 S11PE92 S11Z01 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. 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.
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.

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 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. GOBS95 GOPE92 GOZ01 S11BS95 S11PE92 S11Z01 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 10 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 10 , 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. 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.
The first difference is c 0 , where c 0 is a coefficient used to calculate the saltation flux as in Eq. (A2). In versions before WRF-Chem v3.8, c 0 is equal to 0.5; in WRF-Chem v3.8 and later versions, c 0 is equal to 2.3 (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 0 is 2.3, β 0 is 200, and p m (d) has 12 distributions based on WRF-Chem v3.8 or later versions. For the case 1 experiment, β 0 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 0 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.
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. 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. Data availability. The 6-hourly National Center for Environmental Prediction Final (NCEP/FNL) Operational Global Analysis data at a resolution of 1 • × 1 • can be obtained from https://doi.org/10.5065/D6M043C6 (NCEP/NWS/NOAA/U.S. Department of Commerce, 2000). The observed PM 10 data are collected from the national air quality real time release platform at http://106.37.208.233:20035/ (last access: 23 April 2020; China National Environmental Monitoring Centre, 2020a). The historical data of air quality used in this study can be downloaded from https://quotsoft.net/air/ (last access: 23 April 2020; China National Environmental Monitoring Centre, 2020b). Daily MYD08_D3 files from the MODIS on board the Aqua satellite can be obtained from https://doi.org/10.5067/MODIS/MOD08_D3.006 (Platnick et al., 2015).
Author contributions. MW and YZ conceived the idea and designed the model experiments. YZ performed the simulations, conducted the analysis, and wrote the paper. CZ and SC provided guidelines for the dust simulations in WRF-Chem and helped in the interpretation of the results. XH and ZL provided guidelines for the WRF-Chem simulations and contributed to the model experiments design. YG helped in the PM 10 data processing and usage. Everyone edited the paper.
Competing interests. The authors declare that they have no conflict of interest.