Development of a dynamic dust source map for NMME-DREAM v 1 . 0 model based on MODIS Normalized Difference Vegetation Index ( NDVI ) over the Arabian Peninsula

We developed a time-dependent dust source map for the NMME Dust Regional Atmospheric Model (DREAM v1.0) based on the satellite MODIS Normalized Difference Vegetation Index (NDVI). Areas with NDVI < 0.1 are classified as active dust sources. The updated modeling system is tested for dust emission capabilities over SW Asia using a mesoscale model grid increment of 0.1× 0.1 for a period of 1 year (2016). Our results indicate significant deviations in simulated aerosol optical depths (AODs) compared to the static dust source approach and general increase in dust loads over the selected domain. Comparison with MODIS AOD indicates a more realistic spatial distribution of dust in the dynamic source simulations compared to the static dust sources approach. The modeled AOD bias is improved from −0.140 to 0.083 for the case of dust events (i.e., for AOD> 0.25) and from−0.933 to−0.424 for dust episodes with AOD> 1. This new development can be easily applied to other time periods, models, and different areas worldwide for a local fine tuning of the parameterization and assessment of its performance.


Introduction
The importance of natural particles, namely desert dust, in the weather and climate has been underlined in a great number of studies.Dust is a climatic regulator, as it modifies extensively the radiative balance of the atmospheric column (e.g., Torge et al., 2011;Spyrou et al., 2013;Mahowald et al., 2014).At the same time dust aerosols modify the atmospheric water content (Spyrou, 2018), the way clouds are formed by acting as cloud condensation nuclei (CCN) and ice nuclei (IN), and the precipitation process (Kumar et al., 2011;Solomos et al., 2011;Nickovic et al., 2016).In addition, there is a clear connection between dust particles and human health disorders, as the size of the aerosols produced is small enough to cause respiratory and cardiovascular diseases, as well as pathogenic conditions due to the microorganisms that they can potentially carry (Mitsakou et al., 2008;Esmaeil et al., 2014).
The Arabian Peninsula is one of the most important sources of mineral dust worldwide and together with the Sahara and the Gobi contributes to the formation of a Northern Hemisphere "dust belt" as described by Prospero et al. (2002).Severe dust storms over the Peninsula are quite common, especially during long periods without rain, in the spring and summer (Almazrouia et al., 2012).Particles injected into the atmosphere from arid soils, under favorable Published by Copernicus Publications on behalf of the European Geosciences Union.S. Solomos et al.: Development of an NDVI dynamic dust source map for NMME-DREAM weather conditions (high wind speeds and dry soil), can affect large areas around the sources but also remote locations like the Eastern Mediterranean (Mamouri et al., 2016;Solomos et al., 2017) and the Indian Ocean (Chakraborty et al., 2006).
Due to the multitude and severe effects of dust particles not only on the weather and the ecosystem but on human health as well, the proper description of the production, transport, and eventual deposition of the dust cycle in numerical weather prediction models (NWPs) is essential.In order to be able to accurately describe the dust life cycle in the atmosphere, we need a clear understanding of the areas which can potentially act as "dust sources".The definition of such areas dictates the emission strength and therefore the amount of particles inserted into the atmosphere.A proper representation of dust sources is therefore an essential first step, in studying the impacts of mineral particles in the climate and human societies.Usually the definition of the areas that can act as dust sources is made using global datasets.For example Nickovic et al. (2001) used a subjective correspondence between the Olson World Ecosystems (Olson et al., 1983) and the 13 SSib (simplified simple biosphere; Xue et al., 1991) vegetation types to identify arid and semiarid areas.Similarly, Spyrou et al. ( 2010) used a 30 s global land use/cover database, classified according to the 24 category U.S. Geological Survey (USGS) land use/cover system (Anderson et al., 1976) to define active areas in the SKIRON dust model.Solomos et al. (2011) used the LEAF soil and vegetation sub-model of the Regional Atmospheric Modeling System (RAMS) (Walko et al., 2000) to identify the active dust sources in the RAMS-ICLAMS model.
However, the abovementioned methodologies have some significant drawbacks.The datasets are usually not up to date; therefore, recent land use modifications are not included and not represented.In addition, such "static" databases mean that possible seasonal variations are not taken into account.In order to overcome the above limitations and improving global dust forecasts, Kim et al. (2013) developed a dynamical dust source map for the GOCART dust model by characterizing Normalized Difference Vegetation Index (NDVI) values < 0.15 as active dust spots.Similarly Vukovic et al. (2014) combined MODIS land cover types with pixels having NDVI < 0.1 to identify the seasonal dust sources that caused the severe Phoenix haboob of July 2011 in the US.Such information can be even more relevant at meso-and local scales for determining land use changes and potential dust sources, especially in heterogeneous regions such as the Arabian Peninsula (which has more diverse soil types than, e.g., the Sahara) and the greater SW Asia.In this context, Solomos et al. (2017), used the Landsat-8 NDVI data (assuming also NDVI < 0.1 as active sources) to identify recent changes in land use due to the war in Iraq and Syria resulting in a significantly more realistic simulation of dust properties in the Middle East.
In the current study we present the implementation of a dynamical dust source map in the well-established and widely used the Dust Regional Atmospheric Model (DREAM v1.0) (Nickovic et al., 2001;Pérez et al., 2006).The new development is first tested here for the greater SW Asia but can be extended for use in mesoscale dust modeling applications worldwide.Two experimental simulations are performed for a 1-month period (August 2016) over the greater SW Asia: (1) a control run, where the dust source definition is based on the Ginoux et al. (2001) dataset, and (2) a dynamic source run, where the NDVI values are used to identify the dust sources.The main differences in our approach compared to the previous studies referenced above are that we use a very high-resolution NDVI product (500 m × 500 m) in a regional modeling domain (e.g., Kim et al., 2013 used an 8 km × 8 km NDVI dataset extrapolated to 1 • × 1 • global modeling domain) and our study is not limited to specific test cases (like, for example, Vukovic et al., 2014 andSolomos et al., 2017), but covers an extended time period, as presented below.The model results from both runs are compared to available satellite observations and station measurements inside the modeling domain.In Sect. 1 we describe the methodological steps regarding the model developments and remote-sensing data; Sect. 2 includes the results of the experimental runs and Sect. 3 is a summary and discussion of the study findings.

Model description
The modeling system used in this study is NMME-DREAM v1.0.The meteorological core is the NCEP-NMME atmospheric model (Janjic et al., 2001).DREAM v1.0 is a numerical model created with the main purpose of simulating and predicting the atmospheric life cycle of mineral dust using an Euler-type nonlinear partial differential equation for dust mass continuity (Nickovic et al., 2001(Nickovic et al., , 2016;;Pérez et al., 2006;Pejanovic et al., 2011).In DREAM the concentration approach is used for dust uplift, where surface concentration is used as a lower boundary condition and used for the calculation of surface fluxes, which in turn depends on the friction velocity (Nickovic et al., 2001).This surface concentration is calculated using Eq. ( 11) from Nickovic et al. (2001): where c 1 = 2.4 • 10 −4 Kgr m 5 s 2 is a constant determined from model experiments, u * and u * t are the friction velocity and the threshold friction velocity for dust production, respectively, and δ = a •γ k •β k , where γ k the ratio between the mass available for uplift and the total mass β k the fractions of clay, silt, and sand for each soil class and a is the desert mask (between 0 and 1) calculated from the Ginoux et al. (2001) dataset.Soil moisture and particle size dictate the threshold friction velocity which initializes dust production.Once particles have been lifted from the ground they are driven by the atmospheric model variables and processes.Therefore, turbulent parameters are used in the beginning of the process, when dust is lifted from the ground, and transported by model winds in the later phases, when dust travels away from the sources.The model handles dust in eight size bins, with effective radii of 0. 15, 0.25, 0.45, 0.78, 1.3, 2.2, 3.8,  and 7.1 mm.Dust is treated as a passive tracer and does not interact with radiation or clouds.Dust is eventually settled through rainfall and/or dry deposition processes parameterized according to the scheme of Georgi (1986), which includes deposition by surface turbulent and Brownian diffusion, gravitational settling, and impact on surface elements.
In order to test the use of NDVI for source characterization, the model is set up with a horizontal resolution of 0.1 • ×0.1 • , covering the Arabian Peninsula parts of SW Asia and parts of NE Africa (Fig. 1).In the vertical we use 28 levels stretching from the surface to the top of the atmosphere.August 2016 has been selected as a test period for the model development due to the significant dust activity and variability in wind properties during this month.Oneyear runs for all of 2016 have been conducted to evaluate the performance of the static and dynamic database emission maps.The original classification of dust sources in DREAM is based on Ginoux et al. (2001), which takes into account the preferential sources related to topographic depressions and paleolake sediments.The global mapping of dust sources in Ginoux et al. (2001) is determined from the comparison between the elevation of surface grid points at 1 • ×1 • resolution with the surrounding hydrological basins and with the 1 • ×1 • AVHRR (Advanced Very High Resolution Radiometer) vegetation map (DeFries and Townshend, 1994).Recent studies indicated the contribution of both natural and anthropogenic dust sources to the overall dust emissions detected by the MODIS Deep Blue product (Ginoux et al., 2012) and also the relevance of local geomorphological conditions and sediment supply (Parajuli and Zender, 2017) for the global dust emissions.All these advances in dust emissions are based on static map considerations.
In our work, a numerical procedure has been developed to insert the NDVI satellite information into the model and to update such information each time the NDVI changes during the simulation period.We assume that regions with NDVI values from 0 to 0.1 correspond to bare soil and therefore can be efficient sources ("dust points"; DeFries and Townshend, 1994;Solomos et al., 2017).In general it is not easy to define a global threshold value for all satellite NDVI sensors and all vegetation types worldwide.For example Kim et al. (2013) used a threshold of 0.15 to define global dust sources based on AVHRR retrievals (Tucker et al., 2005;Brown et al., 2006).Here we adopt the 0.1 NDVI threshold due to the bareness of the specific modeling domain since a higher value could overestimate the regional dust sources.The NDVI dataset is at a finer resolution than the model grid (500 m × 500 m), and in order to find the potential for dust production in each model grid box, we calculate the following ratio: where #_of_dust_points is the number of points with NDVI values smaller than 0.1.This approach allows for a dynamic description of dust source areas over the model domain to replace the previously used static database.Moreover, the scaling of satellite data over model grid points allows the use of the same algorithm for different model configurations.Several mountains in the area (e.g., the Sarawat Mountains along the Red Sea coast and the Zagros Mountains in Iraq) could be misclassified as dust sources due to low NDVI values.
In order to exclude such unrealistic emissions from non-soil bare areas or snow-covered areas we have applied a limit of zero dust production above 2500 m over the entire domain.This simple approach has been selected in order to keep our straightforward NDVI mapping independent of vegetation and soil information.The threshold value of 2500 m does not suppress the emissions from lowlands and hillsides (e.g., the coastal areas of the Hejaz Mountains in the Red Sea that have been identified as dust hot spots by Anisimov et al., 2017).
In Fig. 2a we show the static sources in the original model version with a factor of 0 to 1 depending on the source area strength.Accordingly in Fig. 2b  mance of the new methodology, we run the model in two different configurations: (1) using the static Ginoux et al. (2001) dust source database, called DREAM-CTRL run from now on, and (2) using the dynamic NDVI database as described above, called DREAM-NDVI run from now on.Both setups are initialized using the National Centers for Environmental Prediction Global Forecast System (NCEP GFS) analysis files (0.5 • × 0.5 • at 00:00, 06:00, 12:00, and 18:00 UTC), which were used for boundary conditions as well.The two model configurations are identical except for the dust source database.

NDVI description
For the purposes of our study we used the 500 m 16-day averaged NDVI from MODIS (Didan, 2015) for the period of interest.The NDVI is a normalized transform of the nearinfrared to red reflectance ratio, designed to provide a standard for vegetation, and takes values between −1 and +1.Since it is expressed as a ratio, the NDVI has the advantage of minimizing certain types of band-correlated noise (positively correlated) and influences attributed to variations in irradiance, clouds, atmospheric attenuation, and other parameters (Solano et al., 2010).
To create an accurate time-dependent dust source map, we have utilized the NDVI derived from the MODIS-Terra instrument.NDVI is calculated as the normalized difference of reflectance in the red and near-infrared channels (Rouse Jr. et al., 1974;Huete et al., 2002) where X represents surface reflectance as would be measured at ground level (i.e., corrected for atmospheric gas and aerosol effects) in each channel.The 16-day composite is calculated by ingesting two 8-day composite surface reflectance granules, taking into account pixel quality, presence of clouds, and viewing geometry.This procedure can lead to spatial discontinuities, as it is possible that data from different days are used for adjacent pixels, each representing different measurement conditions.If a pixel had no useful measurements during the 16-day period, historic data are used as fill values (Didan et al., 2015).For terrestrial targets, NDVI will take values near 0.8 for vegetated areas and near 0 for barren soil (Huete et al., 1999).The high-resolution dataset was used to calculate the percentage of barren land in each 0.1 • × 0.1 • model grid cells, and this percentage was used to define the effective strength of dust sources in each cell.

Evaluation datasets and metrics
Model evaluation is carried out with two datasets.First, the MODIS monthly aerosol optical depth (AOD) is use to study the spatial distribution of dust in the model domain.For this we use the Level 3 gridded atmosphere monthly product at 1 × 1 resolution, MOD08_ME (Platnick et al., 2017).
Secondly, we evaluate model performance using AERONET AOD retrievals at eight photometric stations.AERONET is a network of sun or sky photometers that derive aerosol optical and microphysical properties at a large number of stations around the world (Holben et al., 1998).For this evaluation, we use Version 3 AOD retrievals that, in comparison with previous versions, improve automatic cloud screening (Giles et al., 2018).Level 2 datasets were used for all stations apart from Kuwait University, where only Level 1.5 data were available.Both model and AERONET AODs were calculated at 532 nm; this was chosen to facilitate future intercomparing with lidar systems that frequently measure at this wavelength (e.g., Pappalardo et al., 2014).AERONET measurements were converted to this wavelength using the 440-870 Ångström exponent and taking into account AOD measurements at 440, 675, and 870 nm; in the cases where the 440 nm AOD was not available, the 500 nm (Mezaira) or 443 nm (KAUST campus) measurement was used instead.
We evaluation model performance using five metrics: mean bias, root mean square error, correlation coefficient, mean fractional bias, and fractional gross error.Concretely, assuming we have n pairs of model values (m i ) and observations (o i ), the mean bias (MB) is defined as where the bar denotes the mean value.The root mean square error (RMSE) is defined as The correlation coefficient (r) is defined as The fractional gross error (FGE) is defined as following Boylan and Russell (2006).Similarly, mean fractional bias (MFB) is defined as following Chang and Hanna (2004).

Results
The test simulation period is 1-31 August 2016 and the results from both simulations are compared to MODIS and AERONET AOD.A 5-day spin-up model run, prior to the experimental period, is used for establishing the dust background over the domain.After finalizing the experimental model configuration, we perform a complete 1-year run (2016) and evaluate the results against AERONET stations.

Dust transport during August 2016
The selected 1-month period is characterized by a significant variability in wind speeds and directions (Fig. 3), which allows the evaluation of the new model version under different conditions.During 1-10 August, east winds prevail over the region and increased dust concentrations are found mostly along the central, east and south coastal areas of the Arabian Peninsula.An anticyclonic circulation is established during 10-15 August over the Arabian Desert and increased dust concentrations are mostly found over the central desert areas.
On 16-26 August the circulation is mainly from northerly directions and thick dust plumes are advected southwards towards the Arabian Sea.The north winds veer east on 26-31 August, and increased dust loads are found over the Gulf on these dates.

Comparison with MODIS and AERONET
The monthly average AOD for August 2016 is shown in Fig. 4 for the two experimental runs (Fig. 4a, b).The DREAM-NDVI run results in a significantly modified spatial distribution of dust presenting increased dust loads over the entire domain and most profoundly over the Red Sea and Gulf regions (Fig. 4b).This dust pattern is closer to the MODIS-observed AOD over the same period that is shown in Fig. 4c.The MODIS AOD in this area is mostly related to dust; however, it must be taken into account that other aerosols not parameterized in the model (e.g., sea salt, sulfates, nitrates) may also contribute to the observed MODIS AOD.The first step is to examine how our methodology compares against the monthly average AOD in our study area.Therefore, the monthly average AOD values produced from our two simulations (DREAM-NDVI run and DREAM-CTRL run) are compared.More specifically, the DREAM-NDVI run reproduces the MODIS-observed AOD pattern that is in general characterized by values of 0.3-0.4 in the NW parts of the Arabian Peninsula and by values of 0.4-0.8 in the SE parts.Significant improvement is also evident over the Red Sea and NE Africa.The DREAM-NDVI run captures the maximum observed AOD values reaching up to 1.6 over the Red Sea and also the southwesterly extension of an AOD tongue of 0.3-0.8towards Sudan.In the eastern parts of the modeling domain, the DREAM-NDVI run again outperforms the DREAM-CTRL run since it reproduces the spatial distribution of AOD 0.4-0.8 over the Arabian Sea and the maximum of 0.8-1.2 at the SE edge of Arabian Peninsula.Over the Gulf, the NDVI run correctly represents the 0.4-0.8AOD but the dust concentration is overpredicted at the Strait of Hormuz and along the Iran-Pakistan coastline.This is mostly due to the prevailing NE winds during the last days of the August 2016 modeling period and due to a possible misclassification of Iran and Pakistan grid points as effective dust sources, thus favoring unrealistic southeasterly transport towards the Gulf of Oman.The DREAM-NDVI AOD is also higher than MODIS AOD over western Saudi Arabia, indicating a possible overprediction of dust sources in this area.
As a second step we run the same model configurations (CTRL and NDVI) for the entire 2016.The modeled dust optical depth is compared with individual AERONET measurements.The model retrievals are interpolated in time to match the AERONET measurement time considering only dust-relevant measurements with an Ångström coefficient < 0.6 (Holben et al., 1998), and the results are shown in Table 1.For completeness we first consider all AERONET stations inside the modeling domain for the evaluation.However the stations that are at the margins of our domain (Cairo_EMA_2, SEDE_BOKER, AgiaMa-rina_Xyliatou, and El_Farafra) are also affected by other dust source areas (e.g., the Sahara) and their statistics are not representative for Arabian and Middle Eastern sources.Instead, the comparison with Arabian Peninsula stations (Eilat, Kuwait_University, KAUST_Campus, and Mezaira) provides more insight on the effects of the new source characterization.As seen in Fig. 5 and also in Table 2, these stations clearly benefited from the experimental run.
In general the two runs present a significant statistical difference and, more remarkably, a reverse of bias (MODEL-AERONET) from negative in the DREAM-CTRL run to positive in the DREAM-NDVI run.The DREAM-NDVI run produces increased AODs that are neither linearly proportional to the DREAM-CTRL run AODs nor uniformly distributed over the domain.When considering only Arabian stations, the statistical metrics in Table 1 and especially the fractional gross error and bias are improved but the RMSE is increased due to the increase in maximum modeled AODs.In order to investigate the sensitivity of our results towards the severity of dust events we further assume two additional air quality states in Table 1: (i) dust events (AOD > 0.25) and (ii) severe dust episodes (AOD > 1).Both cases show an improvement in the bias values over the control simulations.When we consider AOD > 1, the DREAM-NDVI run still underestimates the observed values but with a lower RMSE (0.586 versus 0.983 of the DREAM-CTRL run).This is clearly evident in Fig. 6, where the NDVI run is indeed more realistic for the Arabian stations but still does not reproduce the extreme AOD during severe episodes.For most of the cases, such high AODs should be attributed to dust storms from convective downdrafts (haboobs).These processes are not resolved at mesoscale model resolutions (Solomos et al., 2012(Solomos et al., , 2017;;Vukovic et al., 2014) and thus cannot be represented here.

Summary and discussion
In this study we present the development of a dynamic dust source map for implementation in NMME-DREAM v1.0 over the Arabian Peninsula and the greater areas of the Middle East, SW Asia, and NE Africa.Although the major dust sources worldwide are located in permanent deserts where the NDVI is almost always < 0.1 (e.g., Bodélé Depression, Gobi, Arabian Desert), the dynamical scaling of dust emissions presented here can be important for providing up-to-date evidence of active dust sources over nonpermanent deserts.These may include dried bogs, marshes, and semidesert areas as well as irrigated and nonirrigated farms where land use changes occur throughout the year.Analysis of the modeling results for a 1-year test period (2016) over SW Asia indicated the improved performance of the new parameterization.The DREAM-NDVI run showed a significant increase in dust loads over the greater Arabian Peninsula area and a more realistic representation of the spatial distribution of AOD compared to the corresponding MODIS satellite re- trievals.These findings support the previous results by Kim et al. (2013), who also showed an increase in dust emissions and a more realistic comparison with satellite observations in Saudi Arabia through the introduction of an NDVI-based dynamic source mapping for the GOCART model.Comparison with AERONET measurements also showed significant improvement especially at higher AODs that are also relevant to the model efficiency for air quality purposes (i.e., the model bias is reduced from −0.140 to 0.083 at AOD > 0.25 and from −0.933 to −0.424 at AOD > 1).However, the model statistics are not improved for all AERONET measuring stations and for all air quality states (Table 2), mainly due to a possible misclassification of dust sources in the highlands of Iran and Pakistan.
The main purpose of our work was the development and first testing of this new modeling version.A major advance of our study is the ability to implement the real-time properties of dust sources in air quality simulations (as represented   by the satellite NDVI) and thus capture local or seasonal effects.In general, 1 year is not sufficient for extracting robust statistical results and further analysis is required to examine the performance of the proposed methodology over longer time periods and also over different areas worldwide.For example the simple approach of employing a uniform value of NDVI < 0.1 for determining the active dust sources may not be adequate to represent fine-scale land properties, and further adjustments may be required depending on local-scale characteristics.This new approach for the dynamic characterization of active dust sources based on NDVI can be easily implemented in other atmospheric dust models at different configurations and spatial coverage to improve their performance.

Figure 1 .
Figure 1.DREAM model domain and topography in meters.
we show the new dynamic sources for 1-16 August 2016.The two dust source patterns present remarkable differences especially over western Saudi Arabia and over Iran and Pakistan where the NDVI classification results in stronger emissions.In order to test the perforwww.geosci-model-dev.net/12/979/2019/Geosci.Model Dev., 12, 979-988, 2019

Figure 4 .
Figure 4. Monthly average simulated AOD during August 2016 from the DREAM-CTRL run (a), the DREAM-NDVI run (b) and (c) MODIS.

Figure 5 .
Figure 5. Density scatterplots of modeled and AERONET dust AOD at the stations of Mezaira, KAUST, Eilat, and Kuwait for 2016.

Figure 6 .
Figure 6.Time series of measured and modeled dust AOD for the cases of AERONET AOD > 1.

Table 2 .
Statistical metrics from the annual runs (2016) at AERONET stations.Bold values indicate a correlation coefficient with p < 0.01.