Spatial-Temporal Variation Characteristics of Vertical Dust Flux Simulated by WRF-Chem Model with GOCART and AFWA Dust Emission Schemes (Case Study: Central Plateau of Iran)

Dust and sand storms are among the major threats to central Iran. These phenomena pose irreparable risks to natural ecosystems and human societies, including effects on health. In this study, the spatial and temporal pattern of vertical dust flux (VDF) was used to identify dust sources as well as areas with high potential for dust generation. To simulate VDF, two intense dust storms, from 21 February 2015 and 14 February 2018, were selected using synoptic data and Moderate Resolution Imaging Spectroradiometer (MODIS) images. These dust storms were identified as responsible for a reduction of horizontal visibility to less than 1000 m, using remote sensing tools and Ackerman Dust Index. MODIS images show that these two storms covered most of Central Plateau of Iran. The Weather Research and Forecasting model with chemistry (WRF-Chem) was used to simulate the storms, with either the Global Ozone Chemistry Aerosol Radiation and Transport (GOCART) or Air Force Weather Agency (AFWA) scheme to calculate VDF. Modeled vertical dust fluxes in both events indicate that the Arabian deserts in Saudi Arabia and in southwestern Iran can be identified as main sources of the dust in the central Iranian plateau. The other source of dust is the Hirmand Basin, located in the country of Afghanistan and in the southeast of Iran. The results of VDF simulations indicate that central southeast Iran could be the main dust source of internal origin. Additionally, over seasonal wetlands in Iran, the amount of VDF was simulated to be sometimes over 4000 μg/(m2s), an indication that these areas are sensitive to wind erosion in dry conditions and can be a source of dust. The WRF-Chem results were compared with the horizontal visibility measured in synoptic stations in the area. The results showed that the coefficients of determination of GOCART results with the measured horizontal visibility on 21 February 2015 and 14 February 2018 were 0.72 and 0.76, respectively, while the coefficient values from the simulations with AFWA scheme on 21 February 2015 and 14 February 2018 with the measured horizontal visibility were lower, 0.44 and 0.50, respectively. Modern-Era Retrospective analysis for Research and Applications, version 2 (MERRA2) re-analysis data also showed timing of peak dust levels consistent with the GOCART scheme. Appl. Sci. 2020, 10, 4536; doi:10.3390/app10134536 www.mdpi.com/journal/applsci Appl. Sci. 2020, 10, 4536 2 of 30


Introduction
Dust particles are among the air pollutants that have a significant impact on the radiation budget [1]. They can also reduce air quality and threaten human health [2], and drastically limit horizontal visibility [3,4]. Dust is generated as a result of wind erosion in arid and semi-arid regions. Since Iran is located in the arid belt of the world (approximately 88% of Iran is arid or semi-arid [5]), wind erosion is one of the most important processes of land degradation and landscape change there [6][7][8]. Dust emission occurs in widespread events that are variable in size, time, and particle concentration [9].
Iran is a country with an average annual rainfall of about 243 mm. Most of the country has a dry or semi-arid climate, which has created dry deserts. The deserts of Loot and Kavir are the most important of them. Low rainfall, high evaporation and, as a result, limited soil moisture and poor vegetation are prominent features of these areas. These conditions make the area prone to wind erosion. Wind erosion involves the separation, transfer, and deposition of soil particles by erosive winds [10,11]. Soil drought, lack of vegetation and erosive winds are the main causes of wind erosion, which is more likely to occur in dry and sparsely vegetated areas. Therefore, the dust phenomenon is one of the most important natural hazards of the Central Plateau of Iran. In case of suitable meteorological and hydrological conditions, conditions for dust storms will be promoted.
Numerical prediction models have been designed to simulate dust emissions. Dust storm prediction models have been used operationally in Europe [12], Australia [13], and East Asia [14][15][16][17]. At the U.S. Naval Research Laboratory, the U.S. Navy Mesoscale Prediction System's Coupled Ocean-Atmospheric Mesoscale Prediction System (COAMPS), with an embedded dust aerosol model [18], has produced three-day forecasts for southwest Asia in real-time runs since March 2003. Due to the climatic and environmental importance of soil dust, considerable efforts have been devoted to simulating the production and transport of soil dust aerosols at regional and global scales [18][19][20][21]. Shao [22] emphasizes a mechanistic model of forces acting on dust particles to predict dust flux. Additionally, Liu et al. [18] used satellite data and horizontal visibility to validate the results of the numerical model of dust prediction. One such model that is able to simulate dust particle transport and its feedbacks simultaneously with meteorological fields is the Weather Research and Forecasting model coupled to Chemistry chemistry (WRF-Chem) model [23,24]. Extensive studies have been carried out using WRF-Chem. Nabavi et al. [25] examined the sensitivity of WRF-Chem to dust emissions in order to determine dust sensitive areas in West Asia. The results showed that the tested implementation of WRF-Chem was skillful in some regions (Iraq and Syria). Su and Fang [26] examined the response of dust transport in WRF-Chem to dust flux as well as land surface characteristics in predicting the sources of dust in East Asia. Xiaolan and Hongsheng [27] investigated the effects of soil moisture on sand leakage and dust emission observed in Hur Qian sandy lands in China, and Rizza et al. [28] investigated several dust events in the central Mediterranean region of Apulia (southeastern Italy) in order to assess the sensitivity of WRF-Chem to land surface models. Song et al. [29] used the WRF-Chem model to simulate dust emission and evaluate its relationship with climate change in East Asia. Teixeira et al. [30] used WRF-Chem to investigate its sensitivity to vertical resolution during a 22-23 June 2012 dust event over North Africa; results indicated that the model was able to transport particles from the ground to higher layers. Tang et al. [31] examined the spatial and temporal characteristics of dust events and their contribution to the budget of aerosols in Eastern Asia, using WRF-Chem and a parametric tool for dust analysis. The results showed that dust loadings peaked in the afternoon and in summer. Chen et al. [32] simulated dust concentration in East Asia during summer 2010 using WRF-Chem.
Additionally, Alizadeh-Choobari et al. [33] used WRF-Chem to show that the direct radiative effects of airborne mineral dust reduce the surface temperature, but warm upper air layers. Additionally, Alizadeh-Choobari et al. [34] used WRF-Chem to model the 120-day wind and dust storm activity pattern over Iran's Sistan Basin.
Several schemes are available to simulate dust emissions depending on different input parameters. Each of these schemes parameterizes wind erosion in different ways. One of the wind erosion schemes that are widely used to predict dust particles is the Georgia Institute of Technology-Goddard Global Ozone Chemistry Aerosol Radiation and Transport (GOCART) model [19,35]. Another widely used wind erosion scheme is the Air Force Weather Agency (AFWA) model [25,[36][37][38].
WRF-Chem is an extended version of the WRF model, which includes a chemistry component [23] and has been widely used to study dust impacts on local to regional scales [39][40][41][42][43][44][45] studies have demonstrated the suitability of the WRF-Chem model to simulate the desert dust cycle [25,32]. However, attention must be paid to the chosen land surface model [28] model vertical resolution [30] and particulate matter distribution [46].
The global erodibility map developed by Ginoux et al. [19] is the only one available by default for dust simulations in WRF/Chem. Such a map works well on large-scale regions (i.e., Sahara Desert, Arabian peninsula) [34]. Zhao et al. [43] investigated the sensitivity of the WRF-Chem simulated dust radiative forcing to dust emission and aerosol schemes. Additionally, modeling of the effect of mineral dust on the wind profile and near-surface wind speed was carried out WRF-Chem model [33].
Additionally, Kargar et al. [47] used the WRF-Chem model to numerically simulate extreme dust storms in eastern Iran. Farhadipoue et al. [48] studied the severe dust storms in west and southwest Iran by the WRF-Chem model. Model outputs showed the good performance of WRF-Chem in simulating the spatial distribution of dust.
In the present study, the most severe dust storms in Iran, which reduced horizontal visibility to less than 1000 m, were selected using data from synoptic stations of the region. Subsequently, these dust storms were evaluated using Moderate Resolution Imaging Spectroradiometer (MODIS) images and the Ackerman index. Then, using the WRF-Chem model and the GOCART and AFWA wind erosion schemes, the spatial and temporal pattern of vertical dust flux and of surface dust concentration were studied and compared with time series of visibility at the surface. The goal of this work is to accurately simulate such severe dust storms in order to determine the external and internal dust sources relevant for these events in Iran. The paper is organized as follows. Materials and methods are presented in Section 2. Results are given in Section 3, while Discussion and Conclusions are drawn in Section 4.

Study Area
The Central Plateau is Iran's largest basin, accounting for most of the country's area and consisting of nine sub-basins, including Loot desert, Central Desert, Jazmoorian, Govkhooni, Salt Lake, Siahkoh Desert, Daranjir Desert, Abargho Desert and Tashk Bakhtegan Maharloo. Rainfall in this area is very low, with an annual average less than 100 mm, and locally as low as 25 mm. On the other hand, the potential evaporation rate in this area is high and in many places it reaches more than 4000 mm/year. As such, the potential evaporation can be as high as 40 to 80 times the annual precipitation [49]. The average relative humidity of the Central Plateau plains is around 30% to 40%, and declines to 15% during the warm season. The mean annual temperature varies between 15 and 30 • C, and the reported maximum and minimum temperatures are 51 and −18 • C, respectively [50].
One of the most common natural events in these areas is dust storms, which can be of external or internal origin. Dust storms of foreign origin originate mostly in the Arabian deserts (Iraq, Syria, and Saudi Arabia), the Qareqoom desert in Turkmenistan. The Hirmand Basin in eastern Iran is an internal dust source, and Sistan's 120-day storms could also activate the deserts of Loot and Kavir as sources of internal dust. Most of the dust storms in this region occur in the spring, and the months of June, July, and May show the highest number of dust days [51,52].
In Figure 1, synoptic stations with data on horizontal visibility during dust storms are shown, as well as the geographical location of the main basins.
Appl. Sci. 2020, 10, x FOR PEER REVIEW 4 of 30 internal dust source, and Sistan's 120-day storms could also activate the deserts of Loot and Kavir as sources of internal dust. Most of the dust storms in this region occur in the spring, and the months of June, July, and May show the highest number of dust days [51,52].
In Figure 1, synoptic stations with data on horizontal visibility during dust storms are shown, as well as the geographical location of the main basins.

Selecting Dust Events
Meteorological codes from 36 synoptic stations were used to select the extreme dust events. Meteorological codes were examined from 2006 to 2018, and the storms with the widest (Dustassociated codes recorded at the most stations), most severe (Severe decrease in horizontal visibility at the synoptic station (below 1,000 or below 500 m)), and longest stay (Number of hours or days of dust observed at the stations based on meteorological codes) in the study area were selected. The information on the meteorological codes of synoptic stations is also shown in Table 1. Whirlwind 08 There is dust when observation in around of the station or already existed at the station. Storm 09 The intensity of the storm has declined over the past hour. Dust storm with mild sand 30 The intensity of the storm has not changed over the past hour.
Dust storm with mild sand 31 The intensity of the storm has increased over the past hour.
Dust storm with mild sand 32 The intensity of the storm has declined over the past hour.
Dust storm with severe sand 33 The intensity of the storm has not changed over the past hour.
Dust storm with severe sand 34 The intensity of the storm has increased over the past hour.
Dust storm with severe sand 35

Selecting Dust Events
Meteorological codes from 36 synoptic stations were used to select the extreme dust events. Meteorological codes were examined from 2006 to 2018, and the storms with the widest (Dust-associated codes recorded at the most stations), most severe (Severe decrease in horizontal visibility at the synoptic station (below 1000 or below 500 m)), and longest stay (Number of hours or days of dust observed at the stations based on meteorological codes) in the study area were selected. The information on the meteorological codes of synoptic stations is also shown in Table 1. Table 1. Dust-related weather codes in synoptic station reports.

Name Weather Codes
Suspended dust that has traveled long distances. Suspended dust 06 Dust or sand rising from the station to the Atmosphere.
Dust 07 Evolving Whirlwind at or near the station.
Whirlwind 08 There is dust when observation in around of the station or already existed at the station. Storm 09 The intensity of the storm has declined over the past hour. Dust storm with mild sand 30 The intensity of the storm has not changed over the past hour.
Dust storm with mild sand 31 The intensity of the storm has increased over the past hour.
Dust storm with mild sand 32 The intensity of the storm has declined over the past hour.
Dust storm with severe sand 33 The intensity of the storm has not changed over the past hour.
Dust storm with severe sand 34 The intensity of the storm has increased over the past hour.
Dust storm with severe sand 35 In order to detect and identify dust events, one may use remote sensing technology and synoptic station data. In this research, first, severe dust storms with horizontal visibility less than 1000 m were selected using synoptic data. In this study, 36 synoptic stations located in the study area were used. Table 2 shows information about the synoptic stations. Next, MODIS satellite images and tools for analyzing the remote sensing data were used to map these events. MOD021km Terra Sensor MODIS (36 bands) data were downloaded from website [53]. MODIS Characterization in this study including: Spatial resolution: 1 km, temporal resolution: 5 min, MODIS collection 6.1, possessing level 1B.
For the detection of dust storms, MODIS images were downloaded at desired dates and transformed into HDVI format with ENVI5.1 software (5.1 version, L3harris Geospatial Company, Broomfield, CO, USA, 2013). Geometrical corrections and brightness temperatures of the bands were calculated. After receiving and pre-processing MODIS satellite data, the Brightness Temperature Difference (BTD) method [54] was used to detect the presence of dust. This method considers the brightness temperature difference between band 31 (11 µm) and 32 (12 µm) to detect dust storms; the value of zero is considered as the threshold for dust detection (where values less than 0 K indicate dust).
In this study, because the Aerosol Optical Depth (AOD) of the region on these dates had many missing values due to the presence of clouds, the Ackerman index was used instead.

Determination of Dust Storm Origin Using HYSPLIT
The Hybrid Single Particle Lagrangian Integrated Trajectory (HYSPLIT) model was designed by [55] and has been upgraded over the past few decades. This model is in fact a dual model for calculating dust motion, scattering, and particle simulation. This model is used to simulate dust trajectories and identify areas of origin as well as areas affected by dust storms. HYSPLIT is one of the most widely used models for atmospheric trajectory and dispersion calculations [56,57]. In this study, HYSPLIT calculations were performed online for the identified two events for 24 h integration periods.

WRF-Chem Model
The third version of the Weather Research and Forecasting model coupled to Chemistry (WRF-Chem V3, National Center for Atmospheric research, Boulder, CO, USA, 2011; https://ruc.noaa.gov/wrf/wrfchem/) was used here. WRF-Chem is a non-hydrostatic limited-area numerical model, able to simulate both meteorological fields and dust particle transport and their reciprocal feedbacks. It is both a research and operational model, which simulates and predicts the atmospheric concentrations of different gases and particles, taking into consideration natural mechanisms and human activities (pollutants) that affect dust particles of different sizes, considering the ways they are emitted and deposited.

Dust Schemes in WRF-Chem
From version 3.7 onwards, the WRF-Chem model provides a choice of dust schemes. The most important dust schemes in the WRF-Chem model are GOCART, AFWA and University of Cologne (UoC) [22,58,59]. In this study, according to the volume of the study and also according to the previous research [60], GOCART and AFWA schemas were considered suitable for this region. As follows:

GOCART Dust Scheme
In the GOCART scheme, the particle source (S) function is defined as the ratio of erodible soil subject to the wind erosion process: where S can be interpreted as the probability of sediment accumulation at point i of height z i ; z max and z min are the maximum and minimum elevations in a range of 10 • × 10 • around the grid point, respectively. Dust particle emission fluxes in the GOCART scheme are defined by the following empirical relationship [61]: where F p is (µg/m 2 s) represents the emission flux for size bin p, C G is a dimensionality coefficient (µg s 2 m −5 ), (C G = 1.2 µg s 2 m −5 ), s p is the fraction of each particle size class, S is the source function representing the fraction of alluvium available for wind erosion [41] and U t is a threshold velocity at 10 m, U 10 is horizontal wind speed at 10 m, and U t is the threshold velocity (m s −1 ) below which dust emission does not occur and is a function of particle size, air density and surface moisture [41]. Ginoux et al. [19] suggested that U t ≈ u * t , where u * t is a threshold friction velocity, using the formula of [62] and considering the effect of soil moisture.

Dust Scheme AFWA
The parameterization method proposed by Marticorena and Bergametti [62] for the dust emission was implemented as a scheme or dust module (AFWA) by the United States (US) Air Force Meteorological Center in the WRF-Chem model. In the AFWA scheme, the correction factor is applied to u * t as follows: where f (roughness) is a drag partition correction-if f (roughness) is 1, it means the surface is smooth; Value decreases with increasing amounts of large rocks, cobbles, vegetation, etc. f (moisture) is the soil moisture content, and D p is the particle diameter in bin p. Thus, u * t increases with increasing rockiness, vegetation, and other surface types with high roughness. The f (roughness) value for Southwest Asia is currently set to 1, because the great deserts of Saudi Arabia, Iraq, Syria, Iran, etc. are in this region. Vegetation in these areas is very limited, the land surface is bare and wind erosion systems are active in these areas. For these reasons, the f (roughness)t for these areas is estimated to be 1 (smooth surface). Dust emission is restricted to areas with roughness length Z 0 ≤ 5 m (typically barren lands and sparsely vegetated surfaces) [63].
In the AFWA scheme, the dust flux H caused by the saltation process on the bare soil surface is obtained from the relation introduced by Kawamura [64] as follows: where C is an empirical constant (0.129), ρ a is the density of air parcel, (g/cm 3 ), g is the acceleration due to gravity (cm/s 2 ), u * and u * t (in units of cm/s) are, respectively, the friction velocity and the threshold friction velocity. The friction velocity is considered instead of 10 m horizontal wind speed as in the GOCART scheme.
In this scheme, the vertical mass flux of dust is calculated based on Marticorena and Bergametti [62] with a correction coefficient.
where α is the saltation ratio which increases proportionally to the clay percent of the soil and S is the Dust Source function. Additionally, this equation is only valid for 0% < %(clay) < 20%.

WRF-Chem Model Configuration
In this study, WRF-Chem version 3.9.1 was used to simulate dust storms. Two domains were defined for the model. The first domain was defined in the area with the number of grid points being 90 × 100 and the spacing of 11.2 km. The second domain was defined in the area showed in Figure 2. Initial and boundary conditions were obtained from the National Oceanic and Atmospheric Administration Global Forecast System (NOAA GFS) with horizontal resolution of 0.5 • × 0.5 • United States Geological Survey (USGS) data were used for static geographic data such as altitude, soil properties, vegetation fraction and land use. Table 3 lists the parameterization schemes used to implement the WRF-Chem model. Additionally, the distribution of dust particle diameters in WRF-Chem model is shown in Table 4. To simulate the investigated storms, the model was run starting 24 h before the storm, and the first 3 h of each run were considered as spin-up time. The time interval between networks was 180 s. GFS data with a 6-h time interval was used for initial and boundary conditions. The model output was also saved every 3 h.

Verification of Results
In this section, the results of GOCART and AFWA wind erosion schemas were validated by comparing observed horizontal visibility with simulated surface concentration. One of the most common methods of validation is the use of the Aerosol Robotic Network (AERONET) [65] and Aerosol Optical Depth (AOD) data. Remote sensing AOD data is not available on the mentioned dates because of the presence of clouds in the area. Additionally, there were no AERONET data in these two dates in the region. Therefore, station observations of horizontal visibility must be used. We also used the Modern-Era Retrospective Analysis for Research and Applications, version 2 (MERRA-2) re-analysis, comparing the surface dust concentration extracted from the schemas with that in MERRA-2 re-analysis data. MERRA-2 is the first long-term global reanalysis to assimilate space-based observations of aerosols and represent their interactions with other physical processes in the climate system [66].

Dust Storm Selectin
For this study, using the values of horizontal visibility, two severe dust storms that occurred in To simulate the investigated storms, the model was run starting 24 h before the storm, and the first 3 h of each run were considered as spin-up time. The time interval between networks was 180 s. GFS data with a 6-h time interval was used for initial and boundary conditions. The model output was also saved every 3 h.

Verification of Results
In this section, the results of GOCART and AFWA wind erosion schemas were validated by comparing observed horizontal visibility with simulated surface concentration. One of the most common methods of validation is the use of the Aerosol Robotic Network (AERONET) [65] and Aerosol Optical Depth (AOD) data. Remote sensing AOD data is not available on the mentioned dates because of the presence of clouds in the area. Additionally, there were no AERONET data in these two dates in the region. Therefore, station observations of horizontal visibility must be used. We also used the Modern-Era Retrospective Analysis for Research and Applications, version 2 (MERRA-2) re-analysis, comparing the surface dust concentration extracted from the schemas with that in MERRA-2 re-analysis data. MERRA-2 is the first long-term global reanalysis to assimilate space-based observations of aerosols and represent their interactions with other physical processes in the climate system [66].

Dust Storm Selectin
For this study, using the values of horizontal visibility, two severe dust storms that occurred in winter were selected. The two cases of 21 February 2015 and 14 February 2018 were detected and simulated. Figure 3 represents the horizontal visibility of some synoptic stations in the study area during the two storms. These two events occurred in very dry winter seasons, during which the rainfall in much of the study area (Central Iranian Plateau) through the day of the event was lower than normal.

Detection with MODIS Images
The dust storms of 21 February 2015 and 14 February 2018 were analyzed using the Ackerman index on the MODIS images. The results of the analysis of MODIS images of dust storms are shown in Figure 4.     Back trajectory analysis supports the origin of air reaching Central Iran during these days ( Figure 5). The results show that Arabian deserts affect Iran's atmosphere during dust storms. Figure 5 shows that the storm on 21 February 2015 originated in the Arabian Desert, which can be identified as a storm of outside origin. Additionally, the storm of 14 February 2018, which was observed in the eastern stations of the study area, originated from the southwest of Iran (probably the Arabian deserts). Back trajectory analysis supports the origin of air reaching Central Iran during these days ( Figure 5). The results show that Arabian deserts affect Iran's atmosphere during dust storms. Figure  5 shows that the storm on 21 February 2015 originated in the Arabian Desert, which can be identified as a storm of outside origin. Additionally, the storm of 14 February 2018, which was observed in the eastern stations of the study area, originated from the southwest of Iran (probably the Arabian deserts).

Simulation of Dust Storm on 14 February 2018
Vertical dust flux represents the potential for soil particles to rise and detach; in other words, places with large vertical dust flux can be identified as dust sources. Figure 6 shows the three-hour vertical dust flux for 14 February 2018 storm simulated by the AFWA scheme. One may see that the model is able to reproduce the distribution of dust storm in MODIS image at 09:00 UTC; in particular, the peak of concentration in central Iran and in its southeastern part are well detected.

Spatial Pattern of Vertical Dust Flux Using AFWA Wind Erosion Scheme
The analysis of the vertical dust flux maps on 14 February 2018 using AFWA Wind Erosion Model in the Central Plateau of Iran shows a dust source in the south-western part of the Arabian Desert mainly from 6:00 to 9:00 UTC on 14 February 2018, and another source in Central Iran. This is shown more clearly in Figure 6, where the average spatial pattern of dust flux during the event is shown. As one may see, the dust flux reaches as high as 2,838 µ g/(m 2 s). The most dust rise is in parts of the Gavkhouni wetland and Salt Lake, as shown in Figure 7. During winter, which is generally a wet season, these areas should have a wet surface; however, the results show that if the soil moisture

Simulation of Dust Storm on 14 February 2018
Vertical dust flux represents the potential for soil particles to rise and detach; in other words, places with large vertical dust flux can be identified as dust sources. Figure 6 shows the three-hour vertical dust flux for 14 February 2018 storm simulated by the AFWA scheme. One may see that the model is able to reproduce the distribution of dust storm in MODIS image at 09:00 UTC; in particular, the peak of concentration in central Iran and in its south-eastern part are well detected.

Spatial Pattern of Vertical Dust Flux Using AFWA Wind Erosion Scheme
is reduced, as was the case in 2018, the region may serve as a source of dust as well, potentially becoming a dust hotspot.   The analysis of the vertical dust flux maps on 14 February 2018 using AFWA Wind Erosion Model in the Central Plateau of Iran shows a dust source in the south-western part of the Arabian Desert mainly from 6:00 to 9:00 UTC on 14 February 2018, and another source in Central Iran. This is shown more clearly in Figure 6, where the average spatial pattern of dust flux during the event is shown. As one may see, the dust flux reaches as high as 2838 µg/(m 2 s). The most dust rise is in parts of the Gavkhouni wetland and Salt Lake, as shown in Figure 7. During winter, which is generally a wet season, these areas should have a wet surface; however, the results show that if the soil moisture is reduced, as was the case in 2018, the region may serve as a source of dust as well, potentially becoming a dust hotspot.

Spatial Pattern of Vertical Dust Flux Using GOCART Wind Erosion Model
The GOCART scheme calculated the spatial pattern of vertical dust flux over a 3-h interval using the GFS analysis/forecast as large scale forcing on 14 February 2018, as shown in Figure 8. Appl  The average spatial pattern of vertical dust flux (VDF) simulated with the GOCART model indicates that the Jazmourian Wetland Basin can be a source of dust in the southeast of Iran (Figure 9). Additionally, it indicates that parts of the Loot Desert in Kerman, Sistan, and Baluchestan provinces can be internal dust sources. In fact, the simulated vertical flux of dust by the GOCART scheme reaches 7536 µg/(m 2 s) in the area. The average spatial pattern of vertical dust flux (VDF) simulated with the GOCART model indicates that the Jazmourian Wetland Basin can be a source of dust in the southeast of Iran ( Figure  9). Additionally, it indicates that parts of the Loot Desert in Kerman, Sistan, and Baluchestan provinces can be internal dust sources. In fact, the simulated vertical flux of dust by the GOCART scheme reaches 7,536 µ g/(m 2 s) in the area. Overall, the VDF results using GOCART suggest that central, eastern and south-eastern parts of Iran have the highest potential for dust emission and these areas may be considered as internal Sources.
According to the results of two dust schemes, the vertical flux values simulated by AFWA are lower than those simulated by the GOCART scheme. The number of dust sources shown by the GOCART scheme in the study area is also higher than for the AFWA scheme. Therefore, as a final step, in order to verify and compare the schemes, the correlation between the results of the surface dust concentration simulated by the AFWA and GOCART schemes with the horizontal view of the synoptic stations will be investigated.

Spatial Pattern of Vertical Dust Flux Using AFWA Wind Erosion Scheme
The results of the simulated vertical dust flux of 21 February 2015 are shown in Figure 10. Results show that the main dust source of this storm is the Arabian Desert; focusing on the internal sources, it appears that the Loot Desert, Jazmoorian and Daranjir desert have the highest potential to generate dust in the central Iranian plateau (Figure 11). Overall, the VDF results using GOCART suggest that central, eastern and south-eastern parts of Iran have the highest potential for dust emission and these areas may be considered as internal Sources.
According to the results of two dust schemes, the vertical flux values simulated by AFWA are lower than those simulated by the GOCART scheme. The number of dust sources shown by the GOCART scheme in the study area is also higher than for the AFWA scheme. Therefore, as a final step, in order to verify and compare the schemes, the correlation between the results of the surface dust concentration simulated by the AFWA and GOCART schemes with the horizontal view of the synoptic stations will be investigated.

Spatial Pattern of Vertical Dust Flux Using AFWA Wind Erosion Scheme
The results of the simulated vertical dust flux of 21 February 2015 are shown in Figure 10.
Hence, the AFWA Scheme identified 4 sources of dust in the Loot Desert, Jazmourian, Hirmand and Daranjir Desert. The VDF simulated by this scheme in the Loot Desert reaches 920 μg/(m 2 s) ( Figure 11  Results show that the main dust source of this storm is the Arabian Desert; focusing on the internal sources, it appears that the Loot Desert, Jazmoorian and Daranjir desert have the highest potential to generate dust in the central Iranian plateau (Figure 11).

Spatial Pattern of Vertical Dust Flux Using GOCART Wind Erosion Scheme
The GOCART scheme simulated the vertical flux (μg/m 2 s) of the dust storm on 21 February 2015, with results as shown in Figure 12.
Storm dust emission flux results on 21 February 2015 by GOCART scheme using WRF-Chem model show that Saudi Arabian deserts (to the south-west of the study area), Qareqoom desert (to the north-east) and Hirmand Basin (to the south-east) formed the primary core of the storm; additionally, emission flux simulations showed that the south-eastern regions of Iran as well as northeast Iran have the potential for dust emission. There is high VDF in the Loot Desert and in the Central Plateau; hence, the elevated ground levels in these areas and favorable climatic conditions can cause severe dust emission. The average spatial pattern in Figure 13 also shows that the Jazmourian basin in southeastern Iran may be identified as a critical hotspot. The GOCART scheme illustrates this phenomenon well. Hence, the AFWA Scheme identified 4 sources of dust in the Loot Desert, Jazmourian, Hirmand and Daranjir Desert. The VDF simulated by this scheme in the Loot Desert reaches 960 µg/(m 2 s) ( Figure 11).

Spatial Pattern of Vertical Dust Flux Using GOCART Wind Erosion Scheme
The GOCART scheme simulated the vertical flux (µg/m 2 s) of the dust storm on 21 February 2015, with results as shown in Figure 12.
Storm dust emission flux results on 21 February 2015 by GOCART scheme using WRF-Chem model show that Saudi Arabian deserts (to the south-west of the study area), Qareqoom desert (to the north-east) and Hirmand Basin (to the south-east) formed the primary core of the storm; additionally, emission flux simulations showed that the south-eastern regions of Iran as well as northeast Iran have the potential for dust emission. There is high VDF in the Loot Desert and in the Central Plateau; hence, the elevated ground levels in these areas and favorable climatic conditions can cause severe dust emission. The average spatial pattern in Figure 13 also shows that the Jazmourian basin in southeastern Iran may be identified as a critical hotspot. The GOCART scheme illustrates this phenomenon well.   GOCART for 21 February 2015 shows internal dust sources in the Jazmourian and Loot Desert and Hirmand Basin. In the Jazmourian Basin, the VDF simulated by the model is up to 14,000 µg/m 2 s.

3.5.Validation of WRF-Chem with GOCART and AFWA Wind Erosion Schemes
The results of the simulation of the vertical dust flux by the two schemes are significantly different, so it is necessary to evaluate their results using the surface dust concentration (SDC) simulated by the model compared to the horizontal visibility of the synoptic stations and MERRA2 dataset.

Validation of WRF-Chem by Horizontal Visibility
Surface concentration values simulated by GOCART and AFWA schemes as well as the horizontal visibility of synoptic stations were analyzed using the Pearson's correlation. The results showed that an exponential decay function fitted to the surface dust concentration simulated by GOCART scheme had a coefficient of determination of 0.73 with horizontal visibility values in storm 14 February 2018 and of 0.77 in storm 21 February 2015, which were significantly different from zero at the 95% level ( Figure 14).  Figure 15). These results suggest that the values simulated by the GOCART scheme are better than the AFWA scheme. GOCART for 21 February 2015 shows internal dust sources in the Jazmourian and Loot Desert and Hirmand Basin. In the Jazmourian Basin, the VDF simulated by the model is up to 14,000 µg/m 2 s.

Validation of WRF-Chem with GOCART and AFWA Wind Erosion Schemes
The results of the simulation of the vertical dust flux by the two schemes are significantly different, so it is necessary to evaluate their results using the surface dust concentration (SDC) simulated by the model compared to the horizontal visibility of the synoptic stations and MERRA2 dataset.

Validation of WRF-Chem by Horizontal Visibility
Surface concentration values simulated by GOCART and AFWA schemes as well as the horizontal visibility of synoptic stations were analyzed using the Pearson's correlation. The results showed that an exponential decay function fitted to the surface dust concentration simulated by GOCART scheme had a coefficient of determination of 0.73 with horizontal visibility values in storm 14 February 2018 and of 0.77 in storm 21 February 2015, which were significantly different from zero at the 95% level ( Figure 14).   Figure 16.
Similarly, the values of dust surface concentration were also simulated for 14 February 2018 using the GOCART scheme. Results for modeled dust concentration at a height of 1.8 m above the surface and 3-h time intervals are shown in Figure 17.
The time series of the dust surface concentration with the two schemes are shown in Figures 18  and 19. In order to validate the time series of simulated concentrations by the GOCART and AFWA schemes, dust surface concentration from the MERRA-2 reanalysis was used, as shown in Figure 20. The results of comparing the time series diagram of the three-hour time series model simulated by GOCART schema with the time series from MERRA-2 dataset show that the highest amount of dust concentration in MERRA-2 is at 0 and 3 h on February 14, and this concentration is similar to that simulated by the GOCART scheme. However, the time series simulated with the AFWA scheme ( Figure 19) does not match well with that obtained from MERRA-2.   Figure 16.
Similarly, the values of dust surface concentration were also simulated for 14 February 2018 using the GOCART scheme. Results for modeled dust concentration at a height of 1.8 m above the surface and 3-h time intervals are shown in Figure 17.
The time series of the dust surface concentration with the two schemes are shown in Figures 18  and 19. In order to validate the time series of simulated concentrations by the GOCART and AFWA schemes, dust surface concentration from the MERRA-2 reanalysis was used, as shown in Figure 20. The results of comparing the time series diagram of the three-hour time series model simulated by GOCART schema with the time series from MERRA-2 dataset show that the highest amount of dust concentration in MERRA-2 is at 0 and 3 h on February 14, and this concentration is similar to that simulated by the GOCART scheme. However, the time series simulated with the AFWA scheme ( Figure 19) does not match well with that obtained from MERRA-2.      The time series of the dust surface concentration with the two schemes are shown in Figures 18  and 19. In order to validate the time series of simulated concentrations by the GOCART and AFWA schemes, dust surface concentration from the MERRA-2 reanalysis was used, as shown in Figure 20. The results of comparing the time series diagram of the three-hour time series model simulated by GOCART schema with the time series from MERRA-2 dataset show that the highest amount of dust concentration in MERRA-2 is at 0 and 3 h on February 14, and this concentration is similar to that simulated by the GOCART scheme. However, the time series simulated with the AFWA scheme ( Figure 19) does not match well with that obtained from MERRA-2.

Discussion and Conclusions
In the present study, the WRF-Chem model simulated the dust emission flux during the storms of 14 February 2018 and 21 February 2015. The selection criterion for the storms was observed horizontal visibility of less than 1000 m, which was determined using the data of synoptic stations. In order to validate the synoptic data as well as to investigate the extent of the storms in the study area, MODIS images were used to show the dust on 14 February 2018 and 21 February 2015 in Central Iran. MODIS images confirmed the horizontal visibility observations at synoptic stations. Additionally, the NOAA HYSPLIT model showed that airmasses originating in external dust source regions affected much of central Iran. The results of back trajectories from the NOAA HYSPLIT model show that dust storms on the February 2018 and 2015 involved flow from the Arabian Desert to the Central Plateau. In this study, two schemes were used for WRF-Chem simulation. Different domains were employed for the two events to capture the dust source suggested by the back trajectory analysis. Vertical dust flux maps were simulated using the AFWA and GOCART schemes, and the output fields analyzed at 3-h intervals. The results showed that the WRF-Chem model can simulate the dust phenomenon well and is able to identify highly erodible surfaces with both schemes. This study showed that the Arabian deserts play a key role in the formation of the core of Iran's dust storms. Given that Persian Gulf is in between Iran and the Saudi Arabian deserts, one may understand that conditions are also difficult for shipping due to the reduced horizontal visibility when dust storms occur in these areas. The results also show that Loot desert and central desert are two of the main dust sources in Iran.
These severe storms occurred in winter, which is the cool and rainy season of the region. However, due to climate change that is reducing rainfall in the area, wetland areas are becoming drier even during the cool and rainy season and may become a dust hot spot. Due to the significant decrease in precipitation in the study area and the severe dryness of soil surfaces, the VDF from the Loot desert reached over 4000 µg/(m 2 s). The results of the simulation of the storm of 21 February 2015 in the Central Plateau of Iran by the WRF-Chem model show that the central core of this phenomenon is in the Arabian deserts, along with, in south-eastern Iran, parts of Loot desert and Jazmoorian wetland, which are considered as the internal cores of this storm. Horizontal visibility of synoptic stations and surface concentration simulated by GOCART and AFWA schemes were used to validate and compare the results of the two schemes. The results showed that the values simulated by the GOCART model had a significant correlation with the values of the horizontal visibility at the synoptic stations, but the results of the AFWA model did not show as high a correlation. The number of source regions identified by the GOCART model in the storm 14 February 2018 is eight, including two sources in Loot desert, one source in Jazmourian, one source in Hirmand basin, two sources in Central desert and one source in Dranjir desert. Additionally, the number of sources identified by the GOCART model in the storm of 21 February 2015 is four with two sources in Loot desert, one source in Jazmourian and one source in Hirmand basin.
Several previous studies are consistent with the result of this paper. Cao et al. [67], in a study to determine dust sources in Iran, stated that in south-eastern Iran, Sistan basin can be characterized as an internal dust source. Additionally, Alizadeh-Choobari et al. [34] also pointed to the importance of the Sistan basin in eastern Iran and the 120-day winds in this area, which is one of the most important causes of dust in the eastern and south-eastern regions of Iran. This is consistent with the results of this research (identification of dust centres in Hirmand Basin); Baghbanan et al. [68], studying the spatial pattern of dust storms in Iran, showed that the main source of dust is in southern and south-eastern parts of Iran, while Rezazadeh et al. [69] identified the Arabian deserts and the country of Afghanistan as external dust sources into Iran.
The dust phenomena in Iran, especially in the Central Plateau, is very important for its ecological, economic, and health impacts. The identification of dust sources and planning for their management is of particular importance. Additionally, for the cases analyzed in this study, we identified the Arabian deserts in Saudi Arabia as external dust sources. This illustrates that dust is an international phenomenon and requires international interactions to manage it.