Assessment of natural and anthropogenic aerosol air pollution in the Middle East using MERRA-2, CAMS data assimilation products, and high-resolution WRF-Chem model simulations

Modern-Era Retrospective analysis for Research and Applications v.2 (MERRA-2), Copernicus Atmosphere Monitoring Service Operational Analysis (CAMS-OA), and a high-resolution regional Weather Research and Forecasting model coupled with chemistry (WRF-Chem) were used to evaluate natural and anthropogenic particulate matter (PM) air pollution in the Middle East (ME) during 2015–2016. Two Moderate Resolution Imaging Spectrometer (MODIS) retrievals – combined product Deep Blue and Deep Target (MODIS-DB&DT) and Multi-Angle Implementation of Atmospheric Correction (MAIAC) – and Aerosol Robotic Network (AERONET) aerosol optical depth (AOD) observations as well as in situ PM measurements for 2016 were used for validation of the WRF-Chem output and both assimilation products.

MERRA-2 and CAMS-OA assimilate AOD observations. WRF-Chem is a free-running model, but dust emission in WRF-Chem is tuned to fit AOD and aerosol volume size distributions obtained from AERONET. MERRA-2 was used to construct WRF-Chem initial and boundary conditions both for meteorology and chemical and aerosol species. SO₂ emissions in WRF-Chem are based on the novel OMI-HTAP SO₂ emission dataset.

The correlation with the AERONET AOD is highest for MERRA-2 (0.72–0.91), MAIAC (0.63–0.96), and CAMS-OA (0.65–0.87), followed by MODIS-DB&DT (0.56–0.84) and WRF-Chem (0.43–0.85). However, CAMS-OA has a relatively high positive mean bias with respect to AERONET AOD. The spatial distributions of seasonally averaged AODs from WRF-Chem, assimilation products, and MAIAC are well correlated with MODIS-DB&DT AOD product. MAIAC has the highest correlation (R=0.8), followed by MERRA-2 (R=0.66), CAMS-OA (R=0.65), and WRF-Chem (R=0.61). WRF-Chem, MERRA-2, and MAIAC underestimate and CAMS-OA overestimates MODIS-DB&DT AOD.

The simulated and observed PM concentrations might differ by a factor of 2 because it is more challenging for the model and the assimilation products to reproduce PM concentration measured within the city. Although aerosol fields in WRF-Chem and assimilation products are entirely consistent, WRF-Chem is preferable for analysis of regional air quality over the ME due to its higher spatial resolution and better SO₂ emissions. The WRF-Chem’s PM background concentrations exceed the World Health Organization (WHO) guidelines over the entire ME. Mineral dust is the major contributor to PM (≈75 %–95 %) compared to other aerosol types. Near and downwind from the SO₂ emission sources, nondust aerosols (primarily sulfate) contribute up to 30 % to PM_2.5. The contribution of sea salt to PM in coastal regions can reach 5 %. The contributions of organic matter, black carbon and organic carbon to PM over the Middle East are insignificant. In the major cities over the Arabian Peninsula, the 90th percentile of PM₁₀ and PM_2.5 (particles with diameters less than 10 and 2.5 µm, respectively) daily mean surface concentrations exceed the corresponding Kingdom of Saudi Arabia air quality limits. The contribution of the nondust component to PM_2.5 is <25 %, which limits the emission control effect on air quality. The mitigation of the dust effect on air quality requires the development of environment-based approaches like growing tree belts around the cities and enhancing in-city vegetation cover. The WRF-Chem configuration presented in this study could be a prototype of a future air quality forecast system that warns the population against air pollution hazards.

1 Introduction

Particulate matter (PM) is a complex mixture of sea salt, sulfate, black carbon, organic matter, and mineral dust, suspended in the air. The dramatic increase in the level of air pollution in developing countries over the last decades is forced by rapid economic and population growth, burning of fossil fuels, construction, and agricultural activities (Janssens-Maenhout et al., 2015). However, the primary cause of air pollution in the Middle East (ME) is mineral dust, and it is on the rise (Klingmüller et al., 2016). Along with Asia and Africa, the ME significantly contributes to global dust emissions, which are in the range of 1000–2000 Tg yr⁻¹ (Zender et al., 2004). According to Prospero et al. (2002), the Middle East and North Africa (MENA) regions account for about half of global dust emissions. By integrating surface emissions in MERRA-2 reanalysis, we found that the total global dust emission averaged over the 2015–2016 period is about 1600 Tg yr⁻¹, right in the middle of the Zender et al. (2004) estimate. The dust emission from our simulation domain (see Fig. 1) that covers the ME and nearby areas is about 500 Tg yr⁻¹, contributing ≈30 % to the global dust emission budget. Furthermore, frequent inflows of pollutants from Europe, Africa, and India worsen the air quality over the Arabian Peninsula (Jish Prakash et al., 2015; Kalenderski et al., 2013; Notaro et al., 2013; Reid et al., 2008; Mohalfi et al., 1998; Kalenderski and Stenchikov, 2016; Parajuli et al., 2019). Because of the large amount of dust, the ME is one of the most polluted areas in the world. Located in the center of the northern subtropical dust belt, the Arabian Desert is the third-largest (after the Sahara and the East Asian deserts) region of dust generation, where dust plays a significant role in controlling regional climate (Cahill et al., 2017; Banks et al., 2017; Jish Prakash et al., 2016; Farahat, 2016; Kalenderski and Stenchikov, 2016; Munir et al., 2013; Alghamdi et al., 2015; Lihavainen et al., 2016; Anisimov et al., 2017; Osipov and Stenchikov, 2018).

Figure 1 Simulation domain with marked locations of the air quality monitoring stations (AQMSs) and Aerosol Robotic Network (AERONET) sites.

In addition to natural dust aerosols, the ME receives high concentrations of anthropogenic PM (Karagulian et al., 2015; Al-Taani et al., 2019; Alharbi et al., 2015; Khodeir et al., 2012). The most important anthropogenic aerosol in ME is sulfate, with SO₂ as a precursor; the contributions of other types of aerosols in PM – sea salt, organic matter, and black carbon – are of lesser importance in the ME (Randles et al., 2017). The ME emits about 10 % of the total global anthropogenic SO₂ (Klimont et al., 2013). SO₂ produced in the course of power generation, water desalination, and oil recovery operations (Al-Jahdali and Bisher, 2008) is converted photochemically into sulfate aerosol, which contributes to PM and has significant adverse effects on human health (Lelieveld et al., 2015). Ukhov et al. (2020b) simulated SO₂ transport and distribution over the Middle East using the high-resolution WRF-Chem (Weather Research and Forecasting model coupled with chemistry) model and demonstrated high surface concentrations of SO₂ along the west and east coasts of the Arabian Peninsula.

The impact of aerosols on air quality is characterized by near-surface concentrations of PM, which comprise both PM₁₀ and PM_2.5 (particles with diameters less than 10 and 2.5 µm, respectively). Extended exposure to PM may cause cardiovascular and respiratory disease, lung cancer, and premature mortality on a global scale (Lelieveld et al., 2015). According to the WHO, outdoor air pollution caused 4.2 million premature deaths worldwide in 2016 (WHO, 2018). To protect human health and the environment, the WHO (WHO, 2006) and the National Agencies, e.g., the United States Environmental Protection Agency (US EPA; USEPA, 2010), European Commission (EC; EUEA, 2008), and Kingdom of Saudi Arabia Presidency of Meteorology and Environment (KSA-PME; PME, 2012), issued the air quality regulations for PM that are presented in Table 1. The WHO regulations are the strictest, while the KSA-PME regulations are the softest.

Table 1 Air quality regulations for PM_2.5 and PM₁₀ prescribed by the WHO, US EPA, EC, and KSA-PME (µg m⁻³).

¹ 98th percentile, averaged over 3 years. ² Annual mean, averaged over 3 years. ³ 35 permitted exceedances per year. ⁴ Not to be exceeded more than once per year on average over 3 years.

Global satellite observations of aerosol optical depth (AOD) inform about vertically integrated aerosol loading in an entire atmospheric column. But the near-surface PM concentration cannot be observed from space. These measurements could be conducted only in situ in a limited number of locations. Along with instrumental observations, modern data assimilation products provide valuable information about AOD and near-surface PM concentration even in areas where satellite sensors are unreliable due to factors such as the high reflectivity of land surfaces (Shi et al., 2011). Assimilation products improve the aerosol total column loadings through the assimilation of observed AOD but are not capable of assimilating the aerosol vertical structure and chemical composition. There are two well-known data assimilation products that assimilate atmospheric constituents: MERRA-2 (Randles et al., 2017; Buchard et al., 2017) from the National Aeronautics and Space Administration (NASA) Goddard Space Flight Center (GSFC), and the Copernicus Atmosphere Monitoring Service Operational Analysis (CAMS-OA; Inness et al., 2019a; Flemming et al., 2015; Inness et al., 2015), from the European Centre for Medium-Range Weather Forecasts (ECMWF). These data assimilation products adequately reproduce AOD and PM concentrations in different regions of the world (Provençal et al., 2017; Buchard et al., 2017; Cesnulyte et al., 2014; Cuevas et al., 2015). Provençal et al. (2017), for example, tested PM surface concentrations from the MERRA Aerosol Reanalysis (predecessor of MERRA-2) against observations over Europe. Buchard et al. (2017) evaluated MERRA-2 surface PM_2.5 on the global scale and over the continental United States. Excessive validation of the Monitoring Atmospheric Composition and Climate (MACC) reanalysis (predecessor of CAMS) has been conducted by Cesnulyte et al. (2014), who compared the model AOD with the Aerosol Robotic Network (AERONET) observations. Cuevas et al. (2015) evaluated atmospheric mineral dust from the MACC reanalysis over the MENA region for 2007–2008 using satellite and ground-based observations. MERRA-2 and CAMS-OA are global and have a relatively low spatial resolution (in comparison with the regional models), which diminishes their ability to resolve fine-scale regional spatial features. Like any other model, MERRA-2 and CAMS-OA use emission inventories of anthropogenic pollutants that may be outdated and incomplete, especially in rapidly developing parts of the world, like the ME region (McLinden et al., 2016). SO₂ emissions used in MERRA-2 and CAMS-OA, for example, differ by 45 %–50 % in some ME regions (Ukhov et al., 2020b).

In this study, we evaluate aerosol outputs from MERRA-2, CAMS-OA, and WRF-Chem over the ME against satellite and ground-based AOD observations as well as in situ PM_2.5 and PM₁₀ measurements during the 2015–2016 period, and we assess air pollution over the ME focusing on the following science questions:

1. How accurately do WRF-Chem, MERRA-2, and CAMS-OA capture the abundance of dust aerosol, its volume size, and its spatial distributions over the ME in comparison with AERONET and satellite observations?

2. How accurately do WRF-Chem, MERRA-2, and CAMS-OA capture PM surface concentrations compared with in situ measurements?

3. What are the contributions of dust, sea salt, sulfate, black carbon, and organic matter in PM surface concentrations?

4. What is the overall impact of PM pollution on air quality over the ME region and in the major ME cities?

The paper is organized as follows: Sect. 2 describes the observational datasets used in this study. Section 3 characterizes data assimilation products. In Sect. 4, the WRF-Chem model setup is described. In Sect. 5, the capabilities of WRF-Chem, MERRA-2, and CAMS-OA to simulate dust aerosol abundance over the ME are compared; the PM spatial distributions and PM air pollution in the major ME cities obtained from the WRF-Chem simulations are also discussed. Conclusions are presented in Sect. 6.

2 Observational datasets

To evaluate the data assimilation products and WRF-Chem output, we use Moderate Resolution Imaging Spectrometer (MODIS) AOD retrievals, ground-based AERONET AOD observations, and aerosol volume size distribution retrievals as well as in situ measurements of PM surface concentrations.

2.1 AERONET

AERONET comprises more than 1000 observation sites equipped with CIMEL sunphotometers and PREDE skyradiometers manufactured in France by CIMEL and in Japan by PREDE. They measure direct sun and sky radiances at eight wavelengths (340, 380, 440, 500, 670, 870, 940, and 1020 nm) every 15 min during daylight (Holben et al., 1998). In 2012 we established the KAUST Campus site, which is currently the only permanently operational AERONET site in Saudi Arabia. For this study we have chosen three AERONET sites (KAUST Campus, Mezaira, and Sede Boker; see Fig. 1) that routinely collected data in 2015–2016 and are located within our domain. We utilized level 2.0 (cloud-screened and quality-assured) AERONET AOD data. To facilitate comparison with the model output, the 550 nm AOD is calculated using the following relation:

$$\begin{array}{}\text{(1)}& {\displaystyle \frac{{\mathit{\tau }}_{\mathit{\lambda }}}{{\mathit{\tau }}_{{\mathit{\lambda }}_{\mathrm{0}}}}}={\left({\displaystyle \frac{\mathit{\lambda }}{{\mathit{\lambda }}_{\mathrm{0}}}}\right)}^{-\mathit{\alpha }},\end{array}$$

where α is the Ångström exponent for the 440–675 nm wavelength range provided by AERONET, τ_λ is the optical thickness at wavelength λ, and ${\mathit{\tau }}_{{\mathit{\lambda }}_{\mathrm{0}}}$ is the optical thickness at the reference wavelength λ₀. From here forward, we presume that AOD is given or calculated at 550 nm.

In addition to direct observations of AOD, the AERONET retrieval algorithm provides column-integrated aerosol volume size distribution (AVSD) dV∕dln r (µm³ µm⁻²) on 22 logarithmically equidistant discrete points in the range of radii between 0.05 and 15 µm (Dubovik and King, 2000). We use these retrievals to evaluate the AVSDs calculated by WRF-Chem, CAMS-OA, and MERRA-2.

2.2 MODIS

MODIS instruments on-board the NASA TERRA and AQUA satellites provide aerosol properties over both land and ocean with near-daily global coverage. The standard MODIS AOD aerosol product combines two retrieval algorithms: (1) the MODIS Dark Target (DT) algorithm (Kaufman et al., 1997) is used over the ocean and dark areas with sufficient vegetation; (2) the Deep Blue (DB) algorithm is used over bright desert surfaces of the Sahara and the ME. From this combined product (MODIS-DB&DT v6.1), we use level 2 data for AOD at 550 nm from the daily dataset at 10 km spatial resolution, downloaded from https://ladsweb.modaps.eosdis.nasa.gov/about/purpose (last access: 5 January 2020) (Levy and Hsu, 2015).

Recently, a new MODIS AOD product became available that was obtained using the Multi-Angle Implementation of Atmospheric Correction (MAIAC) algorithm (Lyapustin et al., 2018). This algorithm uses time series analysis and image processing to derive the surface bidirectional reflectance function at fine spatial resolution. MAIAC uses empirically tuned, spatially varying aerosol properties derived from the AERONET climatology and provides AOD at 470 and 550 nm with 1 km spatial resolution over land globally. We include the new MAIAC product (version 6, level 2) in the comparison between simulated and retrieved AODs.

2.3 Surface in situ PM observations

To test the model-produced PM concentrations, we use observations conducted by the air quality monitoring stations (AQMSs) that measure surface concentrations of PM_2.5 and PM₁₀ in Riyadh, Jeddah, and Dammam (megacities of Saudi Arabia; see Fig. 1). Observations are available starting from 2016. The measurements were conducted by the Saudi Authority for Industrial Cities and Technology Zones (MODON). MODON uses an MP101M analyzer to continuously detect PM_2.5 and PM₁₀ concentrations by measuring the absorption of low-energy β radiation that is proportional to the mass of aerosol particles independently of their physicochemical nature (measurement method ISO 10473). The PM_2.5 and PM₁₀ measurement error is ±5 %. The system satisfies the European Standard EN 12341 and US EPA 40CFR part 53 for the continuous monitoring of PM₁₀ and EN 14907 for the continuous monitoring of PM_2.5. The PM measurements are conducted every 15 min, and collected data are transmitted in real-time to servers at MODON for processing and storage. To provide confidence in the operational status of each AQMS, a comprehensive physical audit is conducted quarterly by Ricardo-AEA Ltd, (https://www.ctc-n.org/network/network-members/ricardo-aea-ltd, last access: 5 January 2020).

3 Data assimilation products

MERRA-2 and CAMS-OA assimilate satellite observations to provide aerosol abundance and air quality data globally. MERRA-2 also assimilates AERONET AODs. In contrast, WRF-Chem is a free-running model and does not assimilate observations.

3.1 MERRA-2

MERRA-2 (https://gmao.gsfc.nasa.gov/reanalysis/MERRA-2, last access: 5 January 2020) provides meteorological and atmospheric composition fields on $\mathrm{0.625}{}^{\circ }\times \mathrm{0.5}{}^{\circ }$ latitude–longitude grids and 72 terrain-following hybrid σ–p model layers (Randles et al., 2017; Buchard et al., 2017). The pressure at the model top equals 0.01 hPa. MERRA-2 uses the Goddard Earth Observing System version 5 (GEOS-5) atmospheric model (Rienecker et al., 2008), which is interactively coupled to the Goddard Global Ozone Chemistry Aerosol Radiation and Transport (GOCART) model Chin et al., 2002, 2014; i.e., it takes into account the effects of aerosols on radiation and model dynamics. This model simulates dust and sea salt in five size bins (see Table 2), SO₂, sulfate, organic and black carbon (hydrophobic and hydrophilic), O₃, CO, dimethyl sulfide (DMS), and methane sulfonic acid (MSA). The dust density is 2600 kg m⁻³ for all sizes. Dust and sea salt emissions are calculated in the model depending on the near-surface wind. The dust source function is taken from Ginoux et al. (2001). For anthropogenic emissions, MERRA-2 employs the EDGAR-4.2 (Janssens-Maenhout et al., 2013) emission inventory available on a $\mathrm{0.1}{}^{\circ }\times \mathrm{0.1}{}^{\circ }$ grid. MERRA-2 assimilated AOD at 550 nm from the Advanced Very High-Resolution Radiometer (AVHRR; Heidinger et al., 2014) over the oceans until 2002. Since 2000 MERRA-2 has been assimilating MODIS and Multiangle Imaging Spectroradiometer (MISR; Kahn et al., 2005) data over land and ocean. Both instruments are on the TERRA satellite, which has an equatorial overpass at 10:30 UTC, while the AVHRR has mostly orbited with an afternoon equatorial crossing time. Therefore MERRA-2 continued using AVHRR data over the ocean until 2002, when the AQUA satellite was launched. Since AQUA has an orbit with the Equator overpass at 14:30 UTC, AVHRR data were no longer needed for coverage. We have to mention that MERRA-2 assimilates specially processed MODIS observations, not the standard MODIS-DB&DT aerosol product. It also assimilates AERONET AODs (Randles et al., 2017).

Table 2 Radii ranges (µm) of dust and sea salt bins used in the GOCART model (WRF-Chem, MERRA-2) and in CAMS-OA.

3.2 CAMS-OA

CAMS-OA (https://atmosphere.copernicus.eu/, last access: 5 January 2020) has been conducted in almost real time since July 2012. The CAMS-OA product had a resolution of $\mathrm{0.8}{}^{\circ }\times \mathrm{0.8}{}^{\circ }$ before 21 June 2016 and $\mathrm{0.4}{}^{\circ }\times \mathrm{0.4}{}^{\circ }$ after that, with 60 vertical levels. It employs the ECMWF aerosol data assimilation system developed within the Integrated Forecast System (IFS; Morcrette et al., 2009; Benedetti et al., 2009). The extended version of the Carbon Bond chemical mechanism 5 (CB05; Yarwood et al., 2005) is implemented in the IFS (Flemming et al., 2015). CB05 describes tropospheric chemistry with 54 species and 126 reactions. The chemistry scheme is coupled with the aerosol module.

CAMS-OA simulates five aerosol species: dust, sea salt, sulfate, organic carbon, and black carbon. To calculate dust and sea salt, it uses three dust bins (see Table 2). The dust density is 2600 kg m⁻³ for all bins. Emissions of mineral dust and sea salt depend on simulated near-surface wind speed. Dust emission is parameterized following Marticorena and Bergametti (1995), with the source function adopted from Ginoux et al. (2001). SO₂ oxidation into sulfate aerosol is parameterized using a prescribed latitude-dependent e-folding timescale ranging from 3 d at the Equator to 8 d at the poles. The anthropogenic emissions for the chemical species are taken from the MACCity inventory (Granier et al., 2011), which is available on a $\mathrm{0.5}{}^{\circ }\times \mathrm{0.5}{}^{\circ }$ grid and covers the period 1960–2010. CAMS-OA assimilates MODIS AQUA and TERRA AODs. It has been using observations from Collection 5 since 2009 and Deep Blue since 2015.

4 WRF-Chem

To calculate fine-resolution PM and sulfate fields, we use the Weather Research and Forecasting (WRF) model (Skamarock et al., 2005) coupled with chemistry (WRF-Chem v3.7.1; Grell et al., 2005). The WRF-Chem is used for prediction and simulation of weather, air quality, and dust storms, accounting for the aerosol effect on radiation. WRF-Chem can be configured with one of the few gas-phase chemical mechanisms, photolysis, and aerosol parameterization models. WRF-Chem has been widely used for air quality simulations in different parts of the globe: East Asia (Wang et al., 2010), North America (Kim et al., 2006; Chuang et al., 2011), Europe (Forkel et al., 2012; Ritter et al., 2013), South America (Archer-Nicholls et al., 2015), and the Middle East (Parajuli et al., 2019). To reduce the clock time of our 2-year calculations, we simulated each month of the 2015–2016 period separately. Each simulation starts from the last week of the previous month. This time is considered a spin-up and is excluded from postprocessing. The simulation domain, shown in Fig. 1, is centered at 28^∘ N, 42^∘ E, and a 10 km×10 km horizontal grid (450×450 grid nodes) is employed. The vertical grid comprises 50 vertical levels with enhanced resolution closer to the ground, comprising 11 model levels within the near-surface 1 km layer. The model top boundary is set at 50 hPa.

To improve the representation of the meteorological fields, we apply spectral nudging (Miguez-Macho et al., 2004) above the planetary boundary layer (PBL; >5.0 km) to horizontal wind components (U and V) toward the MERRA-2 wind field. The nudging coefficient for U and V is set to be 0.0001 s⁻¹. We only nudge waves with wavelengths longer than 450 km. This allows us to keep the large-scale motions close to reanalysis and leave the resolved small-scale, high-frequency features unaffected.

The aerosol and chemistry initial and boundary conditions (IC&BC) are calculated using MERRA-2 output using the newly developed Merra2BC interpolation utility (Ukhov and Stenchikov, 2020). To be consistent with aerosol and chemistry IC&BC, we also define the meteorological IC&BC using MERRA-2 output (see Appendix A1).

The following set of physical parameterizations was used in WRF-Chem runs. The Unified Noah land surface model (sf_surface_physics=2) and the Revised MM5 Monin-Obukhov scheme (sf_sfclay_physics=1) are chosen to represent land surface processes and surface layer physics. The Yonsei University scheme is chosen for PBL parameterization (bl_pbl_physics=1). The WRF single-moment microphysics scheme (mp_physics=4) is used for the treatment of cloud microphysics. The New Grell scheme (cu_physics=5) is used for cumulus parameterization. The Rapid Radiative Transfer Model (RRTMG) for both shortwave (ra_sw_physics=4) and longwave (ra_lw_physics=4) radiation is used for radiative transfer calculations. Only the aerosol direct radiative effect is accounted for. More details on the physical parameterizations used can be found at http://www2.mmm.ucar.edu/wrf/users/phys_references.html (last access: 5 January 2020).

4.1 Gas-phase chemistry and aerosols

To calculate the atmospheric chemistry within WRF-Chem, we employ the Regional Atmospheric Chemistry Mechanism (RACM; Stockwell et al., 1997), containing 77 species and 237 reactions, which include 23 photolysis reactions but no heterogeneous chemistry. The RACM chemical module is embedded into WRF-Chem using the Kinetic PreProcessor (KPP; Damian et al., 2002). The role of the KPP is to integrate the system of stiff, nonlinear ordinary differential equations, which represents the specified set of chemical reactions. The photolysis rates are calculated on-line according to Madronich (1987; phot_opt=1). Similar to MERRA-2, the GOCART chemistry module is used to calculate SO₂ to sulfate oxidation (Chin et al., 2002, 2014) by the hydroxyl radical OH, whose abundance is interactively simulated by the RACM.

We use the novel OMI-HTAP SO₂ emission dataset (Liu et al., 2018) based on the combination of distributed SO₂ emissions from residential and transportation sectors, taken from the HTAP-2.2 inventory (Janssens-Maenhout et al., 2015) with the catalog of the strong (>30 kt yr⁻¹) SO₂ point emissions (Fioletov et al., 2016) built using satellite observations by the Ozone Monitoring Instrument (OMI; Levelt et al., 2006; Li et al., 2013). The catalog contains more than 500 point sources of industrial origin, some of which are not present in the widely used EDGAR-4.2 and HTAP-2.2 emission datasets. For example, 14 SO₂ point emissions located in the ME that were previously unaccounted for (mostly in the Arabian Gulf) were detected, most of which are related to the oil and gas industry. OMI-HTAP divides SO₂ emissions into surface and elevated ones. We distribute the surface SO₂ emissions with a constant mixing ratio in the 0–1000 m layer and elevated emissions in the 120–1000 m layer. All other constituents (other PM from biogenic and fossil components, black and organic carbon, etc.), including SO₂ shipping emissions, are taken from the HTAP-2.2 inventory and are treated as surface emissions. OMI-HTAP SO₂ emissions are provided on a $\mathrm{0.1}{}^{\circ }\times \mathrm{0.1}{}^{\circ }$ grid (Liu et al., 2018). We conservatively interpolated them on the WRF-Chem 10 km×10 km grid. See Ukhov et al. (2020b) for details.

To calculate aerosols we employ the GOCART (Chin et al., 2002) aerosol model (chem_opt=301). It is the same microphysical model as that used in MERRA-2 (see Sect. 3.1). Dust and sea salt size distributions in WRF-Chem are approximated by the same five dust and sea salt size bins as those in MERRA-2 (Table 2). However, only the last four “salt” bins in Table 2 are used in WRF-Chem as the first bin appears to be very poorly populated. Dust density is assumed to be 2500 kg m⁻³ for the first dust bin and 2650 kg m⁻³ for dust bins 2–5. Emission of sea salt is calculated according to Gong (2003). Dust emission from the surface is calculated using the GOCART emission scheme Ginoux et al. (2001; dust_opt=1). Dust emission mass flux, F_p (µg m⁻² s⁻¹), in each dust bin $p=\mathrm{1},\mathrm{2},\mathrm{\dots },\mathrm{5}$ is defined by the relation

$$\begin{array}{}\text{(2)}& F_p=\left\{\begin{array}{ll}CS{s}_p{u}_{\mathrm{10}\mathrm{m}}^{\mathrm{2}}\left(u_{\mathrm{10}\mathrm{m}}-u_{\mathrm{t}}\right),& \text{if }{u}_{\mathrm{10}\mathrm{m}}>u_{\mathrm{t}}\\ \mathrm{0},& \text{otherwise}\end{array}\right.\end{array}$$

where C (µg s² m⁻⁵) is a spatially uniform factor which controls the magnitude of dust emission flux; S is the spatially varying topographic source function (Ginoux et al., 2001) that characterizes the spatial distribution of dust emissions; u_10 m is the horizontal wind speed at 10 m height; u_t is the threshold velocity, which depends on particle size and surface wetness; and s_p is a fraction of dust mass emitted into dust bin p, $\sum s_p=\mathrm{1}$.

To avoid natural dust emission in urban areas, we use the built-in WRF-Chem of the U.S. Geological Survey (USGS) 24-category land use dataset (Anderson, 1976). We modify the source function S using the following expression:

$$\begin{array}{}\text{(3)}& S^{\prime }=(\mathrm{1.0}-\mathrm{URBAN}\mathrm{\_}\mathrm{MASK})\cdot S,\end{array}$$

where S^′ is the modified topographic source function, and URBAN_MASK is the USGS “Urban and Built-up Land” mask field. It has the sense of a fraction of urban area in a grid cell and ranges from 0 to 1. Grid cells with URBAN_MASK=1 do not produce natural dust emissions. We do not account for anthropogenic dust emissions within cities, and we therefore potentially underestimate urban dust pollution.

As in our previous studies (Kalenderski et al., 2013; Jish Prakash et al., 2015; Anisimov et al., 2017), we tune dust emissions to fit the AOD from the AERONET stations located within the domain. For this purpose, the factor C from Eq. (2) has been adjusted to achieve the best agreement between simulated and observed AOD at the KAUST Campus, Mezaira, and Sede Boker AERONET sites (see Fig. 1). Assuming that factor C does not depend on time and geographical coordinates, we can only tune the annual average AOD bias. Both simulations and observations represent the total AOD with contributions from all types of aerosols. Because dust dominates all other aerosols in the ME, we choose to tune only the dust emissions. Obtained during test runs, a C value of 0.5 is kept constant in all subsequent production runs. We also tune s_p from Eq. (2) to better reproduce the AVSDs provided by the AERONET inversion algorithm. This tuning and the comparisons of AOD and AVSDs from the assimilation products and WRF-Chem simulations are discussed in detail below.

In situ air quality observations in the Middle East are scarce. It is one of the known problems for air quality research in this area. Things are simplified a bit by the fact that in the ME, dust dominates aerosol pollution. Cloud-Aerosol Lidar and Infrared Pathfinder Satellite Observations (CALIPSO; Vaughan et al., 2004), for example, records dust in 95 % of the profiles (Osipov et al., 2015). The effect of nitrates, ammonia, and organics on AOD and PMs is insignificant; therefore, the employed aerosol–chemical scheme (GOCART-RACM) is adequate for the ME conditions. To support this conclusion, we have conducted laboratory analysis of the chemical composition of soil and dust deposition samples that show little presence of organics and ammonium (Jish Prakash et al., 2016; Engelbrecht et al., 2017). According to Engelbrecht et al. (2017), in 2015, the annual average weight percentages of soluble ions of ammonium (NH₄) and sulfate in deposition samples taken at four sites at the KAUST Campus are 0.05 % and 2.513 %, respectively. This means that available ammonium may neutralize at a maximum of 5 % of sulfate mass. The actual contribution of ammonium sulfate should be lower as some ammonium may also be bound as ammonium nitrate, ammonium phosphate, or ammonium chloride.

4.2 WRF-Chem code modification

We have corrected the source code of the WRF-Chem v3.7.1 with the GOCART aerosol module in several places. These corrections were implemented in the WRF-Chem v4.1.3 official release. We evaluate how they change the results in Ukhov et al. (2020a). Here we only briefly discuss the introduced changes and their effects. Firstly, the diagnostic output of PM concentrations was corrected because contributions of the individual dust and sea salt bins were incorrectly calculated. Therefore, PM_2.5 surface concentrations were erroneously underestimated by 7 %, while PM₁₀ surface concentrations were overestimated by 5 %. Secondly, we found that the contribution of fine dust particles with radii <0.46 µm was omitted in the calculation of AOD; AOD was consequently underestimated by 25 %–30 %. This led to an overestimation of the dust emission flux because we force the simulated AOD to match the AERONET observations. Thirdly, we fixed the dust and sea salt gravitational settling subroutine since the calculations of mass fluxes of settling particles did not initially account for changes in air density. Due to this error, the total mass of dust and sea salt aerosols increased, violating mass conservation.

5 Results

5.1 Regional climate and circulation

The ME is one of the hottest and driest regions on the Earth. Summer in the ME is long and hot with little precipitation. Precipitation mainly occurs in the southwest of the Arabian Peninsula. Winter is mild, with rainfall being mostly associated with cold fronts and cyclones propagating from the eastern Mediterranean (Climate.com, 2018). Emission and transport of dust are driven by winds. Emission and deposition of dust are also sensitive to soil moisture and precipitation (Furman, 2003; Shao, 2008; Yu et al., 2015). However, because the ME is an arid region, the soil moisture and precipitation effects are insignificant.

Figure 2 Seasonally averaged (2015–2016) wind barbs (m s⁻¹) at 10 m, sea level pressure (contours), and erodibility function (shading; Ginoux et al., 2001). (a) Winter (DJF), (b) Summer (JJA).

Figure 2 shows contours of sea level pressure, topographic source function S (Ginoux et al., 2001), and seasonally (2015–2016) averaged wind speed barbs at 10 m height over the ME during winter (DJF) and summer (JJA) from WRF-Chem simulations. Over northeast Africa in winter (see Fig. 2a), the strong pressure gradient between the Red Sea trough and the stationary high-pressure system over Egypt predominantly generates moderate northeasterly winds (up to 10 m s⁻¹). Therefore in winter, dust storms occur more frequently in the west of the Arabian Peninsula. Over the central and eastern Arabian Peninsula and the eastern part of the ME, winds are relatively weak and do not have a clear direction. However, cold fronts generated by Mediterranean cyclones can cause dust storms and dust transport to central regions of the Arabian Peninsula.

In summer (see Fig. 2b) the high-pressure system over the eastern Mediterranean and low-pressure system over the Arabian Gulf promote moderate north-northwesterly winds known as shamals (Yu et al., 2016; Hamidi et al., 2013), which dominate over the central part of the Arabian Peninsula. Shamals are the primary drivers of dust storm events over this area (Yu et al., 2016; Shao, 2001; Middleton, 1986; Goudie and Middleton, 2006; Notaro et al., 2015). Shamals bring dust to the Arabian Gulf and the north and central part of Saudi Arabia from the Tigris–Euphrates basin of Syria and Iraq (Anisimov et al., 2018).

Figure 3 Seasonally averaged (2015–2016) wind speed at 10 m from WRF-Chem, MERRA-2, and CAMS-OA during winter (DJF) and summer (JJA).

Figure 3 shows wind speed seasonally averaged for 2015–2016 at 10 m from MERRA-2, CAMS-OA, and WRF-Chem during winter (DJF) and summer (JJA). WRF-Chem spatial distributions of wind speed agree well with MERRA-2 and CAMS-OA, but due to the higher spatial resolution, WRF-Chem better resolves the fine-scale spatial structures of the 10 m wind field over complex terrain. All panels have similar seasonal variations of wind speed. In winter, maximum winds are stronger over the southeast of the domain. In the central and northern parts of the domain winds are weak. In summer, wind speed increases in the northern and central parts of the ME. Somali jets produce strong (10–15 m s⁻¹) winds in the Arabian Sea along the coasts of Somalia and Oman.

Table 3 Pearson correlation coefficient R and root mean square difference (RMSD; m s⁻¹) for the seasonally averaged 2015–2016 wind components U and V at 10 m.

^∗ wrt – with respect to.

To conduct the statistical analysis, we interpolated the seasonally averaged 2015–2016 zonal and meridional wind components (U and V) at 10 m from WRF-Chem and CAMS-OA on a MERRA-2 grid and calculated the Pearson correlation coefficient (R) and root mean square differences (RMSDs) between each pair (see Table 3), respectively. RMSD is calculated using the same formula as the root mean square error (RMSE). The procedure of calculation of these parameters is given in Appendix A2. Pearson correlation coefficients provided in Table 3 are close to 1. The highest correlation is achieved between MERRA-2 and CAMS-OA. WRF-Chem's correlation coefficient with respect to MERRA-2 is smaller but exceeds that of the WRF-Chem–CAMS-OA pair. The WRF-Chem and MERRA-2 wind fields are close partly because WRF-Chem boundary conditions are built using MERRA-2 reanalysis, and the large-scale winds are nudged (see Sect. 4) to the ones from MERRA-2 over the PBL.

The RMSDs (see Table 3) are lower in winter than in summer. All RMSDs are in the range of 0.45–0.85 m s⁻¹. The lowest RMSDs are between MERRA-2 and CAMS-OA. Notably, the correlation coefficients for the meridional component V are higher, and the RMSDs are lower when compared with the zonal wind component U. This is because the northern winds are stable since they are maintained by the large-scale processes. In contrast, the zonal wind component, which is affected by small-scale processes like sea breezes, is variable. The results of the statistical analysis in Table 3 and the clear similarity of the spatial patterns (among all products) of the averaged 10 m wind fields presented in Fig. 3 suggest that WRF-Chem captures the magnitude and spatial distribution of the 10 m wind. Thus, we conclude that WRF-Chem with the selected set of physical parameterizations satisfactorily simulates both the large- and mesoscale atmospheric processes in the ME.

5.2 AOD

In this section, we evaluate the ability of WRF-Chem, CAMS-OA, and MERRA-2 to reproduce the aerosol content in the atmosphere accurately. This content is characterized by AOD. In the ME, mineral dust contribution to the total AOD is dominant (≈87 %; Kalenderski and Stenchikov, 2016; Osipov et al., 2015). The treatment of optically active dust within the model is therefore vitally important. AOD is calculated based on aerosol concentrations and aerosol optical properties, which depend upon aerosol size distribution. We therefore evaluate how well WRF-Chem and assimilation products reproduce aerosol volume size distribution.

5.2.1 Aerosol volume size distributions

Dust particles are emitted into the lower atmospheric layer with some predominant size distribution (Martin and Kok, 2017; Kok, 2011). Emitted dust is processed by the atmosphere to produce the atmospheric dust size distribution that is retrieved by the AERONET inversion algorithm (Dubovik and King, 2000) and reported as column-integrated AVSD. Strictly speaking, AERONET AVSD incorporates contributions from all types of aerosols. But the size distribution of emitted dust has the strongest effect on column-integrated AVSD because dust dominates all other aerosols in the ME. Therefore, we have to tune the dust emission parameters in the first place.

Equation (2) assumes that emission mass fluxes into five dust size bins are controlled by the s_p fractions. In WRF-Chem the default values of s_p fractions for the five dust bins (see Table 2) are 0.1, 0.25, 0.25, 0.25, and 0.25. We found that with these default s_p fractions, WRF-Chem underestimated the volume of fine dust particles in the first bin, $\mathrm{0.1}\mathrm{\mu }\mathrm{m}<r<\mathrm{1}\mathrm{\mu }\mathrm{m}$, compared with AERONET AVSD, whereas the volume of the second bin, $\mathrm{1}\mathrm{\mu }\mathrm{m}<r<\mathrm{1.8}\mathrm{\mu }\mathrm{m}$, was overestimated. In combination with fitting the observed AOD by tuning of factor C, this led to an increase in the total emitted dust mass since fine particles are optically more efficient per unit mass than coarse particles. To achieve a better agreement between the simulated and AERONET AVSDs, we adjusted the fractions s_p to be 0.15, 0.1, 0.25, 0.4, and 0.1. A similar approach was implemented in Khan et al. (2015) using the Modal Aerosol Dynamics Model for Europe (MADE) and Secondary Organic Aerosol Model (SORGAM) chemistry–aerosol scheme. This s_p modification is in line with (Adebiyi and Kok, 2020) as it effectively decreased emission of dust particles with radii r<2.5 µm and increased emission of coarse particles with radii r>2.5 µm (see Appendix A3). We use the updated s_p values in all our WRF-Chem simulations.

Figure 4 Seasonally averaged 2015–2016 AVSDs (µm³ µm⁻²) obtained from MERRA-2, CAMS-OA, WRF-Chem, and the AERONET inversion algorithm at the (a) KAUST Campus, (b) Mezaira, and (c) Sede Boker AERONET sites. Winter (DJF), spring (MAM), summer (JJA), and autumn (SON).

Figure 4 shows seasonally averaged 2015–2016 volume size distributions obtained from MERRA-2, CAMS-OA, AERONET, and WRF-Chem with updated s_p fractions. The comparison is conducted for the KAUST Campus, Mezaira, and Sede Boker AERONET sites (see Fig. 1) since only these sites have information on AVSDs during the 2015–2016 period. The effect of s_p modification could be seen in Fig. 4 by comparing AVSDs from WRF-Chem with an updated set of s_p and MERRA-2 that uses the default s_p set. A direct comparison of AVSDs from the WRF-Chem runs with the updated and default s_p sets is shown in Appendix A3.

Both MERRA-2 and WRF-Chem use the GOCART aerosol scheme with the same five dust bins, and they approximate the shape of the AERONET AVSD relatively well. CAMS-OA uses only three dust bins (see Table 2) and fails to reproduce the AERONET AVSD even qualitatively. It overestimates the volume of particles with radii of 0.55–0.9 µm and underestimates the volume of particles with radii of 0.9–20 µm. With the latest system upgrade in 2019, this weakness of CAMS-OA has been corrected by introducing a new dust scheme (Nabat et al., 2012).

The volume size distributions from the model and assimilation products demonstrate pronounced seasonal variability with the increased amount of dust in the atmosphere during spring and summer. Since the KAUST Campus and Mezaira sites are located in the vicinity of strong dust sources, the coarse mode at these sites is more pronounced than at the Sede Boker site, which is farther from the strong dust emission sources.

The fine mode in the AERONET AVSD is more pronounced at the KAUST Campus site in comparison with the other AERONET sites due to its proximity to strong SO₂ sources located along the west coast of Saudi Arabia (Ukhov et al., 2020b). This proximity leads to a higher contribution of fine sulfate particles to the fine mode. The smaller volume of fine particles in the WRF-Chem and MERRA-2 simulated AVSD (see Fig. 4) is in part because the simulated AVSDs show only dust omitting the contributions of sulfate and sea salt. Sea salt particles and droplets are relatively large and mostly contribute to the coarse mode.

Figure 5 shows the contributions of dust, sea salt, and sulfate aerosols into the AVSD at the KAUST Campus AERONET site in the WRF-Chem simulation averaged for two summer seasons (JJA) in the 2015–2016 period. In WRF-Chem, sulfate aerosol is computed using a bulk approach. For calculating aerosol optical properties, it is assumed that sulfate aerosol comprises two log-normal modes: nuclei and accumulation. According to the WRF-Chem source code, the nuclei mode median radius μ_nuc is 0.005 µm, and the geometric width σ_nuc is 1.7; the accumulation mode median radius μ_acc is 0.035 µm, and the geometric width σ_acc is 2.0. The nuclei mode comprises 25 % of the sulfate aerosol mass and accumulation mode −75 %. It is assumed that sulfate aerosol density is 1800 kg m⁻³, and sea salt density is 2200 kg m⁻³. Figure 5 demonstrates that the contribution of the sulfate nuclei mode in the aerosol volume is almost negligible, while the sulfate accumulation mode adds in the volume of aerosol particles with radii <1 µm. The contribution of the sea salt aerosol into AVSD in WRF-Chem simulations is very little.

Figure 5 Summer (JJA) averaged 2015–2016 AVSD (µm³ µm⁻²) at the KAUST Campus AERONET site obtained from the AERONET inversion algorithm and from WRF-Chem.

5.2.2 Comparison with AERONET AOD

The comparison of the daily averaged AOD time series and corresponding scatter plots calculated using WRF-Chem, MERRA-2, CAMS-OA, MODIS-DB&DT, and MAIAC data with AERONET AOD observations conducted at KAUST Campus, Mezaira and Sede Boker during the 2015–2016 period is presented in Fig. 6. Because AERONET conducts observations only during the daylight, we interpolated WRF-Chem, MERRA-2, and CAMS-OA AODs to the AERONET measurement times and then conducted time averaging to make simulated and observed AODs consistent. AODs from MODIS-DB&DT and MAIAC are provided as a daily average. Although MODIS routinely provides observations only twice a day during daylight, up to four observations might be collected on some days due to overlap of the TERRA and AQUA orbits at some locations.

Figure 6 Daily averaged AOD at three AERONET sites (KAUST Campus, Mezaira, Sede Boker) and corresponding scatter plots computed for WRF-Chem, AERONET, MERRA-2, CAMS-OA, MODIS-DB&DT, and MAIAC: (a) 2015, (b) 2016.

Table 4 Pearson correlation coefficient R and mean bias calculated for daily averaged AOD time series from WRF-Chem, CAMS-OA, MERRA-2, MODIS-DB&DT, and MAIAC with respect to AERONET AOD observations.

The scatter plots show that the model and assimilation products are capable of reproducing the magnitude and temporal evolution of the observed AERONET AOD at all sites. During both years, KAUST Campus and Mezaira sites show higher AOD in summer and lower AOD in winter. To quantify the capability of the WRF-Chem, MERRA-2, and CAMS-OA models as well as the MODIS-DB&DT and MAIAC products to reproduce the AERONET AOD, we calculate the Pearson correlation coefficient R and mean bias (see Appendix A2) with respect to the AERONET AOD observations for the 2015–2016 period (see Table 4). The correlation coefficients are the highest for MERRA-2 and MAIAC. MAIAC shows better correlation than MERRA-2 during 2015 (0.88–0.96), but MERRA-2 is better correlated with AERONET (0.85–0.91) than MAIAC in 2016. CAMS-OA, despite the fact that it does not assimilate AERONET, shows better correlations (0.65–0.87) than MODIS-DB&DT (0.56–0.84). However, CAMS-OA overestimates AOD, particularly during acute dust events, and has a relatively high positive mean bias. The R coefficient for the WRF-Chem AOD is 0.43–0.85. MERRA-2 and WRF-Chem have the lowest mean bias in comparison with the other models and products. MODIS-DB&DT mean bias is positive in 2015 and 2016, while MAIAC mean bias is negative for KAUST Campus and Mezaira and positive for Sede Boker during both years.

We have to mention here that the satellite retrievals and MERRA-2 use AERONET observations for calibration. WRF-Chem is tuned to reduce the annual mean bias with respect to AERONET observations. CAMS-OA does not assimilate AERONET AODs. In WRF-Chem, we did not tune the temporal correlation between the model and AERONET AOD. In this sense, the correlation coefficient between WRF-Chem and AERONET AOD provides an independent evaluation of the model performance (see Table 4). It is expected that the temporal correlation for the assimilation products and satellite retrievals will be higher than for the free-running WRF-Chem.

5.2.3 Comparison of spatial AOD distributions

We also examine how well MERRA-2, CAMS-OA, MAIAC, and WRF-Chem reproduce spatial patterns and seasonal variability of the AOD in comparison with the conventional MODIS-DB&DT retrievals. The seasonally and annually averaged 2015–2016 AOD fields from WRF-Chem, CAMS-OA, MERRA-2, and the two MODIS retrievals DB&DT and MAIAC are presented in Fig. 7. The seasonally averaged AODs from WRF-Chem, MERRA-2, and CAMS-OA are shown at their original spatial resolution and were calculated using only daytime (06:00–14:00 UTC or 09:00–17:00 LT) output. The AODs were sampled under all-sky conditions, which in the ME does not make much of a difference as cloud fraction is low. For statistical comparison, we interpolated AOD fields (preserving the area's average AODs) on the MERRA-2 grid and calculated the Pearson correlation coefficient R root mean square error RMSE and mean bias with respect to MODIS-DB&DT AOD (see Table 5). When conducting statistical analysis, the grid cells with undefined pixels in MODIS-DB&DT and MAIAC retrievals were excluded.

Figure 7 Seasonally averaged 2015–2016 AOD. The right column is annual mean AOD. Rows: (a) WRF-Chem, (b) MERRA-2, (c) CAMS-OA, (d) MODIS-DB&DT, and (e) MAIAC. Winter (DJF), spring (MAM), summer (JJA), and autumn (SON). White dots are undefined pixels. Black triangles denote locations of the KAUST Campus, Mezaira, and Sede Boker AERONET stations.

Table 5 Pearson correlation coefficient (R), root mean square error (RMSE), and mean bias calculated for seasonally and annually averaged 2015–2016 AOD geographic distributions from CAMS-OA, MAIAC, MERRA-2, and WRF-Chem with respect to MODIS-DB&DT AOD.

The statistical scores provided in Table 5 show that the annual mean AOD from MAIAC has the highest correlation (R=0.796) but also the highest RMSE (0.123) and the biggest bias (−0.095) with respect to MODIS-DB&DT AOD. MERRA-2 annual mean AOD has R=0.663 with respect to MODIS-DB&DT AOD, CAMS-OA has R=0.650), and WRF-Chem has R=0.609 with RMSEs of 0.116 for all of them. WRF-Chem, MERRA-2, and CAMS-OA demonstrate similar AOD patterns, but WRF-Chem and MERRA-2 underestimate and CAMS-OA overestimates MODIS-DB&DT AOD during all seasons, with annual mean biases of $=-\mathrm{0.009}$, −0.042, and 0.039, respectively. MAIAC underestimates AOD in comparison with MODIS-DB&DT, which is consistent with the MAIAC and MODIS-DB&DT AOD comparison with AERONET AOD (see Table 4 and Fig. 6).

Based on the comparison of WRF-Chem AOD with the AOD from MODIS and AERONET observations, we conclude that spatial and temporal WRF-Chem's AOD distribution is in good agreement with the available satellite and ground-based observations, i.e., annual mean correlation R exceeds 0.6 (see Table 5), and correlation with AERONET is 0.43–0.85 (see Table 4).

5.3 PM air pollution

To test the ability of the data assimilation products and models to characterize PM air pollution in the ME, we compare surface daily mean PM_2.5 and PM₁₀ concentrations from WRF-Chem, MERRA-2, and CAMS-OA with daily averaged measurements conducted by the three AQMSs (see Figs. 8 and 9). The AQMSs are installed in Jeddah, Riyadh, and Dammam (Fig. 1), the Saudi Arabian megacities. PM measurements conducted by MODON (see Sect. 2.3) are available starting from 2016. The modeled PM_2.5 and PM₁₀ concentrations were sampled from the model fields at the exact AQMS locations. The following formulas were used to calculate PM_2.5 and PM₁₀ surface concentrations using WRF-Chem and MERRA-2 output:

$$\begin{array}{}\text{(4)}& \begin{array}{rl}{\mathrm{PM}}_{\mathrm{2.5}}& ={\mathrm{DUST}}_{\mathrm{1}}+{\mathrm{DUST}}_{\mathrm{2}}\times \mathrm{0.38}+{\mathrm{SEAS}}_{\mathrm{1}}+{\mathrm{SEAS}}_{\mathrm{2}}\\ & +{\mathrm{SEAS}}_{\mathrm{3}}\times \mathrm{0.83}+\mathrm{sulfate}\\ & +\left({\mathrm{OC}}_{\mathrm{1}}+{\mathrm{OC}}_{\mathrm{2}}\right)\times {\mathrm{OC}}_{\mathrm{mfac}}+{\mathrm{BC}}_{\mathrm{1}}+{\mathrm{BC}}_{\mathrm{2}}\\ {\mathrm{PM}}_{\mathrm{10}}& ={\mathrm{DUST}}_{\mathrm{1}}+{\mathrm{DUST}}_{\mathrm{2}}+{\mathrm{DUST}}_{\mathrm{3}}+{\mathrm{DUST}}_{\mathrm{4}}\times \mathrm{0.74}\\ & +{\mathrm{SEAS}}_{\mathrm{1}}+{\mathrm{SEAS}}_{\mathrm{2}}+{\mathrm{SEAS}}_{\mathrm{3}}+{\mathrm{SEAS}}_{\mathrm{4}}+\mathrm{sulfate}\\ & +\left({\mathrm{OC}}_{\mathrm{1}}+{\mathrm{OC}}_{\mathrm{2}}\right)\times {\mathrm{OC}}_{\mathrm{mfac}}+{\mathrm{BC}}_{\mathrm{1}}+{\mathrm{BC}}_{\mathrm{2}},\end{array}\end{array}$$

where ${\mathrm{DUST}}_{\mathrm{1},\mathrm{2},\mathrm{3},\mathrm{4}}$, ${\mathrm{SEAS}}_{\mathrm{1},\mathrm{2},\mathrm{3},\mathrm{4}}$, OC_1,2, BC_1,2, and sulfate are respectively the concentrations of the dust and sea salt in the first four bins, organic and black carbon (hydrophobic and hydrophilic), and sulfate ion (${\mathrm{SO}}_{\mathrm{4}}^{\mathrm{2}-}$). As was mentioned in Sect. 4.1, SEAS₁ is not present in the WRF-Chem output. So for WRF-Chem we assume SEAS₁=0. The factor OC_mfac=1.8 accounts for the conversion of organic carbon into organic matter.

Figure 8 Left: WRF-Chem daily averaged PM_2.5 surface concentrations (µg m⁻³) with MODON observations, MERRA-2, and CAMS-OA at Jeddah, Riyadh, and Dammam. The dash-dotted line corresponds to the 25 µg m⁻³ WHO daily average guideline. Right: stacked bars show the decomposition of the PM_2.5 annual mean surface concentrations into dust and nondust components. The dash-dotted line corresponds to the 10 µg m⁻³ WHO annual guideline. Numbers on the right-hand side of the WRF-Chem, CAMS-OA, and MERRA-2 bars show the contribution (%) of the dust and nondust components to the total PM_2.5 concentration. (a) 2015, (b) 2016.

Figure 9 Left: WRF-Chem daily averaged PM₁₀ surface concentrations (µg m⁻³) with the MODON observations as well as MERRA-2 and CAMS-OA at Jeddah, Riyadh, and Dammam. The dash-dotted line corresponds to the 50 µg m⁻³ WHO daily guideline. Right: stacked bars show the decomposition of the PM₁₀ annual mean surface concentrations into dust and nondust components. The dash-dotted line corresponds to the 20 µg m⁻³ WHO annual guideline. Numbers on the right-hand side of the WRF-Chem, CAMS-OA, and MERRA-2 bars show the contribution (%) of the dust and nondust components to the total PM₁₀ concentration. (a) 2015, (b) 2016.

CAMS-OA PM_2.5 and PM₁₀ were calculated using the following relations (https://confluence.ecmwf.int/display/CUSF/PM10+and+PM25+global+products, last access: 5 January 2020):

$$\begin{array}{}\text{(5)}& \begin{array}{rl}{\mathrm{PM}}_{\mathrm{2.5}}& ={\mathrm{DD}}_{\mathrm{1}}+{\mathrm{DD}}_{\mathrm{2}}+{\mathrm{SS}}_{\mathrm{1}}/\mathrm{4.3}+\mathrm{0.5}\times {\mathrm{SS}}_{\mathrm{2}}/\mathrm{4.3}\\ & +\mathrm{0.7}\times \left({\mathrm{OM}}_{\mathrm{1}}+{\mathrm{OM}}_{\mathrm{2}}+\mathrm{sulfate}\right)+{\mathrm{BC}}_{\mathrm{1}}+{\mathrm{BC}}_{\mathrm{2}}\\ {\mathrm{PM}}_{\mathrm{10}}& ={\mathrm{DD}}_{\mathrm{1}}+{\mathrm{DD}}_{\mathrm{2}}+{\mathrm{DD}}_{\mathrm{3}}\times \mathrm{0.4}+{\mathrm{SS}}_{\mathrm{1}}/\mathrm{4.3}\\ & +{\mathrm{SS}}_{\mathrm{2}}/\mathrm{4.3}+{\mathrm{OM}}_{\mathrm{1}}+{\mathrm{OM}}_{\mathrm{2}}+\mathrm{sulfate}+{\mathrm{BC}}_{\mathrm{1}}\\ & +{\mathrm{BC}}_{\mathrm{2}},\end{array}\end{array}$$

where ${\mathrm{DD}}_{\mathrm{1},\mathrm{2},\mathrm{3}}$, SS_1,2, sulfate, BC_1,2, and OM_1,2 are surface concentrations of dust in three bins, sea salt in two bins, sulfate, black carbon, and organic matter (hydrophobic and hydrophilic). The size ranges of dust and sea salt bins from CAMS-OA are presented in Table 2.

Table 6 Contributions (%) of dust and nondust components to PM_2.5 for Jeddah, Riyadh, and Dammam during 2015–2016.

¹ For WRF-Chem and MERRA-2: $\mathrm{dust}={\mathrm{DUST}}_{\mathrm{1}}+{\mathrm{DUST}}_{\mathrm{2}}\times \mathrm{0.38}$, $\mathrm{BC}={\mathrm{BC}}_{\mathrm{1}}+{\mathrm{BC}}_{\mathrm{2}}$, sulf=sulfate, $\mathrm{OM}=({\mathrm{OC}}_{\mathrm{1}}+{\mathrm{OC}}_{\mathrm{2}})\times {\mathrm{OC}}_{\mathrm{mfac}}$, $\mathrm{salt}={\mathrm{SS}}_{\mathrm{1}}+{\mathrm{SS}}_{\mathrm{2}}+{\mathrm{SS}}_{\mathrm{3}}\times \mathrm{0.83}$. ² For
CAMS-OA: $\mathrm{Dust}={\mathrm{DD}}_{\mathrm{1}}+{\mathrm{DD}}_{\mathrm{2}}$, $\mathrm{Sulf}=\mathrm{0.7}\times \mathrm{sulfate}$, $\mathrm{BC}={\mathrm{BC}}_{\mathrm{1}}+{\mathrm{BC}}_{\mathrm{2}}$, $\mathrm{OM}=\mathrm{0.7}\times ({\mathrm{OM}}_{\mathrm{1}}+{\mathrm{OM}}_{\mathrm{2}}$), $\mathrm{Salt}={\mathrm{SS}}_{\mathrm{1}}/\mathrm{4.3}+\mathrm{0.5}\times {\mathrm{SS}}_{\mathrm{2}}/\mathrm{4.3}$.
Abbreviations of the aerosols' names correspond to those given in Sect. 5.3.

Table 7 Contributions (%) of dust and nondust components to PM₁₀ for Jeddah, Riyadh, and Dammam during 2015–2016.

¹ For WRF-Chem and MERRA-2: $\mathrm{dust}={\mathrm{DUST}}_{\mathrm{1}}+{\mathrm{DUST}}_{\mathrm{2}}+{\mathrm{DUST}}_{\mathrm{3}}+{\mathrm{DUST}}_{\mathrm{4}}\times \mathrm{0.74}$, sulf=sulfate, $\mathrm{BC}={\mathrm{BC}}_{\mathrm{1}}+{\mathrm{BC}}_{\mathrm{2}}$, $\mathrm{OM}=({\mathrm{OC}}_{\mathrm{1}}+{\mathrm{OC}}_{\mathrm{2}})\times {\mathrm{OC}}_{\mathrm{mfac}}$, $\mathrm{salt}={\mathrm{SS}}_{\mathrm{1}}+{\mathrm{SS}}_{\mathrm{2}}+{\mathrm{SS}}_{\mathrm{3}}$
+SS₄. ² For CAMS-OA: $\mathrm{Dust}={\mathrm{DD}}_{\mathrm{1}}+{\mathrm{DD}}_{\mathrm{2}}+{\mathrm{DD}}_{\mathrm{3}}\times \mathrm{0.4}$, Sulf=sulfate, $\mathrm{BC}={\mathrm{BC}}_{\mathrm{1}}+{\mathrm{BC}}_{\mathrm{2}}$, $\mathrm{OM}={\mathrm{OM}}_{\mathrm{1}}+{\mathrm{OM}}_{\mathrm{2}}$, $\mathrm{Salt}={\mathrm{SS}}_{\mathrm{1}}/\mathrm{4.3}+{\mathrm{SS}}_{\mathrm{2}}/\mathrm{4.3}$.
Abbreviations of the aerosols' names correspond to those given in Sect. 5.3.

The histograms in the right-side panels in Figs. 8 and 9 show the annual mean PM concentrations from WRF-Chem, MERRA-2, and CAMS-OA split into the dust and nondust components. The dashed and dash-dotted horizontal lines correspond to KSA-PME limits and WHO air quality guidelines for daily (in the left-side panels) and annual mean (in the right-side panels) PM concentrations. We also calculated the separate contributions of sulfate, sea salt, organic matter, and black carbon to the nondust PM_2.5 and PM₁₀ (see Tables 6 and 7, respectively).

The sporadic peaks in the observations, which are not captured by the model and assimilation products, are due to unaccounted factors, such as nearby traffic, construction work, and local anthropogenic or natural emissions, which are not present in the emission inventories, or they are due to meteorological fluctuations that are not resolved in the models. Talking about extreme dust pollution cases, we analyzed dust surface concentrations using WRF-Chem output during the dust storm that took place in the Jeddah region on 8 July 2016. The calculated surface concentrations in all dust bins ${\mathrm{DUST}}_{\mathrm{1},\mathrm{2},\mathrm{3},\mathrm{4},\mathrm{5}}$ at the peak of the storm were 55, 58, 63, 111, and 11 µg m⁻³, respectively. The sum of all dust bins yields a total dust concentration of 298 µg m⁻³.

5.3.1 PM_2.5

Figure 8 shows that the daily averaged PM_2.5 concentrations observed by MODON AQMSs at all locations never drop below the WHO limit of 25 µg m⁻³. During the severe dust events, this limit is exceeded in 2016 by a factor of 10–15. The less restrictive KSA-PME limit of 35 µg m⁻³ is exceeded by a factor of 7–11 during the dust outbreaks. Annually averaged MODON measurements are 8–18 times higher than the 10 µg m⁻³ WHO limit and 5–12 times higher than the 15 µg m⁻³ KSA-PME limit for annual mean PM_2.5 concentrations.

Both data assimilation products and WRF-Chem underestimate annual mean PM_2.5 concentrations by a factor of ≈3 in Jeddah and Riyadh and slightly overestimate (though WRF-Chem slightly underestimates) PM_2.5 in Dammam in comparison with observed concentrations during 2016. The CAMS-OA annual mean surface PM_2.5 concentrations in Jeddah and Riyadh are higher than those from WRF-Chem and MERRA-2, providing the best fit for MODON observations, at least on an annual mean (during 2016) basis.

Annual mean PM_2.5 concentrations from WRF-Chem and MERRA-2 exceed the WHO limit of 10 µg m⁻³ by a factor of ≈4–7 and ≈6–10, respectively, in all locations. The KSA-PME limit of 15 µg m⁻³ for annual average PM_2.5 concentrations is exceeded by a factor of ≈2.5–4.5 and ≈4–6.5, respectively, for WRF-Chem and MERRA-2.

In Jeddah and Dammam, WRF-Chem and MERRA-2 show similar relative contributions of nondust components to PM_2.5 (30 %–34 % in Jeddah and 12 %–14 % in Dammam), but in MERRA-2 sea salt is a major contributor to nondust PM_2.5, while in WRF-Chem it is sulfate (see Table 6). This difference between WRF-Chem and MERRA-2 is mainly because MERRA-2 generates more sea salt but also because MERRA-2 underestimates SO₂ emissions located in the Arabian Gulf and along the west coast of Saudi Arabia (Ukhov et al., 2020b) and hence underestimates sulfate concentrations, as discussed in Sect. 4.1. In Riyadh, the contribution of the nondust component to PM_2.5 is ≈9 %–12 % for both MERRA-2 and WRF-Chem. In CAMS-OA, the contribution of nondust particulates to PM_2.5 in Jeddah and Dammam is ≈7 %–10 %, and the contribution of sea salt is little. According to Table 6, in all considered cities, the contribution of black carbon (BC) to PM_2.5 is not significant for all models. In MERRA-2, the contribution of organic matter (OM) to PM_2.5 is more substantial (but still minor) in comparison with WRF-Chem and CAMS-OA. In general, among all models the contribution of dust to PM_2.5 in Jeddah is 65 %–90 %, while in Riyadh and Dammam this contribution is 85 %–95 % (see Table 6).

5.3.2 PM₁₀

Daily averaged MODON measurements almost continuously exceed the WHO guideline of 50 µg m⁻³ at all locations (see Fig. 9). In Riyadh and Dammam, PM₁₀ concentration is higher than in Jeddah, where the KSA-PME limit of 340 µg m⁻³ for daily averaged PM₁₀ is exceeded in 2016 by a factor of about a dozen. In Dammam, this limit is more frequently exceeded, especially during the summer period. During acute dust events in Dammam, daily averaged PM₁₀ concentrations can exceed the WHO guideline limit by a factor of more than 10–20. Annually averaged MODON measurements are 7–11 times higher than the 20 µg m⁻³ WHO guideline and 2–3 times higher than the 80 µg m⁻³ KSA-PME limits for annual mean PM₁₀ concentrations.

In contrast with MERRA-2 and CAMS-OA, WRF-Chem compares better with PM₁₀ observations by MODON in all locations. MERRA-2 overestimates by a factor of ≈1.2–1.8 and CAMS-OA underestimates by a factor of ≈1.5–2 annual mean PM₁₀ MODON observations in all locations. This is in agreement with Cuevas et al. (2015), who stated that MACC (the predecessor of CAMS-OA) underestimates PM₁₀ daily and monthly means all year long, and with our findings in Sect. 5.2.1, where we have shown that CAMS-OA underestimates the volume of particles with radii of 0.9–20 µm. Annual mean PM₁₀ concentrations from WRF-Chem and MERRA-2 exceed the WHO limit of 20 µg m⁻³ by factors of ≈6–15 and ≈10–20, respectively, in all locations. The KSA-PME limit of 80 µg m⁻³ for annual average PM₁₀ concentrations is exceeded by factors of ≈1.5–4 and ≈1.5–5, respectively, for WRF-Chem and MERRA-2.

According to Table 7, MERRA-2 shows the highest contribution of sea salt to PM₁₀ in the coastal cities of Jeddah (≈21 %–25 %) and Dammam (≈6 %–7 %). MERRA-2 demonstrates the lowest (≈1 %–2 %) contribution of sulfate to PM₁₀, while WRF-Chem and CAMS-OA contribute a similar amount of sulfate to PM₁₀ in Jeddah (≈7 %) and Riyadh (≈2.4 %). MERRA-2 also shows the lowest contribution (≈0.1 %–0.2 %) of black carbon (BC) to PM₁₀ in all considered cities. The CAMS-OA organic matter (OM) contribution to PM₁₀ is in 2–8 times greater than the WRF-Chem and MERRA-2 contributions. CAMS-OA demonstrates the lowest (0.1 %–0.2 %) contribution of sea salt to PM₁₀. The contribution of dust to PM₁₀ in Jeddah is 70 %–90 %, while in Riyadh and Dammam this contribution is 90 %–96 %. Minimal contribution (≈3.5 %–4 %) of nondust components to PM₁₀ is observed among all models in Riyadh.

5.3.3 Spatial patterns of PM air pollution

In this section we use the WRF-Chem output averaged for 2015–2016 to discuss the spatial patterns of aerosol pollution and partial contributions from natural and anthropogenic aerosols to PM over the ME.

Figure 10 WRF-Chem's average 2015–2016 surface concentrations (µg m⁻³): (a) PM_2.5; (b) PM₁₀; (c) PM_2.5 ∕ PM₁₀ ratio; (d) black carbon and organic matter, $({\mathrm{OC}}_{\mathrm{1}}+{\mathrm{OC}}_{\mathrm{2}})\times {\mathrm{OC}}_{\mathrm{mfac}}+{\mathrm{BC}}_{\mathrm{1}}+{\mathrm{BC}}_{\mathrm{2}}$, (µg m⁻³); (e) sea salt, ${\mathrm{SEAS}}_{\mathrm{2}}+{\mathrm{SEAS}}_{\mathrm{3}}+{\mathrm{SEAS}}_{\mathrm{4}}$, (µg m⁻³); (f) sulfate (µg m⁻³) and locations of AERONET stations; (g) dust PM_2.5 ∕ (total PM_2.5) ratio; (h) dust PM₁₀ ∕ (total PM₁₀) ratio; (i) sulfate ∕ (PM_2.5 total nondust) ratio. Abbreviations of the aerosols' names correspond to those given in Sect. 5.3.

Figure 10a, b, c show the spatial distributions of the PM_2.5 and PM₁₀ surface concentrations and the PM_2.5 ∕ PM₁₀ ratio. The left limits of the color bars for PM_2.5 and PM₁₀ are set to be equal to the corresponding WHO annual guideline concentrations. Over the whole domain, the annual mean surface concentrations of PM_2.5 and PM₁₀ exceed WHO guidelines of 10 and 20 µg m⁻³, respectively. The regions of high surface concentrations coincide with the main dust sources, which span from northern Iraq to Oman, and include Sudan, Egypt, Libya, and Turkmenistan. PM concentrations in these regions exceed even the less restrictive KSA-PME air quality limit for annual mean PM_2.5 and PM₁₀ concentrations by a factor of more than 5.

In the entire domain, the maximum, minimum, and mean values of the PM_2.5 ∕ PM₁₀ ratio (see Fig. 10c) are 0.73, 0.20, and 0.31, respectively. As expected, lower PM_2.5 ∕ PM₁₀ ratios (0.2–0.3) are observed over the dust source regions (i.e., along the eastern Arabian Peninsula, Iraq, and northern Africa), where both coarse and fine particles are generated. However, large particles cannot be transported as far from source regions as small particles due to the shorter lifetime of large particles compared with small particles with respect to deposition processes. The higher values (0.4–0.6) of the PM_2.5 ∕ PM₁₀ ratio are observed over southeastern Europe, Turkey, Ethiopia, and western parts of the Arabian Peninsula that are farther from the main dust sources.

Figure 10d shows the sum of surface concentrations of organic matter and black carbon, $({\mathrm{OC}}_{\mathrm{1}}+{\mathrm{OC}}_{\mathrm{2}})\times {\mathrm{OC}}_{\mathrm{mfac}}+{\mathrm{BC}}_{\mathrm{1}}+{\mathrm{BC}}_{\mathrm{2}}$). Their maximum, minimum, and mean concentration values are 31.8, 0.2, and 1.3 µg m⁻³, respectively. Their contribution to aerosol pollution over the Arabian Peninsula in WRF-Chem simulations is insignificant. Figure 10e shows the surface concentration of sea salt calculated as a sum of concentrations in each bin ${\mathrm{SEAS}}_{\mathrm{2}}+{\mathrm{SEAS}}_{\mathrm{3}}+{\mathrm{SEAS}}_{\mathrm{4}}$. Over the seas and coastal areas, the average concentration of sea salt is 3–12 µg m⁻³. In summer, strong winds in Somali jets (see Fig. 3b) intensify sea salt emission over the Arabian Sea, creating high surface concentrations of sea salt (27–42 µg m⁻³) along the coasts of Somalia and Oman. Due to prevailing northern winds, the transport of sea salt from the Mediterranean Sea to Egypt and Libya is observed.

The relatively high sulfate surface concentration (see Fig. 10f) is observed in the vicinity of the strong SO₂ sources located along the west and east coast of Saudi Arabia and over the Arabian Gulf as well as downwind from those sources (see Ukhov et al., 2020b, for details). Figure 10f also denotes the locations of the AERONET stations, as in Fig. 1. The sulfate concentration at the KAUST Campus site is higher than at the Mezaira and Sede Boker AERONET sites (see Sect. 5.2.1), so it experiences a more pronounced contribution of sulfate particulates to the fine mode of the AVSD (see Figs. 4a and 5). Due to the prevailing northern winds along the Red Sea, sulfate aerosols originating from SO₂ emissions along the Red Sea coast spread along the west coast of the Arabian Peninsula towards Yemen. The drift of sulfate from the Arabian Gulf into the interior of the eastern part of the Arabian Peninsula is also seen. The sulfate annual mean background surface concentration over the US for the period 2003–2012, obtained in Buchard et al. (2016), is 2–3 µg m⁻³, similar to the background concentrations we see in the ME. But downwind or in the vicinity of strong SO₂ point emissions, sulfate concentrations are 3–4 times higher. Similar sulfate surface concentrations for the period 2000–2016 over the US were obtained in van Donkelaar et al. (2019), where the concentrations reach ≈10 µg m⁻³ over the eastern part of the US during summer. In Al-Taani et al. (2019) the average 1980–2016 sulfate concentration computed for the UAE is 2.5–3 times lower. This difference is caused by the fact that Al-Taani et al. (2019) took the sulfate fields from MERRA-2 reanalysis, which underestimates the SO₂ emissions as shown in Ukhov et al. (2020b).

Contributions of dust to PM_2.5 and PM₁₀ calculated as ratios of dust PM_2.5 to total PM_2.5 and dust PM₁₀ to total PM₁₀ are shown in Fig. 10g and h, respectively. Due to relatively low dust surface concentrations over the eastern Mediterranean, Turkey, and southeastern Europe, the contribution of dust to PM_2.5 and PM₁₀ is 20 %–50 % and 50 %–80 %, respectively. In other areas that are closer to the dust source regions, the contribution of dust to PM is above 80 %.

Figure 10i shows the ratio between the concentration of sulfate aerosol with respect to the total concentration of PM_2.5 nondust aerosols. The maximum, minimum, and mean values of this ratio are 0.84, 0.07, and 0.52, respectively. This ratio is low over the seas, where sea salt is prevalent but consistently exceeds 0.6 over land. Sulfate, therefore, is the primary anthropogenic pollutant over land. In the central and southern parts of Saudi Arabia and over Iran, sulfate contributes 60 %–90 % to the total PM_2.5 nondust aerosol concentration. Over the other parts of the Arabian Peninsula, the northern part of Sudan, Libya, and Egypt, sulfate contributes approximately 40 %–60 % to the total PM_2.5 nondust aerosol concentration.

5.3.4 PM air pollution in the major ME cities

To evaluate the air quality in the ME's major cities, we calculate for their locations the average 2015–2016 daily PM_2.5 and PM₁₀ surface concentrations, their 90th percentiles, and the contribution of the dust and nondust components to PM (see Fig. 11). We calculate the number of days during the 2015–2016 period when the daily PM_2.5 and PM₁₀ surface concentrations exceed the US EPA air quality limit of 35 and 150 µg m⁻³, respectively.

Figure 11 Annual mean 2015–2016 PM surface concentrations (µg m⁻³) calculated for the major ME cities and PM decomposition into dust and nondust (sulfate, sea salt, black carbon, and organic matter) components (stacked bars). Abbreviations of the aerosols' names correspond to those given in Sect. 5.3. Hatched bars denote 90th percentiles (µg m⁻³) calculated using daily mean PM concentrations. WHO guidelines and KSA-PME air quality limits for annually averaged PM are shown by dash-dotted and dashed lines. Numbers over the stacked bars correspond to the number of days during 2015–2016 when daily averaged PM surface concentration exceeded the US EPA air quality limit. (a) PM_2.5; the daily averaged US EPA air quality limit is 35 µg m⁻³, and the annual WHO guideline and KSA-PME limit are 10 and 15 µg m⁻³, respectively. (b) PM₁₀; the daily averaged US EPA air quality limit is 150 µg m⁻³, and the annual WHO and KSA-PME limits are 20 and 80 µg m⁻³, respectively.

Figure 11 shows that the annually averaged PM_2.5 and PM₁₀ exceed the WHO air quality guidelines by factors of 2–9 and 3–20, respectively in all major cities of the ME except Ankara. The KSA-PME air quality limit for annual mean PM_2.5 is exceeded by a factor of up to 6 and by a factor of up to 5 for PM₁₀. Due to the lack of strong dust sources nearby, air quality conditions in the cities in the eastern Mediterranean are more favorable when compared with those in the Arabian Peninsula. In these cities, the air pollution shifts from natural to anthropogenic as the contribution of nondust aerosols to PM_2.5 increases up to 30 %–45 %, in contrast with the cities located in the Arabian Peninsula, where this contribution is up to 8 %–25 %. Sulfate aerosol is the major contributor to nondust PM_2.5.

The cities at the eastern coast of the Arabian Peninsula have the highest 90th percentiles of daily mean PM concentrations. For example, in Dammam, Abu Dhabi, Doha, and Kuwait City, the 90th percentiles of daily mean surface concentration of PM₁₀ and PM_2.5 are in the range of 400–740 and 130–180 µg m⁻³, respectively. This is above the KSA-PME air quality limits for daily mean PM₁₀ and PM_2.5.

In eastern Mediterranean cities, within the 2015–2016 period, the US EPA air quality daily mean limits are exceeded on 40–75 d for PM₁₀ and 60–100 d for PM_2.5. In the cities of the Arabian Peninsula, Iraq, and Iran, the US EPA PM daily mean limits are exceeded on 94–627 d for PM₁₀ and 213–640 d for PM_2.5 during the same period.

6 Conclusions

This study evaluates MERRA-2 and CAMS-OA data assimilation products and high-resolution WRF-Chem simulations aimed at assessing the impact of aerosols on PM air pollution over the Middle East for the 2015–2016 period. It also compares the new MODIS AOD retrieval, MAIAC, with the conventional MODIS-DB&DT and AERONET AOD over the Middle East's major dust source regions. MERRA-2 and WRF-Chem use the five-bin dust aerosol model and demonstrate a better agreement with the AERONET-retrieved size distribution than CAMS-OA, which uses a three-dust-bin microphysical model. CAMS-OA overestimates the volume of fine dust particles with radii of 0.55–0.9 µm and underestimates the volume of coarse dust particles with radii of 0.9–20 µm in comparison with the AERONET aerosol volume size distribution.

At all considered AERONET sites, WRF-Chem, CAMS-OA, MERRA-2, MODIS-DB&DT, and MAIAC are capable of reproducing the magnitude and temporal evolution of the AERONET AOD time series during the whole period. MAIAC and MERRA-2 have the highest correlation to AERONET AOD. CAMS-OA tends to overestimate AERONET AOD, especially during severe dust events. The MODIS-DB&DT and MAIAC AOD mean biases with respect to AERONET observations are of the same magnitude (slightly larger than that of MERRA-2), but the MODIS-DB&DT AOD is biased positively, and the MAIAC AOD is biased negatively except for Sede Boker for both years. The AOD fields from WRF-Chem and assimilation products exhibit similar spatial patterns, but WRF-Chem, MAIAC, and MERRA-2 underestimate and CAMS-OA overestimates MODIS-DB&DT AOD. MAIAC has the highest spatial correlation to the conventional MODIS-DB&DT AOD, followed by MERRA-2, CAMS-OA, and WRF-Chem.

The capability of WRF-Chem, MERRA-2, and CAMS-OA to reproduce PM air pollution over the Middle East was tested against in situ measurements. These PM measurements are conducted in the industrial regions of Jeddah, Riyadh, and Dammam, which complicates one-to-one comparison with the output from global and regional models. Annual mean PM concentrations from WRF-Chem and MERRA-2 exceed the WHO limit by a factor of almost 20. The KSA-PME limit for annual average concentrations is also exceeded by a factor of more than 6. MERRA-2 and WRF-Chem underestimate the observed annual mean PM_2.5 concentrations during 2016 in Jeddah and Riyadh almost 3 times. CAMS-OA and MERRA-2 overestimate and WRF-Chem underestimates observed annual mean PM_2.5 in Dammam. CAMS-OA underestimates (by a factor of 1.5–2) annual mean PM₁₀ observations in all locations primarily due to its deficient dust size distribution. The CAMS-OA annual mean PM_2.5 fits the PM observations better than other products.

The PM composition analysis over rural areas shows that in WRF-Chem, the annual average PM_2.5/PM₁₀ ratio over the ME is about 0.3. It decreases to 0.25 over the major dust source regions, i.e., in the eastern Arabian Peninsula, Iraq, and northern Africa. In most parts of the Middle East, dust is the major contributor to PM. The sulfate aerosol contribution to PM_2.5 is essential in the areas where strong SO₂ sources are present, i.e., on the west and east coasts of Saudi Arabia and over the Arabian Gulf. In these areas sulfate surface concentration reaches 8–11 µg m⁻³, while the “clean” background level is 2–4 µg m⁻³. High sulfate content along the west coast of Saudi Arabia is consistent with the increased volume of the fine mode at the KAUST Campus AERONET site in comparison with the Mezaira and Sede Boker sites. In WRF-Chem, sulfate is the major nondust pollutant over the Middle East. Sulfate aerosols contribute 60 %–90 % to the total PM_2.5 nondust aerosols over the central and southern parts of Saudi Arabia and Iran. Over the other parts of the Arabian Peninsula, northern Sudan, Libya, and Egypt, sulfate contributes approximately 40 %–60 % to the total PM_2.5 nondust aerosol concentration.

The analysis of the annually averaged PM_2.5 and PM₁₀ surface concentrations in the Middle East's major cities shows a very high PM pollution level. In Dammam, Abu Dhabi, Doha, and Kuwait City, the 90th percentile of PM₁₀ and PM_2.5 annual mean surface concentrations exceeds 400–740 and 130–180 µg m⁻³, respectively, which is above the KSA-PME air quality limit. In the eastern Mediterranean, dust concentration drops, and nondust aerosols' contribution to PM_2.5 increases up to 30 %–45 %. In the cities located on the Arabian Peninsula, the contribution of the nondust component to PM_2.5 is 8 %–25 %, which limits the effect of the emission control on air quality. In the eastern Mediterranean cities, the daily mean surface PM concentrations exceed the US EPA air quality daily mean limit on 40–75 d for PM₁₀ and 60–100 d for PM_2.5 during the 2015–2016 period. In the major cities over the Arabian Peninsula, Iraq, and Iran, the US EPA air quality daily mean limit is exceeded on 94–627 d for PM₁₀ and 213–640 d for PM_2.5. In Jeddah and Dammam, WRF-Chem and MERRA-2 show similar relative contributions of the nondust component to PM_2.5 (30 %–34 % in Jeddah and 12 %–14 % in Dammam). But in MERRA-2, sea salt is a major contributor to nondust aerosol concentration, while in WRF-Chem, it is a sulfate. This difference is because MERRA-2 both generates more sea salt and underestimates SO₂ emissions and, consequently, sulfate concentrations. In CAMS-OA, the contribution of nondust particulates to PM_2.5 in Jeddah and Dammam is ≈7 %–10 %, and the contribution of sea salt is little. In Riyadh, the contribution of the nondust component to PM_2.5 is ≈9 %–12 % for both MERRA-2 and WRF-Chem. In Jeddah, Dammam, and Riyadh, the contribution of black carbon to PM_2.5 is insignificant for all products. MERRA-2 shows the highest contribution of sea salt and the lowest contribution of black carbon and sulfate to PM₁₀ in all locations. CAMS-OA demonstrates the lowest contribution of sea salt to PM₁₀. The minimum contribution of nondust components to PM₁₀ is observed in Riyadh among all models.

Thus, in this study, we found that MERRA-2 and CAMS-OA assimilation products as well as WRF-Chem output, despite some intrinsic uncertainties, could be used for evaluation of the PM air pollution over the ME. All products show the dominant contribution of mineral dust to PM. However, in the Arabian coastal areas, where SO₂ emissions are high, contributions of both sulfate and sea salt could be significant. The broad effect of natural aerosols on air quality in the ME puts stricter requirements on anthropogenic pollution control. The impact of dust could be alleviated by employing architectural solutions that are specific to desert areas, increasing in-city vegetation cover, and providing air quality forecasts to alarm the population of hazardous air quality. The developed WRF-Chem modeling framework can be used to simulate other pollutants like NO_x and O₃. The results of the current research could serve as a basis for an improved air quality forecast system that interactively calculates high-resolution, radiative, dynamical, atmospheric chemistry and aerosol processes, driven by natural and anthropogenic emissions. This system will be especially valuable for the prediction of extreme pollution events. It will also improve understanding of the impact of anthropogenic and natural pollution on air quality and human health in the ME region.

Appendix A

A1 Meteorological boundary and initial conditions

To be consistent with IC&BC for chemical species and aerosols, we developed a procedure to build meteorological IC&BC from MERRA-2 reanalysis for all required by WRF-Chem meteorological parameters. In particular, the following 3D parameters were processed: pressure (Pa), geopotential height (m), temperature (K), meridional and zonal wind components (m s⁻¹), and relative humidity (%). The 2D parameters include the following: surface pressure (Pa); sea level pressure (Pa); meridional and zonal wind components at 10 m (m s⁻¹); temperature at 2 m (K); relative humidity at 2 m (%); skin temperature (K); ice mask (0/1); terrain height (m); land/sea mask (1/0); soil temperature at 0–10 (cm), 10–40 (cm), 40–100 (cm), and 100–200 (cm); soil moisture at 0–10 (cm), 10–40 (cm), 40–100 (cm), and 100–200 (cm); snow depth (m); and snow water equivalent (kg m⁻²).

A2 Statistics

We calculated the following statistical parameters to quantify the level of agreement between estimations and observations.

Pearson correlation coefficient (R):

$$\begin{array}{}\text{(A1)}& R={\displaystyle \frac{\sum _{i=\mathrm{1}}^N\left(F_i-\stackrel{\mathrm{‾}}{F}\right)\left(O_i-\stackrel{\mathrm{‾}}{O}\right)}{\sqrt{\sum _{i=\mathrm{1}}^N{\left(F_i-\stackrel{\mathrm{‾}}{F}\right)}^{\mathrm{2}}\sum _{i=\mathrm{1}}^N{\left(O_i-\stackrel{\mathrm{‾}}{O}\right)}^{\mathrm{2}}}}};\end{array}$$

root mean square error (RMSE):

$$\begin{array}{}\text{(A2)}& \mathrm{RMSE}=\sqrt{{\displaystyle \frac{\mathrm{1}}{N}}\sum _{i=\mathrm{1}}^N{\left(F_i-O_i\right)}^{\mathrm{2}}};\end{array}$$

mean bias:

$$\begin{array}{}\text{(A3)}& \mathrm{bias}={\displaystyle \frac{\mathrm{1}}{N}}\sum _{i=\mathrm{1}}^N\left(F_i-O_i\right);\end{array}$$

F_i is the estimated value, O_i is the observed value, $\stackrel{\mathrm{‾}}{F}=\frac{\mathrm{1}}{N}\sum _{i=\mathrm{1}}^N{F}_i$ and $\stackrel{\mathrm{‾}}{O}=\frac{\mathrm{1}}{N}\sum _{i=\mathrm{1}}^N{O}_i$ are their averages, and N is the number of data.

A3 Comparison of AERONET and WRF-Chem volume size distributions

The GOCART dust emission formula (2) calculates dust mass flux into the atmosphere within five dust bins. In this formula the factor C controls the total dust emission mass flux, and the s_p fractions split this flux into five different dust bins. We assume that ∑s_p is 1. To match the observed AERONET AVSD, we changed the default $s_p=\mathit{\{}\mathrm{0.1},\mathrm{0.25},\mathrm{0.25},\mathrm{0.25},\mathrm{0.25}\mathit{\}}$ to $\mathit{\{}\mathrm{0.15},\mathrm{0.1},\mathrm{0.25},\mathrm{0.4},\mathrm{0.1}\mathit{\}}$. This means that 15 % of the total dust mass flux comes as clay and 85 % as silt.

In the original formulation the fractions s_p are not normalized, and ∑s_p is 1.1. This is not essential as the total flux is multiplied by the factor C that is tuned to fit the observed AOD. So we can normalize the original s_p fractions by dividing them by 1.1 and multiplying factor C by 1.1. This will not change any results in Eq. (2) but gives the s_p set of $\mathit{\{}\mathrm{0.09},\mathrm{0.2275},\mathrm{0.2275},\mathrm{0.2275},\mathrm{0.2275}\mathit{\}}$ that is normalized to 1 consistently with our approach. Figure A1 compares the AVSDs calculated with the updated and default s_p fractions for summer (JJA) of 2015.

Figure A1 Volume size distributions at KAUST Campus, Mezaira and Sede Boker AERONET sites averaged for JJA of 2015 from WRF-Chem (bars) and from AERONET (solid line): (a) updated s_p fractions $=\mathit{\{}\mathrm{0.15},\mathrm{0.1},\mathrm{0.25},\mathrm{0.4},\mathrm{0.1}\mathit{\}}$, (b) default s_p fractions $=\mathit{\{}\mathrm{0.1},\mathrm{0.25},\mathrm{0.25},\mathrm{0.25},\mathrm{0.25}\mathit{\}}$.