Simulation of a Severe Sand and Dust Storm Event in March 2021 in Northern China: Dust Emission Schemes Comparison and the Role of Gusty Wind

: Northern China experienced a severe sand and dust storm (SDS) on 14/15 March 2021. It was difﬁcult to simulate this severe SDS event accurately. This study compared the performances of three dust-emission schemes on simulating PM 10 concentration during this SDS event by implementing three vertical dust ﬂux parameterizations in the Comprehensive Air-Quality Model with Extensions (CAMx) model. Additionally, a statistical gusty-wind model was implemented in the dust-emission scheme, and it was used to quantify the gusty-wind contribution to dust emissions and peak PM 10 concentration. As a result, the LS scheme (Lu and Shao 1999) produced the minimum errors for peak PM 10 concentrations, the MB scheme (Marticorena and Bergametti 1995) underestimated the PM 10 concentrations by 70–90%, and the KOK scheme (Kok et al. 2014) overestimated PM 10 concentrations by 10–50% in most areas. The gusty-wind model could reasonably reproduce the probability density function of 2-min wind speeds. There were 5–40% more dust-emission ﬂux and 5–40% more peak PM 10 concentrations generated by the gusty wind than the hourly wind in the dust-source regions. The increase of peak PM 10 concentration caused by gusty wind in the non-dust-source regions was higher than in the dust-source regions, with 10–50%. Implementing the gusty-wind model could help improve the LS scheme’s performance in simulating PM 10 concentrations of this severe SDS event. More work is still needed to investigate the reliability of the gusty-wind model and LS scheme on various SDS events.


Introduction
Sand and dust storms (SDSs) are extremely hazardous weather systems characterized by low visibility, strong wind, and high particle-matter concentration that frequently occur in arid and semi-arid regions such as East Asia [1][2][3].There were about ten SDSs events per year in East Asia in the last 20 years, and the variations of the frequency and intensity SDSs events were closely related to surface conditions and the intensity of the polar vortex [2,4].SDSs events not only cause loss of life and agriculture destruction but also influence weather and climate by disturbing radiation balance [1,[5][6][7].SDSs were shown to be an important source of PM 10 pollution in East Asia [8,9].Many methods, including observations and numerical models, have been used to comprehend SDSs processes [7,[10][11][12][13][14].
Studies have been performed to improve the model performance on the simulation of SDS events.Parameterization of dust-emission flux is one of the main challenges in dustmodel development [10].Several dust-emission schemes have been developed to estimate dust emission.These dust-emission schemes can reasonably well simulate the SDSs events.However, the dust emissions vary by several orders of magnitude when using different dust-emission schemes [10,[15][16][17].Kang et al. [15] reported that the vertical dust fluxes differed greatly under the same horizontal dust flux condition because the sandblasting efficiency varied among these dust-emission schemes.Other studies also indicated that sandblasting efficiency was crucial for simulating dust-emission fluxes [16,17].
Gusty wind also plays an important role in modeling the dust-emission flux.The large variations in gusty-wind speed and the coherent structure of gusty wind could generate more dust-emission flux than the steady wind [18,19].Wang and Zheng [20] predicted the sand-transport rates, respectively, under steady and unsteady wind, using a one-dimension model.Results showed that the sand-transport rates by unsteady wind could be 10-200% higher than the steady wind [20].However, the gusty wind was neglected and simplified as the hourly wind in most dust-emission schemes [11].Studies about gustywind parameterization [21][22][23] and the probability density function of gusty wind [24] made it possible to consider the gusty wind in dust-emission schemes.
Northern China experienced a severe SDS in March 2021, leading to excessive particulate pollution and strong wind [25][26][27].Yin et al. [26] analyzed the climatological causes of this SDS event.They concluded that the strong gusty wind in the dust-source area corresponded to the increases in the PM 10 concentration in Beijing, China.This study implements three vertical dust-flux parameterizations in the Comprehensive Air-Quality Model with Extensions (CAMx) to investigate their performances on this severe SDS event.
A gusty-wind prediction model was developed to simulate the probability density function of gusty-wind speed and implemented in the dust-emission scheme to quantify the role of gusty wind in this SDS event.Atmospheric rivers intertwined strongly with exceptional SDS events [28,29]; hence, we also analyzed the possible influence of atmospheric rivers on this exceptional SDS event.

Data
The weather phenomena, 10 m wind speed, 2 m temperature, visibility, and PM 10 concentrations were obtained from China Meteorological Administration.Two-minute average 10 m wind speed and max instantaneous 10 m wind speed were used to analyze the probability density function of gusty wind in Northern China.The instantaneous wind speed was a 3 s-sample wind speed.Hourly PM 10 concentrations in six sites, including three sites near the dust-source regions and three sites in the transport pathway, were used to test the model performance on dust emission and transportation.PM 10 concentrations over 6500 µg/m 3 were recorded as 6500 µg/m 3 , due to the limitation of measurements in some stations, such as Hohhot and Yinchuan.The bias, root mean square error (RMSE), and the index of agreement (IOA), which were explained in detail in Baker et al. [30], were used to evaluate the simulation of wind and PM 10 concentrations.Moderate-resolution imaging spectroradiometer (MODIS) aerosol optical depth (AOD) data, at 550 nm from the Deep Blue algorithm over land, was also used to evaluate the simulated AOD.

Model Configuration
The Comprehensive Air-Quality Model with Extensions (CAMx) was used to model the transport and diffusion of particle matters [31].The CAMx model was an off-line chemistry transport model.The Weather Research and Forecasting model (WRF 4.0) [32] was used to generate meteorological fields for CAMx and the dust-emission model.The first 24 h were disregarded as spin-up for each SDSs episode modeling.Boundary and initial conditions were derived from the National Centers for Environmental Prediction (NCEP) final (FNL) global reanalysis data, with the resolution of 1 • .The physical parameterization options chosen for the WRF modeling included Dudhia shortwave radiation, RRTM longwave scheme, YSU planetary-boundary layer scheme, and the Unified Noah land-surface scheme.The model domain covers most of East Asia.The model has a horizontal grid spacing of 24 km, with 30 vertical levels from the surface to 50 hpa.Figure 1 shows the model domain.

Dust-Emission Schemes
A dust-emission scheme was used to generate windblown dust emissions for CAMx, based on a revised approach used in the global ECHAM/MESSy atmospheric chemistryclimate model [33,34].The scheme calculates dust-emission flux based on surface friction velocity and threshold surface friction-velocity.The threshold friction-velocity depends on the size of soil particles, air density, and soil moisture, which was explained in detail in our previous study [9].Three vertical dust-emission parameterizations were implemented in the dust-emission scheme to perform an inter-comparison of different dust-emission parameterizations.
The first vertical dust-flux parameterization (MB scheme) was proposed by Marticorena and Bergametti [35].The relationship between vertical and horizontal dust fluxes is related to the clay content of soils less than 20%.
where F is the vertical dust emission flux, Q is the horizontal dust flux, and f c is the clay content in %.
The second vertical dust-flux parameterization (LS scheme) was proposed by Lu and Shao [5].In the LS scheme, the saltation bombardment efficiency ratio F/Q is linearly proportional to the fraction of dust contained in the soil and inversely proportional to the soil surface-hardness parameter.
where C α and C β are constant values; g is the gravitational acceleration; p is the plastic pressure of soil surface in N/m 2 ; f s is the sand fraction in soil; ρ b and ρ p are the bulk-soil and soil-particle density, respectively; and u * is surface friction velocity.The bulk-soil density in our model refers to Liu et al. [25], and other parameters in the LS scheme refer to Kang et al. [15].The third vertical dust-flux parameterization (KOK scheme) was defined in the study of Kok et al. [11].In the KOK scheme, soil erodibility is a critical parameter for calculating dust-emission flux.
where C d0 , C e , and C α are dimensionless coefficients; ρ a is the air density; u * t is the threshold friction velocity; u * st is the standardized threshold friction velocity, which is defined as the value of u * t at standard atmospheric density at sea level (1.225 kg/m 3 ); u * st0 is standardized threshold friction velocity for the highly erodible soil, which is set to 0.16 m/s; and f bare is the fraction of the surface that consists of bare, erodible soil.The parameters setting refers to the study of Kok et al. [11].
The land cover, land-surface information, and dust size distribution in the model can be found in our previous study [9].We implement these three schemes in the CAMx wind-blowing dust scheme to simulate the SDSs.

Gusty Wind Scheme
A statistical model for predicting wind-speed probabilities [24] was used to characterize the gusty wind during the SDS episode.Efthimiou et al. [24] reported that the statistical behavior of near-surface wind speeds could be adequately represented by the Beta distribution, which was even slightly better than the widely used Weibull distribution.The probability density function (PDF) for the time-averaged horizontal wind speed is given as follows: where V is the wind speed, V max is the max wind speed over a period of time, V a is the mean wind speed over a period of time, I is the wind speed fluctuation intensity, and σ 2 V is variance of wind speed.
According to Harper et al. [36] and Suomi et al. [37], the standard deviation of the wind speed can be estimated as follows: A simple parameterization is used to parameterize V max [38] as follows: where a σ and b g are the two factors (normalized gust wind speed and gust factor) for predicting the standard deviation of wind speed and max wind speed by using the mean wind speed.These two factors can be obtained from the statistics of the observed wind speed.The processes mentioned above allow us to use the WRF output hourly wind speed and the statistical factors to parameterize the distribution of gusts.We implement this gusty-wind model in the dust-emission scheme to quantify the contribution of gusty wind to dust-emission flux.

Sand and Dust Storm Event Description
In Mongolia and Northern China, a super-strong Mongolian cyclone (min SLP = 978 hPa, Figure 2a) triggered a severe SDS event on 14/15 March 2021 [26].The minimum visibility was much lower than 0.5 km, and the wind speed exceeded 15 m/s in Mongolia and Northwestern China during this SDS episode (Figure 2c).The sand dust particles were transported to most areas to the north of the Yangtze River by the north wind and caused low visibility (<7 km, Figure 2c) and high PM 10 concentrations [26].It was the strongest SDS event in the recent decade [26].The destructive cooling and warming in early and late winter led to bare and loose land surfaces in dust-source areas in Mongolia.Temperatures in Inner Mongolia were significantly warmer than climatic conditions by 5-8 degrees Celsius in the first half of March.Lack of precipitation, excessive snowmelt, and strong evaporation resulted in dry soil.The loose and dry soil in Mongolia provided favorable conditions in dust sources [26,27].The atmospheric river (Figure 2b) from the South China Sea to Northeast China may enhance the intensity of the Mongolian cyclone by providing more water vapor for latent heat release [39,40].When the super-strong Mongolian cyclone developed, the strongest SDS event occurred in Mongolia and China.

Characteristics of Gusty Wind during the Sand and Dust Storm Event
There were strong gusts during this SDS episode that were caused by northerly advection, downward transport of westerly momentum, and strong instability [26].The instantaneous wind speed exceeded 20 m/s in Inner Mongolia along the China-Mongolia border and reached 15 m/s in Beijing (Figure 3a).The gusty-wind speeds were also greater than 20 m/s in Mongolia on 14 March 2021 [26].The strong gusty winds shook, blew, lifted, and transported sand and dust particles into North China.We used the 2-min averaged wind speed to calculate the standard deviation and max wind speed over a 1-h period.The normalized gusty-wind speeds (defined in Equation ( 8)) ranged from 2.3 to 2.7 at most stations over Northern China (Figure 3c).In other words, the gust-wind speeds tended to be about 2.5 standard deviations away from the average wind speed.The correlation coefficients of the standard deviation (σ V ) with the difference be- tween the max and the mean 2-min wind speed were around 0.9 (Figure 3d), which lent high confidence in estimating the standard deviation by using normalized gust-wind speeds.
The gust factor is another parameter for the gusty-wind model.The max 2-min wind speed can be estimated by Equation (9).The gust factor ranged from 1.2 to 1.3 for stations in Inner Mongolia (Figure 3e).The correlation coefficients were greater than 0.9 (Figure 3f) between hourly max 2-min wind speed and hourly mean wind speed.The max 2-min wind speed could be used to estimate the standard deviation.However, it could not be used in Equation (6) to estimate the wind PDF, because the max 2-min wind speed was a measured value, rather than an extreme value as a direct consequence of statistical uncertainties associated with limited measuring times [24].Therefore, we used the max instantaneous wind speed to approximate the extreme wind speed in Equation (6).The max instantaneous wind can also be estimated by Equation ( 9) by replacing the mean wind speed with the max 2-min wind speed.The hourly max instantaneous to hourly max 2-min wind speed ratio ranged from 1.4 to 1.6 in most stations (Figure 3b).
In this study, a σ was set to be 2.5, and the gust factor was set as 1.25 for the max wind speed in Equation (8) and 1.5 × 1.25 for max wind speed in Equation (6).In order to evaluate the gusty-wind model scheme, the frequency of measured 2-min wind speed at an hourly average wind speed of 10 m/s was used to verify the model.The gusty wind could reproduce the PDF of 2-min wind speed (Figure 4).As shown in Table 1, the WRF model performance on the wind was satisfied and comparable with results from other studies [30].The model slightly overestimated the wind speeds in dust-source regions.Generally, jointing of gusty-wind model and WRF model could reasonably estimate the PDF of wind speed.

Comparison of the Simulated PM 10 Concentrations Generated by the Three Schemes
Three abovementioned vertical dust-emission parameterizations were used to simulate the SDS episode.Figure 5 showed the distributions of total dust-emission flux and the peak PM 10 concentrations during the SDS episode simulated by the three schemes.The MB scheme generated the lowest dust emission flux and PM 10 concentration, and the KOK scheme generated the highest.The dust-emission flux and peak PM 10 concentration generated by the LS scheme were about 6-10 times and 3-6 times than by the MB scheme, respectively.The distribution of peak PM 10 concentration, using the LS scheme, was similar to the distribution of visibility (Figure 2c).The peak PM 10 concentration was higher than 5000 µg/m 3 in the area where the min visibility was lower than 1 km and was higher than 200 µg/m 3 in the area where the min visibility was lower than 7.5 km.The KOK scheme generated much higher PM 10 concentrations, especially in Qinghai, where the min visibility was about 4 km.The MB scheme generated much lower PM 10 concentrations, especially in the Taklimakan desert, where the min visibility was lower than 1 km.(b,d,f)) during the sand and dust episode using three dust-emission schemes ((a,b), using MB scheme; (c,d), using LS scheme; (e,f), using KOK scheme).PM 10 concentrations in six sites were used to evaluate the performance of the three schemes (Figure 6 and Table 2).Three sites (Hohhot, Yinchuan, and Lanzhou) were located near the dust-source areas.The peak PM 10 concentrations exceeded 6500 µg/m 3 in Hohhot and Yinchuan.The LS and KOK schemes reproduced high PM 10 concentrations, but the MB scheme underestimated PM 10 concentrations.The LS scheme produced the minimum error for the peak PM 10 concentrations in Lanzhou, with about −609 µg/m 3 .PM 10 concen-trations in Beijing, Zhengzhou, and Jinan were used to evaluate dust transportation.The MB scheme underestimated the peak PM 10 concentrations at all three sites by large margins, about 70-90%.The LS scheme produced the minimum errors for peak PM 10 concentrations, and the KOK scheme overestimated the peak PM 10 concentrations at all three sites.The IOA was around 0.9 for the three schemes in four sites, indicating the model could capture the variation of PM 10 concentration.The AOD values reached 4.5 over Northern China (Figure 7a) at 10:30 a.m. on 15 March 2021, and the southern boundary of the area with a high AOD value (>2.5) in Southern Hebei and Central Shanxi.All three schemes could reproduce the distribution of AOD over Northern China.Due to the lack of anthropogenic emissions, the model did not generate the distribution of AOD in Southern Asia.The AOD values simulated by the three schemes (Figure 7b-d) were as varied as the simulated PM 10 concentrations.MB scheme underestimated the AOD by a large margin, with the maximum value being only around 1. The KOK scheme overestimated the AOD, with the highest value exceeding 6.The LS scheme performed better on simulating AOD than the other two schemes.However, it still underestimated the AOD in the area around Hetao plain and overestimated the AOD in Northwest China.

The Role of Gusty Wind
The gusty-wind model (Section 2.4) could reproduce the PDF of 2-min average wind, which could be considered as gusty disturbances [41].We used the PDF of 2-min average wind to calculate the PDF friction velocity in one hour and then generated the dustemission flux.
Figure 8 shows the percentage changes of simulated dust-emission flux and max PM 10 concentration caused by implementing the gusty-wind model.In most parts of the dust-source area, the gusty wind generated 5-40% more dust-emission flux than the hourly averaged wind.In the edge area of the dust source, the value reached 200%.These results were consistent with Wang and Zheng [20], where a one-dimensional numerical model of unsteady sand saltation was used to quantify the effects of the unsteady wind on sand-transport rate.Implementing the gusty-wind model could lead to a 5-40% increase in the peak PM 10 concentration in the dust-source area and a 10-50% increase in East-Central China.Therefore, the gusty wind could generate greater dust-emission flux, and more dust particles could be transported to downwind regions.The increase in PM 10 concentration was greater and more widespread than the increase in dust-emission flux.The gusty model generated the biggest percentage increase in dust-emission flux and the peak PM 10 concentrations with the KOK scheme and the smallest with the MB scheme because the vertical flux parameterization was not related to wind speed in the MB scheme.In Section 3.3, the LS scheme performed better when simulating PM 10 concentrations than the other two schemes for this SDS episode.Implementing the gusty wind in the LS scheme could improve the model's performance in regard to simulating the PM 10 concentrations (Table 3).In Section 3.3, the LS scheme underestimated the peak and the mean PM 10 concentrations.The model with the gusty-wind scheme reduced the modeling biases for the peak and the mean PM 10 concentrations.For instance, the error was reduced from −6.7% to −1.4% for mean PM 10 concentration in Beijing.Northern China experienced another SDS event on 27-29 March.This SDS event was much weaker than the SDS event on 14 to 15 March, with lower peak PM 10 concentrations [25,27].In order to investigate the reliability of the gusty-wind model with the LS scheme, the same model configurations were used to simulate the SDS event on 27-29 March.The evaluation results (Table 4) indicated that the LS scheme generated the lowest biases in most cities.However, the biases were much larger for the SDS event on 27-29 March than the SDS event on 14/15 March.The gusty-wind model with the LS scheme overestimated the peak PM 10 concentration by 10~154%.More work is needed to improve model's reliability on various SDS events.

Conclusions
Northern China experienced an exceptional SDS event triggered by a super-strong Mongolian cyclone and loose land surfaces on 14/15 March 2021.There was an atmospheric river that could enhance the intensity of the Mongolia cyclone through the latent heat release.This study compared the performance of three vertical dust-flux parameterizations on the severe sand and dust storm event during 14/15 March 2021 by implementing these three schemes in the CAMx model.There was a tenfold difference between the dustemission flux generated by these three schemes, while the difference between the peak PM 10 concentrations was up to six times.The MB scheme [35] underestimated the PM 10 concentrations by 70-90%, while the KOK scheme [11] overestimated PM 10 concentrations in most areas by 10-50%, especially for Qinghai.The LS scheme [5] produced the minimum errors for peak PM 10 concentrations.
A statistical gusty-wind model was implemented in the dust-emission model.The gusty wind model could reasonably reproduce the probability density function of 2-min wind speeds.In the dust-source regions, compared with the hourly wind, dust-emission flux and peak PM 10 concentrations generated by the gusty wind increased by 5-40%.The percent increase of dust-emission flux due to gusty wind could reach 200% on the edge of the dust-source area.The increase of peak PM 10 concentration caused by gusty wind in the non-dust-source region was higher than that in the dust-source region, with 10-50%.Considering the underestimations of the LS scheme, implementing the gusty-wind model could improve the LS scheme's performance on PM 10 concentration simulation for the severe SDS event on 15 March 2021.
Our study was based on an extreme SDS event, and it was difficult to conclude that a particular scheme was adaptable for modeling and forecasting SDSs in Northern China.For instance, our model did not simulate the SDS event on 27 to 29 March as accurately as it did for the SDS event on 14 to 15 March.More studies for various SDS events are needed to improve the model performance.The gusty-wind model only considers the variation of wind speed.Therefore, further works are expected in developing the gusty-wind model to consider the effects of the coherent structure of gusty wind.

Figure 1 .
Figure 1.Model domain and the PM 10 monitoring stations (the black dots) used to evaluate the model.

Figure 2 .
Figure 2. Atmospheric circulations (a) 500 hPa geopotential height and sea level pressure, (b) 500 hPa geopotential height and integrated water vapor (mm)) at 8:00 a.m. on 15 March 2021 Beijing Time (BJT), and the minimum visibility (color dots) and wind fields at the onset of minimum visibility during the SDS event (c).

Figure 3 .
Figure 3. Distribution of max instantaneous wind speed (a) and statistical factors of gusty wind during the period from 13 to 15 March 2021.(b) Ratio of hourly max instantaneous to hourly max 2-min wind speed, (c) the normalized gust wind speed (Equation (8)) of 2-min wind speed, (d) the correlation coefficient of the standard deviation (σ V ) with the difference between the max and the mean 2-min wind speed in 1-h period, (e) the gust factor (Equation (9)) of 2-min wind speed, and (f) the correlation coefficient of the max with the mean 2-min wind speed in 1-h period.

Figure 4 .
Figure 4. Probability density function (PDF) of measured 2-min wind speed (solid line) and PDF generated by the gusty-wind model (dash line) at an hourly average wind speed of 10 m/s.

Figure 6 .
Figure 6.Comparison of simulated PM 10 concentrations with observations (OBS: the observed PM 10 concentrations, the concertation over 6500 µg/m 3 was not recorded because of the limitation of measurements in Hohhot and Yinchuan; MB, LS, and KOK were the simulated PM 10 concentrations, using MB, LS, and KOK schemes, respectively.Time is BJT).

Figure 8 .
Figure 8. Percentage changes of simulated dust-emission flux (a,c,e) and peak PM 10 concentration (b,d,f) caused by the implementation of gusty-wind model ((a,b) use MB scheme, (c,d) use LS scheme, and (e,f) use KOK scheme).

Table 1 .
Statistic parameters of WRF model performance on wind and temperature during 14-15 March 2021 (Northwest China, including Xinjiang, Gansu, Ningxia, Shaanxi, and Inner Mongolia in China; Northern China, including Beijing, Hebei, Shanxi, Shandong, and Henan in China).

Table 2 .
Statistic parameters of model performance on PM 10 concentrations during SDS event.(TheIOA for Hohhot was not provided, because there were 12 missing values which exceeded 6500 µg/m 3 ).

Table 3 .
Simulation biases of peak PM 10 and mean PM 10 , using LS scheme (LS) and implementing gusty-wind model (LS-GUST).

Table 4 .
Simulated and observed peak PM 10 concentrations during SDS event on 27/28 March 2021 (values in brackets are relative errors).