ModIs Dust AeroSol (MIDAS): a global fine-resolution dust optical depth data set

Monitoring and describing the spatiotemporal variability in dust aerosols is crucial for understanding their multiple effects, related feedbacks, and impacts within the Earth system. This study describes the development of the ModIs Dust AeroSol (MIDAS) data set. MIDAS provides columnar daily dust optical depth (DOD) at 550 nm at a global scale and fine spatial resolution (0.1^∘ × 0.1^∘) over a 15-year period (2003–2017). This new data set combines quality filtered satellite aerosol optical depth (AOD) retrievals from MODIS-Aqua at swath level (Collection 6.1; Level 2), along with DOD-to-AOD ratios provided by the Modern-Era Retrospective analysis for Research and Applications version 2 (MERRA-2) reanalysis to derive DOD on the MODIS native grid. The uncertainties of the MODIS AOD and MERRA-2 dust fraction, with respect to the AEronet RObotic NETwork (AERONET) and LIdar climatology of vertical Aerosol Structure for space-based lidar simulation (LIVAS), respectively, are taken into account for the estimation of the total DOD uncertainty. MERRA-2 dust fractions are in very good agreement with those of LIVAS across the dust belt in the tropical Atlantic Ocean and the Arabian Sea; the agreement degrades in North America and the Southern Hemisphere, where dust sources are smaller. MIDAS, MERRA-2, and LIVAS DODs strongly agree when it comes to annual and seasonal spatial patterns, with colocated global DOD averages of 0.033, 0.031, and 0.029, respectively; however, deviations in dust loading are evident and regionally dependent. Overall, MIDAS is well correlated with AERONET-derived DODs (R=0.89) and only shows a small positive bias (0.004 or 2.7 %). Among the major dust areas of the planet, the highest R values (>0.9) are found at sites of North Africa, the Middle East, and Asia. MIDAS expands, complements, and upgrades the existing observational capabilities of dust aerosols, and it is suitable for dust climatological studies, model evaluation, and data assimilation.

1 Introduction

Among tropospheric and stratospheric aerosol species, dust aerosol is the most abundant component in terms of mass, contributing more than half of the global aerosol amount (Textor et al., 2006; Zender et al., 2011). Preferential sources of dust aerosols are located in areas where precipitation is low, thus favoring aridity, whereas a significant contributing factor is the accumulation of alluvial sediments. Such regions comprise deserts, dry lake beds, and ephemeral channels (e.g., Middleton and Goudie, 2001; Prospero et al., 2002; Ginoux et al., 2012). Previous studies (Prospero et al., 2002; Ginoux et al., 2012) have shown that the major portion of the global dust burden originates from the Sahara, which hosts the most intense dust source of the planet, the Bodélé Depression located in the northern Lake Chad Basin. In North Africa, large amounts of mineral particles are also emitted in the western Sahara, while other noticeable sources of smaller spatial extension are located in the eastern Libyan Desert, in the Nubian Desert (Egypt), and Sudan (Engelstaedter et al., 2006).

One of the major dust sources of the planet, following North Africa, is the Middle East, with several active regions (Pease et al., 1998; Hamidi et al., 2013; Yu et al., 2013) in which wind-blown dust is emitted from alluvial plains (Tigris–Euphrates river system) and sandy deserts (Rub' al Khali desert). Important dust sources are also recorded in the Asian continent, particularly in the Taklamakan Desert (Ge et al., 2014), in the Gobi Desert (Chen et al., 2017), in the central parts of Asia (Karakum Desert; Li and Sokolik, 2018), in the Sistan Basin (Alizadeh Choobari et al., 2013), and in the Thar Desert (Hussain et al., 2005). In North America, mineral particles emitted from the Mojave and Sonoran deserts (Hand et al., 2017) have mainly natural origins, while in the Chihuahuan Desert as well as in the Southern Great Plains the anthropogenic interference on soil can favor emissions of dust particles and, subsequently, their entrainment in the atmosphere (Hand et al., 2016). Overall, the major portion of the global dust budget arises from the deserts of the Northern Hemisphere (Ginoux et al., 2012) while mineral aerosols are also emitted in Australia (Ekström et al., 2004), Southern Africa (Bryant et al., 2007; Vickery et al., 2013), and South America (Gassó and Torres, 2019) but to a lesser extent. At a global scale, most of the entrained dust loads in the atmosphere originate from tropical and subtropical arid regions; however, about 5 % of the global dust budget consists of particles emitted from high-latitude sources (Bullard et al., 2016).

Dust plays a key role in several aspects of the Earth system, such as climate (e.g., Lambert et al., 2013; Nabat et al., 2015) and weather (Pérez et al., 2006; Gkikas et al., 2018, 2019). The key role of dust consists in the perturbation of the Earth–atmosphere system radiation budget (Sokolik and Toon, 1996; Haywood and Bucher, 2000) by mineral particles, the productivity of oceanic waters (Jickells et al., 2005) and terrestrial ecosystems (Okin et al., 2004), and the effects on human health (Kanatani et al., 2010; Kanakidou et al., 2011; Pérez García-Pando et al., 2014; Du et al., 2016; Querol et al., 2019). Dust is characterized by a pronounced temporal and spatial variability due to the heterogeneity of the emission, transport, and deposition processes governing its life cycle (Schepanski, 2018). A variety of atmospheric circulation mechanisms, spanning from local to planetary scales, are responsible for the uplifting of erodible particles from bare soils (Koch and Renno, 2005; Knippertz et al., 2007; Klose and Shao, 2012; Fiedler et al., 2013) and their subsequent transport (Husar et al., 2001; Prospero and Mayol-Bracero, 2013; Yu et al., 2015; Flaounas et al., 2015; Gkikas et al., 2015), accumulation, and removal (Zender et al., 2003; Ginoux et al., 2004) from the atmosphere.

Given the scientific importance of dust in the Earth system and the numerous socioeconomic impacts (Stefanski and Sivakumar, 2009; Weinzierl et al., 2012; Kosmopoulos et al., 2018), there is a need to monitor and forecast dust loads at different spatiotemporal scales. Contemporary satellite observations, available over long-term periods, have proven to be a powerful tool in such efforts as they provide wide spatial coverage, relatively high sampling frequency, and considerably high accuracy. Spaceborne retrievals have been widely applied in aerosol research for the description of dust load features and their evolution (e.g., Kaufman et al., 2005; Liu et al., 2008; Peyridieu et al., 2013; Rashki et al., 2015; Gkikas et al., 2013, 2016; Marinou et al., 2017; Proestakis et al., 2018). Even more accurate aerosol observations, although locally restricted, are derived by ground-based platforms consisting of Sun photometers, lidars, and in situ instruments. Based on these measurements, the columnar optical and microphysical properties of mineral particles have been analyzed extensively (Giles et al., 2012), altitude-resolved information of optical properties has provided insight about the dust vertical distribution (Mamouri and Ansmann, 2014), and a comprehensive description of dust optical, microphysical, and chemical properties has been achieved from surface and aircraft in situ instruments (Rodríguez et al., 2012; Liu et al., 2018). Finally, through the deployment of atmospheric–dust models (e.g., Pérez et al., 2011; Haustein et al., 2012), global (e.g., Ginoux et al., 2004) and regional (e.g., Basart et al., 2012) displays of the dust burden are provided.

Traditionally, observations have been utilized to evaluate and eventually constrain model performance. Observations are increasingly used in data assimilation (DA) schemes for aerosol forecast initialization (Di Tomaso et al., 2017) and the development of reanalysis data sets (Benedetti et al., 2009; Lynch et al., 2016; Gelaro et al., 2017). The most exploited reanalysis data sets in dust-related studies are the Modern Era Retrospective analysis for Research and Applications aerosol reanalysis (MERRAero; Buchard et al., 2015) and its evolution to the Modern-Era Retrospective analysis for Research and Applications version 2 (MERRA-2; Gelaro et al., 2017), and the Copernicus Atmosphere Monitoring Service Reanalysis (CAMSRA; Inness et al., 2019) and its predecessor monitoring atmospheric composition and climate (MACC; Inness et al., 2013). Current reanalysis data sets provide information about dust aerosols at high temporal resolution and decadal timescales. However, even though AOD observations are assimilated, the performance of the simulated outputs is partly model driven, and their resolution is relatively coarse.

The overarching goal of the present study is to describe the development of the ModIs Dust AeroSol (MIDAS) data set, providing dust optical depth (DOD) over a 15-year period (2003–2017). The powerful element of this product is its availability at fine spatial resolution (0.1^∘ × 0.1^∘) and on a daily basis and the provision of full global coverage (i.e., both over land and ocean). Ginoux et al. (2012) analyzed DOD at the same spatial resolution and for a long-term period but restricted the analysis to continental surfaces as the scientific focus was put on the identification of natural and anthropogenic dust sources. Voss and Evan (2020) combined satellite (MODIS; Advanced Very High Resolution Radiometer – AVHRR) aerosol retrievals and MERRA-2 winds to analyze DOD at a coarse spatial resolution (1^∘ × 1^∘) for extended time periods. Cloud–Aerosol LIdar with Orthogonal Polarization (CALIOP)-based vertical dust backscatter and extinction profiles, along with the respective column integrated DODs at 1^∘ × 1^∘ spatial resolution, are distributed via the LIdar climatology of vertical Aerosol Structure for space-based lidar simulation (LIVAS) database (Amiridis et al., 2015). Taking advantage of the spectral signature of dust at thermal infrared (TIR) wavelengths, DOD is also provided by an Infrared Atmospheric Sounding Interferometer (IASI; Vandenbussche et al., 2013; Capelle et al., 2018; Clarisse et al., 2019) and Spinning Enhanced Visible and Infrared Imager (SEVIRI; Ackerman, 1997) instruments aboard the polar-orbiting Meteorological Operational (METOP) satellites and the geostationary Meteosat Second Generation (MSG) satellite, respectively. In this case, the conversion of DOD from TIR to the mid-visible spectrum range is subjected to several assumptions related to size and other properties. Dust observations from IASI are provided at a global scale twice per day, and those of SEVIRI cover the hemisphere centered at the prime meridian over the Equator every 15 min (Schepanski et al., 2012). Thanks to their high sampling frequency, fine spatial resolution, and long-term availability, the aforementioned data sets have been used for the identification of dust source activation across North Africa (Schepanski et al., 2007; Vandenbussche et al., 2020). Based on the current status described above, MIDAS data set expands, complements, and upgrades existing observational capabilities of dust aerosols, being suitable for research studies related to climatology, model evaluation, and data assimilation.

For the development of the fine-resolution MIDAS DOD, a synergy of MODIS-Aqua (Sect. 2.1), MERRA-2 (Sect. 2.2), LIVAS (Sect. 2.3), and AEronet RObotic NETwork (AERONET; Sect. 2.4) aerosol products has been deployed by exploiting the strong capabilities of each data set. Based on the applied methodology (Sect. 3.1), the DOD is calculated by the product of MODIS-Aqua Level 2 AOD and the colocated DOD-to-AOD ratio from MERRA-2. The uncertainty of the DOD is calculated by combining the uncertainties of MODIS AOD and MERRA-2 dust fraction (MDF), using AERONET and LIVAS, respectively, as a reference (Sect. 3.2). We thoroughly compare MDF against the LIVAS dust portion in Sect. 4.1. The MIDAS DOD is evaluated against AERONET in Sect. 4.2 and compared with MERRA-2 and LIVAS DODs in Sect. 4.3. In Sect. 4.4, we provide the annual and seasonal global geographical distributions of the MIDAS DOD as a demonstration of the developed product. Finally, the main findings are summarized and the conclusions are drawn in Sect. 5.

2 Data sets

2.1 MODIS

The MODerate resolution Imaging Spectroradiometer (MODIS) is a passive sensor measuring the top of atmosphere (TOA) reflectance in order to retrieve aerosol optical depth (AOD), among other aerosol optical properties, at various wavelengths spanning from the visible to the near-infrared spectrum range. MODIS, mounted on the NASA's twin polar satellites of Terra and Aqua, has acquired high-quality aerosol data since 2000 and 2002, respectively, while, thanks to its wide swath (∼2330 km), providing near-global observations almost on a daily basis. The derivation of AOD is achieved through the implementation of two retrieval algorithms based on the Dark Target (DT) approach, valid over oceans (Remer et al., 2002, 2005, 2008) and vegetated continental areas (Levy et al., 2007a, b, 2010) but relying on different assumptions and bands, and the Deep Blue (DB) approach (Hsu et al., 2004; Sayer et al., 2013), providing retrievals over all cloud-free and snow-free land surfaces, including arid and semi-arid surfaces. MODIS data sets are organized into various collections, depending on the version of the retrieval algorithms, and into a number of levels, depending on their spatial and temporal resolution. For our purposes, we are utilizing Collection 6.1 (C061) MODIS-Aqua Level 2 (L2) retrievals over the period 2003–2017, which are reported at 5 min swath granules (Levy et al., 2013) and are accessible from the Level 1 and Atmosphere Archive and Distribution System (LAADS) Distributed Active Archive Center (DAAC) (https://ladsweb.modaps.eosdis.nasa.gov/, last access: 17 December 2020). All the updates applied in the latest version of MODIS DB and DT retrievals, with respect to Collection 6, are provided in the relevant technical documents available at the Atmosphere Discipline Team Imager Products web page (https://atmosphere-imager.gsfc.nasa.gov/documentation/collection-61, last access: 17 December 2020).

Each MODIS swath is composed of 203 × 135 retrievals, with increasing pixel size from the nadir view (10 km × 10 km) towards the edge of the satellite scan (48 km × 20 km) and to which a quality assurance (QA) flag is assigned (Hubanks et al., 2018). More specifically, these bit values represent the reliability of the algorithm output and are equal to 0 (no confidence), 1 (marginal), 2 (good), and 3 (very good). MODIS AOD retrievals are acquired based on different algorithms according to the underlying surface type. In order to fill observational gaps attributed to the assumptions or limitations of the applied MODIS algorithms, the DT ocean (QA ≥ 1), DT land (QA = 3), and DB land (QA ≥ 2) AOD retrievals are merged based on the normalized difference vegetation index (NDVI) and the highest accuracy criterion (Sayer et al., 2014). This merged AOD is stored in the scientific data set (SDS) named AOD_550_Dark_Target_Deep_Blue_Combined, which is extracted and processed for the needs of the current work. Finally, two quality filtering criteria are applied to the raw MODIS AODs for eliminating observations which may be unreliable. AODs associated with cloud fraction (CF) higher than 0.8, and those with no adjacent retrievals, are masked out, following the recommendations of previous studies (Anderson et al., 2005; Zhang and Reid, 2006; Hyer et al., 2011; Shi et al., 2011). The first criterion is associated with the potential cloud contamination on AODs while the second one discards suspicious retrievals from the data set.

2.2 MERRA-2

The Modern-Era Retrospective Analysis for Research and Applications version 2 (MERRA-2), developed by the NASA Global Modeling and Assimilation Office (GMAO), is the first atmospheric reanalysis spanning over the new modern satellite era (1980 onward) in which aerosol–radiation interactions and the two-way feedbacks with atmospheric processes are taken into account (Gelaro et al., 2017). The key components of MERRA-2 (Buchard et al., 2017) are the Goddard Earth Observing System (GEOS-5; Molod et al., 2015), which is radiatively coupled to the Goddard Chemistry Aerosol Radiation and Transport model (GOCART; Chin et al., 2002; Colarco et al., 2010) and the 3D variational (3DVar) Gridpoint Statistical Interpolation (GSI) analysis system (Wu et al., 2002).

The GOCART aerosol module simulates emission, sinks, removal mechanisms (dry deposition and gravitational settling, large-scale wet removal and convective scavenging), and the chemical processes of five aerosol species, namely dust, sea salt, sulfate, and black and organic carbon. Their optical properties are based on the updated Optical Properties of Aerosols and Clouds (OPAC) database (Hess et al., 1998), incorporating dust nonspherical shape (Meng et al., 2010; Colarco et al., 2014), and are calculated according to Colarco et al. (2010). For coarse particles (i.e., dust and sea salt), five noninteracting size bins are considered for which the emissions are driven by the wind speed, based on the parameterizations of Marticorena and Bergametti (1995), for dust and on the modified version of Gong (2003) for sea salt. Both hydrophobic and hydrophilic black carbon (BC) and organic carbon (OC) emitted from anthropogenic activities (i.e., fossil fuel combustion) and natural processes (i.e., biomass burning) are considered. Regarding sulfate aerosols (SO₄), these are either primarily emitted or are formed by the chemical oxidation of sulfur dioxide gas (SO₂) and dimethyl sulfide (DMS). Until 2010, daily emissions of eruptive and degassing volcanoes were derived from the AeroCom Phase II project (Diehl et al., 2012; http://aerocom.met.no/, last access: 17 December 2020) and afterwards only a repeating annual cycle of degassing volcanoes was included in MERRA-2. The hygroscopic growth of sea salt, sulfate, and hydrophilic carbonaceous aerosols is determined by the simulated relative humidity (RH), and the subsequent modification of particles' shape and composition is taken into account in computations of particles' fall velocity and optical parameters (Randles et al., 2017). A detailed description of the emission inventories, along with the global climatological maps representative for the period 2000–2014, are given in Randles et al. (2017).

MERRA-2 is a multidecadal reanalysis in which a variety of meteorological and aerosol observations are jointly assimilated (Gelaro et al., 2017). The former group of observations consists of ground-based and spaceborne atmospheric measurements/retrievals summarized in Table 1 of Gelaro et al. (2017), while the full description is presented in McCarty et al. (2016). For aerosol data assimilation, the core of the utilized satellite data is coming from the MODIS instrument multichannel radiances in addition to observational geometry parameters, cloud fraction, and ancillary wind data. Over oceans, AVHRR radiances, from January 1980 to August 2002, are used as well, and over bright surfaces (albedo > 0.15), the non-bias-corrected AOD retrieved for the Multiangle Imaging SpectroRadiometer (MISR; Kahn et al., 2005) is assimilated from February 2000 to June 2014. Apart from spaceborne radiances and retrievals, the Level 2 (L2) quality-assured AERONET measurements (1999–October 2014; Holben et al., 1998) are integrated in the MERRA-2 assimilation system (Goddard Aerosol Assimilation System – GAAS), which is presented in Randles et al. (2017; Sect. 3). The cloud-free MODIS radiances (DT algorithm; Collection 5) and AVHRR radiances (above oceanic regions) are used for the derivation of bias-corrected AODs, via a neural net retrieval (NNR), adjusted to the log-transformed AERONET AODs. It must be clarified that only the MERRA-2 AOD is directly constrained by the observations, while the model's performance (background forecast) and data assimilation structure (parameterization of error covariances) are responsible for the aerosol speciation, among other aerosol diagnostics (Buchard et al., 2017).

In the present study, we use the columnar MERRA-2 total and dust AOD at 550 nm in order to calculate the contribution, in optical terms, of mineral particles to the overall load. The computed dust-to-total AOD ratio (i.e., MDF) is evaluated against LIVAS and then used for the derivation of MIDAS DOD. MERRA-2 products (M2T1NXAER files; V5.12.4; aerosol diagnostics) have been downloaded from the NASA Goddard Earth Sciences (GES) Data and Information Services Center (DISC) server (https://disc.gsfc.nasa.gov/, last access: 17 December 2020) and are provided as hourly averages at 0.5^∘ × 0.625^∘ lat–long spatial resolution.

2.3 LIVAS

The European Space Agency (ESA) LIVAS database (Amiridis et al., 2015) contains a pure dust satellite-based product spanning from 2007 to 2015, which has been derived from the Cloud–Aerosol Lidar with Orthogonal Polarization (CALIOP) sensor on board the Cloud–Aerosol Lidar and Infrared Pathfinder Satellite Observation (CALIPSO) satellite. This active sensor has acquired altitude-resolved observations of aerosols and clouds since mid-June 2006 (Winker et al., 2010). CALIPSO, flying in the A-Train constellation (Stephens et al., 2002), provides almost simultaneous observations with Aqua, thus making their synergistic implementation for aerosol research feasible and powerful. CALIOP, an elastic backscatter two-wavelength polarization-sensitive Nd:YAG lidar in a near-nadir-viewing geometry (since 28 November 2007; 3^∘ off nadir), emits linearly polarized light at 532 and 1064 nm and detects the copolar components at 532 and 1064 nm and the cross-polar component at 532 nm, relative to the laser polarization plane (Hunt et al., 2009). Based on the attenuated backscatter profiles (Level 1B) and the implementation of retrieval algorithms (Winker et al., 2009), aerosol and/or cloud profiles and layer products are provided at various processing levels (Tackett et al., 2018). CALIOP Level 2 (L2) aerosol and cloud products are provided at a uniform spatial resolution along horizontal (5 km) and vertical (60 m) dimensions. Detectable atmospheric features are first categorized to aerosols or clouds and are then further discriminated into specific subtypes according to Vaughan et al. (2009). For aerosols, in Version 3 used here, six subtypes are considered, consisting of clean marine, dust, polluted continental, clean continental, polluted dust, and smoke (Omar et al., 2009). Based on the aerosol subtype classification, specific extinction-to-backscatter ratios (lidar ratio – LR) are applied for the provision of extinction coefficient profiles along the CALIPSO orbit track (Young and Vaughan, 2009).

In this study, we use the CALIOP pure dust product available in the aforementioned LIVAS database (hereafter called the LIVAS data set), which has been developed according to the methodology described in Amiridis et al. (2013) and updated in Marinou et al. (2017). The aforementioned technique relies on the incorporation of aerosol backscatter coefficient profiles and a depolarization ratio, providing a strong evidence of dust presence, due to mineral particles' irregular shape (Freudenthaler et al., 2009; Burton et al., 2015; Mamouri and Ansmann, 2017), thus allowing the separation of the dust component from aerosol mixtures. The LIVAS data set is obtained by applying appropriate regionally dependent LR values (see Fig. S1; Marinou et al., 2017; Proestakis et al., 2018 and references within), instead of the raw universal CALIOP dust LR (40 sr; Version 3), which are multiplied with the dust backscatter coefficient profiles at 532 nm in order to calculate the corresponding extinction coefficient profiles. After a series of strict quality screening filters (Marinou et al., 2017), the columnar total, dust, and nondust optical depths, as well as the DOD-to-AOD ratio over the period 2007–2015, are aggregated at 1^∘ × 1^∘ grid cells covering the whole globe. The performance of the LIVAS pure DOD product has been assessed against AERONET over North Africa and Europe (Amiridis et al., 2013) revealing a substantial improvement when the abovementioned methodological steps are applied. The LIVAS pure DOD product has been utilized in a variety of research studies, such as the assessment of dust outbreaks (Kosmopoulos et al., 2017; Solomos et al., 2018) and phytoplankton growth (Li et al., 2018), the 4D description of mineral loads over long-term periods (Marinou et al., 2017; Proestakis et al., 2018), the evaluation of dust models (Tsikerdekis et al., 2017; Georgoulias et al., 2018; Konsta et al., 2018), as well as the evaluation of new satellite products (Georgoulias et al., 2016).

2.4 AERONET

Ground-based observations acquired from the AEronet RObotic NETwork (AERONET; Holben et al., 1998) have been used as reference in this work in order to evaluate the accuracy of the MIDAS DOD product. The evaluation analysis has been performed by utilizing the almucantar (inversion) retrievals, providing information for the total aerosol amount (AOD) and for other microphysical (e.g., volume size distribution) and optical (e.g., single scattering albedo) properties (Dubovik and King, 2000; Dubovik et al., 2006). In the present study, focus is put on the aerosol optical properties retrieved at four wavelengths (440, 675, 870, and 1020 nm), utilizing spectral AODs and sky (diffuse) radiances as inputs. More specifically, we used Version 3 (V3) AERONET data (Giles et al., 2019; Sinyuk et al., 2020) of AOD (for total and coarse aerosols), the Ångström exponent (α), and single scattering albedo (SSA). For the amount (AOD) and size (α)-related optical parameters, only quality assured retrievals (i.e., Level 2; L2) are used, whereas for the SSA, the L2 and Level 1.5 (L1.5) observations are merged in order to ensure maximum availability. Unfavorable atmospheric conditions or restrictions on solar geometry result in a reduced amount of inversion outputs compared to the availability of Sun-direct measurements or the spectral deconvolution algorithm (SDA; O'Neill et al., 2003) retrievals. Even though the aforementioned AERONET data provide information about aerosol size (i.e., Ångström exponent) or coarse AOD (from SDA retrievals), the optimum approach for identifying dust particles and discriminating them from other coarse particles (i.e., sea salt) requires the use of SSA, as discussed in the next paragraph.

Through the combination of the selected optical properties from almucantar retrievals, we achieved the spectral matching between ground-based and spaceborne observations and the determination of DOD on AERONET retrievals. Regarding the first part, the α_{440–870 nm} and AOD_870 nm values are applied in the Ångström formula in order to interpolate the AERONET AOD at a common wavelength (i.e., 550 nm) with MODIS. For the evaluation of DOD, a special treatment of AERONET retrievals is required in order to identify conditions where dust particles either only exist in or clearly dominate over other aerosol species. The vast majority of previous studies (e.g., Fotiadi et al., 2006; Toledano et al., 2007; Basart et al., 2009) have relied on the combination of AOD and α for aerosol characterization, associating the presence of mineral particles with low α levels and considerable AODs. Here, we are keeping records where the α_{440–870 nm}≤0.75 and SSA_675 nm−SSA_440 nm>0, without taking into account the aerosol optical depth. The first criterion ensures the predominance of coarse aerosols, while the second one serves as an additional filter for discriminating dust from sea salt particles, taking advantage of the specific spectral signature of SSA (i.e., decreasing absorptivity for increasing wavelengths in the visible spectrum) in pure or dust-rich environments (Giles et al., 2012).

Then, from the coarse AODs at 440, 675, and 870 nm we calculate the corresponding α, which is applied in order to obtain the AERONET coarse AOD at 550 nm. This constitutes the AERONET-derived DOD, assuming that the contribution of fine dust particles (particles with radii less than the inflection point in the volume size distribution) is small. Likewise, through this consideration, any potential contamination from small-size particles of anthropogenic or natural origin (e.g., biomass burning), which is likely far away from the sources, is tempered or avoided.

3 Methods

3.1 Derivation of dust optical depth on MODIS swaths

The core concept of our approach is to derive DOD on MODIS L2 retrievals, provided at a fine spatial resolution, via the synergy with the MERRA-2 products. More specifically, the MERRA-2 dust fraction (MDF) to total AOD_550 nm (Eq. 1) is multiplied with the MODIS AOD_550 nm in order to calculate DOD_550 nm at swath level (Eq. 2).

$$\begin{array}{}\text{(1)}& {\displaystyle }\mathrm{MDF}={\displaystyle \frac{{\mathrm{AOD}}_{\mathrm{DUST};\mathrm{MERRA}-\mathrm{2}}}{{\mathrm{AOD}}_{\mathrm{TOTAL};\mathrm{MERRA}-\mathrm{2}}}}\text{(2)}& {\displaystyle }{\mathrm{DOD}}_{\mathrm{MODIS}}={\mathrm{AOD}}_{\mathrm{MODIS}}\times \mathrm{MDF}\end{array}$$

To achieve that, the data sets are colocated temporally and spatially. MERRA-2 outputs are provided at a coarse spatial resolution (0.5^∘ × 0.625^∘), in contrast to MODIS-Aqua observations (10 km × 10 km). MODIS swaths are composed by 203 × 135 retrievals, and for each one of them, we compute the nearest distance from the MERRA-2 grid points, considering the closest hourly time step to the MODIS overpass time. Then, the MDF is used to calculate the DOD from the AOD on MODIS swath native grid. Our approach purposely avoids the inclusion of additional optical properties providing information on aerosol size (α), available from MODIS, and absorptivity (Aerosol Index; Torres et al., 1998), available from Ozone Monitoring Instruments (OMIs), that are characterized by inherent limitations. Previous evaluation studies (Levy et al., 2013; Sayer et al., 2013) have shown that size parameters acquired by MODIS are highly uncertain, particularly over land and at low AOD conditions. In addition, since early 2008, the OMI sensor has lost half of its swath due to the row anomaly issue (Torres et al., 2018), thus hampering the MODIS–OMI colocation when it is attempted at a fine spatial resolution.

3.2 Uncertainty estimation

As expressed in Eq. (2), the MIDAS DOD results from the product of MODIS AOD and MDF. The uncertainty of the DOD product (Δ(DOD)) accounts for the corresponding uncertainties in the AOD and MDF, which are calculated using AERONET and LIVAS, respectively, as a reference. The mathematical expression of the Δ(DOD), given in Eq. (3), as follows, results from the implementation of the product rule on Eq. (2).

$$\begin{array}{}\text{(3)}& \mathrm{\Delta }\left(\mathrm{DOD}\right)=\mathrm{\Delta }\left(\mathrm{AOD}\right)\times \mathrm{MDF}+\mathrm{AOD}\times \mathrm{\Delta }\left(\mathrm{MDF}\right).\end{array}$$

The term Δ(AOD) defines the expected error (EE) confidence envelope within which ∼68 % of the MODIS–AERONET AOD differences are expected to fall. This term varies, depending on the applied MODIS aerosol retrieval algorithm.

For each of the two DT retrieval algorithms, we use the corresponding linear equations expressing Δ(AOD) with respect to AERONET AOD over the ocean (Levy et al., 2013; Eq. 4) and land (Levy et al., 2010; Eq. 5). For the DB AOD land retrievals (Eq. 6), we use the formula for prognostic uncertainty estimates given in Sayer et al. (2013) but with updated coefficients a and b for C061 data varying between vegetated and arid surface types (https://atmosphere-imager.gsfc.nasa.gov/sites/default/files/ModAtmo/modis_deep_blue_c61_changes2.pdf, last access: 17 December 2020). More specifically, over vegetated land, a and b are equal to 0.079 and 0.67, respectively, while the corresponding values over barren soils are equal to 0.12 and 0.61. The land cover classes have been extracted from the International Geosphere–Biosphere Programme (IGBP) database available via the MCD12C1 data (https://lpdaac.usgs.gov/products/mcd12c1v006/, last access: 17 December 2020). For the merged (DB and DT) land AOD, the uncertainty is estimated via the square root of the quadrature sum of the DT land and DB land uncertainties divided by two (Eq. 7).

$$\begin{array}{}\text{(4)}& {\displaystyle }\mathrm{\Delta }\left({\mathrm{AOD}}_{\mathrm{DT}-\mathrm{Ocean}}\right)=\pm \left(\mathrm{0.10}\times \mathrm{AOD}+\mathrm{0.04}\right)\text{(5)}& {\displaystyle }\mathrm{\Delta }\left({\mathrm{AOD}}_{\mathrm{DT}-\mathrm{Land}}\right)=\pm \left(\mathrm{0.15}\times \mathrm{AOD}+\mathrm{0.05}\right)\text{(6)}& {\displaystyle }\mathrm{\Delta }\left({\mathrm{AOD}}_{\mathrm{DB}-\mathrm{Land}}\right)=\pm \left({\displaystyle \frac{a+b\times \mathrm{AOD}}{\mathrm{AMF}}}\right)\text{(7)}& {\displaystyle }\begin{array}{rl}\mathrm{\Delta }\left({\mathrm{AOD}}_{\mathrm{DTDB}-\mathrm{Land}}\right)& =\\ & \pm {\displaystyle \frac{\sqrt{{\left[\mathrm{\Delta }\left({\mathrm{AOD}}_{\mathrm{DT}-\mathrm{Land}}\right)\right]}^{\mathrm{2}}+{\left[\mathrm{\Delta }\left({\mathrm{AOD}}_{\mathrm{DB}-\mathrm{Land}}\right)\right]}^{\mathrm{2}}}}{\mathrm{2}}}.\end{array}\end{array}$$

Before proceeding with the calculation of the Δ(DOD), a few key aspects must be highlighted for the sake of clarity. In Eqs. (4) and (5), the AOD uncertainty is defined as a diagnostic error, since it is calculated utilizing AERONET as a reference. Here, we are using the same equations, replacing AERONET AODs with those given by MODIS. This relies on the fact (results not shown here) that their averages from a global perspective are almost unbiased; however, at regional level, small negative or positive offsets (lower than 0.05 in absolute terms) are recorded in the vast majority of AERONET sites, thus supporting our argument. For the ocean AOD uncertainty, the defined EE margins (Levy et al., 2013) have been modified in order to sustain symmetry by keeping the upper bound (i.e., thus including more than 68 % of expected MODIS–AERONET pairs within the EE). Sayer et al. (2013) estimated the uncertainty of DB AOD by taking into account the geometric air mass factor (AMF) resulting from the sum of the reciprocal cosines of the solar and viewing zenith angles (Eq. 6).

The LIVAS dust fraction is our reference for estimating the MDF uncertainty. The analysis is performed at 1^∘ × 1^∘ spatial resolution, considering only grid cells in which both MERRA-2 and LIVAS DODs are higher than or equal to 0.02. According to this criterion, more than 450 000 LIVAS–MERRA-2 colocated pairs have been found that are sorted (ascending order) based on MDF (ranging from 0 to 1) and are then grouped in equal size bins containing 20 000 data for each subsample. For every group, we computed the median MDF (x axis) and the 68th percentile of the absolute MERRA-2–LIVAS dust fraction (y axis), and then we found the best polynomial fit (Eq. 8).

$$\begin{array}{}\text{(8)}& \begin{array}{rl}\mathrm{\Delta }\left(\mathrm{MDF}\right)& =\pm \left(\mathrm{2.282}\times {\mathrm{MDF}}^{\mathrm{4}}-\mathrm{6.222}\times {\mathrm{MDF}}^{\mathrm{3}}\right.\\ & \left.+\mathrm{4.700}\times {\mathrm{MDF}}^{\mathrm{2}}-\mathrm{0.969}\times \mathrm{MDF}+\mathrm{0.199}\right).\end{array}\end{array}$$

Depending on the selected MODIS algorithm, the appropriate combination between AOD (Eqs. 4, 5, 6, and 7) and MDF (Eq. 8) uncertainties is applied to calculate the Δ(DOD; Eq. 3) on each daily measurement (i.e., DOD) at each grid cell. These pixel-level DOD uncertainties are averaged over the entire study period and for each season, and the obtained findings will be discussed along with the global spatial patterns (Sect. 4.5) of dust optical depth in order to provide a measure of the reliability of the derived MIDAS DOD product.

4 Results

In the following sections, a series of analyses, including the evaluation of MDF with respect to LIVAS (Sect. 4.1), the evaluation of MIDAS DOD versus AERONET observations (Sect. 4.2), and an intercomparison among MIDAS, LIVAS, and MERRA-2 DODs (Sect. 4.3) is presented. All the aforementioned steps are performed in order to justify the validity of the applied methodology and to understand its limitations. In the last section (Sect. 4.4), the global annual and seasonal DOD patterns are presented as a demonstration of the MIDAS data set, and the obtained spatiotemporal features are briefly discussed. The detailed climatological study is provided in a companion paper.

4.1 Evaluation of MERRA-2 dust fraction versus LIVAS

The evaluation of the MDF is a critical step in our analysis since it is used as the scaling factor of the MODIS AOD for the derivation of the MIDAS DOD. For this reason, the corresponding columnar parameter provided by LIVAS (see Sect. 2.3) is used as reference. It must be highlighted that the only existing evaluation studies of MERRA-2 aerosol products have been performed either for specific aerosol species or limited time periods (e.g., Buchard et al., 2017; Veselovskii et al., 2018) or for the total load (e.g., Mukkavilli et al., 2019; Sun et al., 2019) in specific regions showing the ability of MERRA-2 to reproduce the integrated aerosol fields. Nevertheless, the speciation of the suspended particles, which is, to a large extent, determined by the model physics assumptions (Gelaro et al., 2017), has not been thoroughly evaluated. Therefore, the present analysis complements and further expands the existing works providing insight about the performance of MERRA-2 in terms of discriminating among aerosol types (particularly for dust) and subsequently estimating their contribution to the total atmospheric load.

Figure 1 depicts the geographical distributions of the dust-to-total AOD ratio, based on MERRA-2 (a) and LIVAS (b) averaged over the time frame (2007–2015) of the LIVAS data set. The corresponding maps of mean bias, fractional bias (FB), fractional gross error (FGE), and correlation coefficient (R) are given in Fig. 2. For consistency, we regridded the MERRA-2 data to 1^∘ × 1^∘ spatial resolution and selected the closest output to the CALIOP overpass time only during daytime hours when aerosol retrievals obtained by passive sensors at visible wavelengths are assimilated. At first glance, the spatial patterns are very similar, particularly in areas where the presence of dust is predominant. Across the dust belt (Prospero et al., 2002), the most evident deviations (MDF underestimation by ∼0.1 or 10 %) are recorded at the borders of Afghanistan and Pakistan (Dasht-e Margo and Kharan deserts) and in the Taklamakan Desert (Fig. 2a). However, the FB (Fig. 2b) and FGE (Fig. 2c) metrics (Yu et al., 2006), which are less affected by outliers compared to the bias, are close to zero (ideal score) in most of the aforementioned regions, thus indicating a very good performance of MERRA-2. In terms of temporal covariation (Fig. 2d), moderate R values (0.5–0.6) are obtained in land areas where the presence of dust is predominant, with the exception of the western parts of Sahara where the correlation levels are slightly higher than zero. Due to the complex and highly variable nature of the emission processes and, therefore, the poorer behavior of the model, the correlation tends to be smaller over the main dust sources throughout the year (right column in Fig. S2). In downwind regions of the Northern Hemisphere, particularly over the main transport pathways (i.e., Atlantic Ocean, Mediterranean and Arabian seas, and East Asia), correlation substantially increases (up to 0.9). In addition, FB and FGE metrics reveal a good performance by MERRA-2, which, however, downgrades for increasing distances from the sources due to the weaker dust contribution to the total aerosol load. An exception is observed for the mean bias along the tropical Atlantic Ocean, where the MDF is overestimated by up to 10 % in the eastern parts, in contrast to longitudes westward of 45^∘ W, where zero biases or slight underestimations (∼5 %; Caribbean Sea) are obtained.

Figure 1 Annual geographical distributions of dust contribution to total aerosol optical depth at 1^∘ × 1^∘ spatial resolution based on (a) MERRA-2 at 550 nm and (b) LIVAS at 532 nm, during daytime conditions, over the period 2007–2015.

Figure 2 Annual geographical distributions illustrating the assessment of MDF versus LIVAS dust fraction, during daytime conditions, at 1^∘ × 1^∘ spatial resolution, according to the primary skill metrics of (a) mean bias, (b) fractional bias, (c) fractional gross error, and (d) correlation coefficient representative for the period 2007–2015.

A discrepancy between the LIVAS and MERRA-2 dust portion is found in the Mojave, Sonoran, Chihuahuan desert areas extending between southwestern US and northern Mexico. As shown in Fig. 1, the dust contribution given by LIVAS in those areas is more widespread and stronger in contrast to MERRA-2, which simulates lower dust amounts over the sources (Mojave Desert) and the surrounding regions (maximized during December, January, and February, DJF, and March, April, and May, MAM; Fig. S2). According to the evaluation metrics in those areas (Fig. 2), the MDF underestimation ranges between 20 % and 50 %, negative FB (down to −1), and high FGE values (locally exceeding 1) are evident while the correlation levels are low, particularly over Mexico. In the Southern Hemisphere, the deficiency of MERRA-2 is pronounced along the western coasts of South America and in the Patagonian and Monte deserts, particularly during June, July, and August (JJA) and September, October, and November (SON; Fig. S2), with both being situated in Argentina. Similar results are found in Southern Africa, while in Australia a contrast between its western, eastern, and central parts, with slight MDF underestimations and overestimations of up to 20 % in absolute terms, respectively, is recorded (Fig. 2a). Nevertheless, the agreement between MERRA-2 and LIVAS in temporal terms is supported by the moderate-to-high R values over the hot spot regions (Fig. 2d). Outside of the main dust-affected regions, an obvious discrepancy is found in the eastern Canada and northeastern Russia, where MDF yields very low values (<20 %) in contrast to LIVAS reaching values of dust fraction up to 50 %. Due to their geographical position, the occurrence of dust loads might not be frequent there; however, their contribution to the total load can be significant under low AOD conditions, which are mainly recorded in the region. This might indicate a poor representation by MERRA-2. However, potential cloud contaminations in the lidar signals, affecting LIVAS reliability, must also be taken into account.

The obtained discrepancies are mainly driven by the partial representation of dust sources in MERRA-2, resulting in potentially underestimated dust emissions and, subsequently, lower dust contribution to the total burden. Dust is originated either from natural (arid lands, salt lakes, and glacial lakes) or from anthropogenic sources (Ginoux et al., 2012). Nevertheless, dust sources in MERRA-2 are based on Ginoux et al. (2001), accounting mostly for natural dust emission areas. This could partly explain the higher LIVAS dust contribution levels that are also evident in the seasonal distributions where the interannual variability in dust fraction is illustrated (Fig. S2). Interestingly, most of the MDF underestimations (i.e., bluish colors in Fig. 2a) are recorded in mountainous areas. Depending on the homogeneity of the atmospheric scene over regions characterized by complex topography, variations in the optical paths of subsequent CALIPSO L2 profiles considered in the LIVAS product may result in unrealistic DOD and AOD values. Previous evaluation studies (e.g., Omar et al., 2013) have shown that CALIOP underestimates AOD with respect to ground-based AERONET retrievals, particularly over desert areas (Amiridis et al., 2015), which is attributed, primarily, to the incorrect assumption of the lidar ratio (LR; Wandinger et al., 2010) and, secondarily, to the inability of the lidar to detect thin aerosol layers (particularly during daytime conditions due to the low signal-to-noise ratio). The former factor is related to aerosol type, and for Saharan dust particles, the necessary increase in LR (from 40 to 58 sr) substantially improved the level of agreement with AERONET and MODIS (Amiridis et al., 2013). Similar adjustments (increments) to the raw LR values, which are highly variable (Müller et al., 2007; Baars et al., 2016), considered in the CALIOP retrieval algorithm have been applied in the LIVAS product over other source areas of mineral particles (see Sect. 2.3; Fig. S1). An additional factor that must be taken into account is the number of MERRA-2–LIVAS pairs that are used for the metrics calculation. The corresponding global geographical distribution (Fig. S3i), calculated over the period 2007–2015, shows that in areas where the model–satellite agreement is good (Fig. 2), the number of common samples is high (>100) in contrast to regions with a low number of common samples (<50) and where the computed metrics are degraded.

In order to complete the evaluation of the MDF versus LIVAS, the dependency of the level of agreement on the spatial representativeness within the 1^∘ × 1^∘ LIVAS grid cell has also been investigated. Figure S3ii displays the long-term averaged geographical distribution of the number of CALIOP L2 profiles (up to 24) aggregated for the derivation of the LIVAS 1^∘ × 1^∘ grid cell. According to the global map, the maximum number is recorded in the latitudinal band extending from 45^∘ S to 45^∘ N, while the impact of extended clouds around the Equator is apparent. Outside this zone, the number of profiles used is mainly less than 14 and decreases towards the poles due to the enhanced cloudiness. The same evaluation metrics presented in Fig. 2 have also been computed at a planetary scale for individual classes of CALIOP L2 profiles (Fig. S4) aggregated for the derivation of the LIVAS 1^∘ × 1^∘ grid cells. Overall, about 3.4 million pairs (ALL tick named ALL in Fig. S4a) have been found over the period 2007–2015 and are almost equally distributed for bins spanning from 8 to 20, while the number of colocated data is higher in the lowermost (≤7) and uppermost (≥21) tails of the distribution. The FB (Fig. S4c), FGE (Fig. S4d), and correlation (Fig. S4e) results reveal that the consistency between MDF and LIVAS gradually improves for higher grid cell representativeness. At a global scale, MERRA-2 overestimates the dust fraction by up to 1.5 % with respect to LIVAS (Fig. S4b; ALL sample). Among the bin classes, the MDF–LIVAS differences are mostly positive and lower than ∼3 % and decrease further when at least 12 CALIOP profiles are aggregated for the derivation of the LIVAS grid cell.

4.2 Evaluation of MIDAS DOD versus AERONET

In the present section, we provide an evaluation of the MIDAS DOD against the corresponding AERONET product (Sect. 2.4). An illustration of the MODIS–AERONET colocation method (an example from aerosol optical depth without applying the criteria for DOD) is shown in Fig. S5. The obtained results at global and station level are presented in Figs. 3 and 4, respectively. As expected, the number of coincident spaceborne and ground-based DODs collected at 436 AERONET stations (red circles in Fig. 3a) is low (10 478 pairs) due to the limited amount of almucantar retrievals and the implementation of filters for the determination of DOD in AERONET data. According to the global scatterplot metrics (Fig. 3b), a very good performance of the MIDAS DOD is revealed since both data sets are well correlated (R=0.89), with MIDAS only slightly overestimating DOD compared to AERONET (0.004 or 2.7 %). Only AERONET AODs associated with α lower to or equal than 0.75 are kept for the evaluation procedure. While this threshold is higher compared to previously applied cut-off levels (e.g., Dey et al., 2004; Tafuro et al., 2006; Reid et al., 2008; Kim et al., 2011; Gkikas et al., 2016), our global scatterplot metrics are very similar when reducing α from 0.75 to 0.25 (results not shown here).

Figure 3 (a) AERONET sites where at least one pair of ground-based and spaceborne retrievals has been recorded according to the defined colocation criteria during the period 2003–2017. (b) Density scatterplot between MIDAS (y axis) and AERONET (x axis) dust optical depth at 550 nm. The solid and dashed lines stand for the linear regression fit and equal line (y=x), respectively.

Figure 4 Scatterplot metrics between MIDAS and AERONET DOD_550 nm, at station level, during the period 2003–2017. (a) Number of concurrent MIDAS–AERONET observations, (b) correlation coefficient, (c) root mean square error, and (d) bias defined as spaceborne minus ground-based retrievals. The obtained scores are presented only for sites with at least 30 MIDAS–AERONET matchups.

The evaluation analysis was also performed for each station individually. Figure 4 only depicts sites with at least 30 coincident MIDAS–AERONET observations, thus making meaningful the comparison at station level. This criterion is satisfied in 86 stations, which comprise 77 % (or 8095) of the total population of coincident DODs overall and are mostly located over dust sources and in areas affected by dust transport. Figure 4a shows the station-by-station variability in the number of common MIDAS–AERONET observations ranging from 100 to 457 (Banizoumbou, Niger) across North Africa and the Middle East, whereas in the remaining sites it is mainly lower than 70. Between the two data sets, very high R values (up to 0.98) are found in North Africa, the Middle East, outflow regions (Cabo Verde, Canary Islands, and the Mediterranean) and at distant areas (Caribbean Sea) affected by long-range transport. Across the Sahel, maximum root mean square error (RMSE) levels (up to 0.26) are recorded (Fig. 4c) due to the intense loads and strong variability in the Saharan dust plumes. Regarding biases, positive deviations of up to 0.08 are computed in most AERONET sites in the area, while the largest negative offsets (down to −0.14) are recorded at the stations of Ilorin and Djougou (near to the coasts of the Gulf of Guinea), in agreement with Wei et al. (2019b). Several reasons may explain the obtained MIDAS–AERONET differences over the abovementioned stations, taking into account that the MDF is generally well reproduced. The first one is related to the MODIS retrieval algorithm itself and, more specifically, to the applied aerosol models, surface reflectance, and cloud-screening procedures (Sayer et al., 2013). The second factor is the omission of fine DOD in AERONET data, which would likely reduce the positive biases. However, its contribution to the total dust AOD is difficult and probably impossible to be accurately quantified. Similar tendencies are found for RMSE and bias in the Middle East, where the satellite and ground-based DODs are well correlated in general. In the Mediterranean, the temporal covariation between the two data sets is quite consistent (R>0.8), with the MIDAS DOD being slightly underestimated probably due to the MDF underestimation (mainly recorded in JJA; Fig. S2).

In Asia, few stations are available with sufficient number of MIDAS–AERONET matchups in which the slight positive and the negative DOD biases (Fig. 4d) are generally consistent with those of AOD (results not shown here). This indicates that the MODIS AOD offsets are transferred to MIDAS DOD, which is also affected by the MDF underestimation (Fig. 2a). In the southwestern United States, the evaluation scores at 12 AERONET sites show a moderate performance of the derived DOD (R – 0.28–0.94; bias – −0.034–0.003; RMSE – 0.02–0.04) attributed to the inability of MERRA-2 to adequately reproduce the contribution of dust particles to the total aerosol load in optical terms. Finally, our assessment analysis in the Southern Hemisphere for stations located in Argentina, Namibia, and Australia reveals MIDAS–AERONET deviations, spanning from −0.03 (Cordoba-CETT) to 0.02 (Gobabeb), and correlations, ranging from 0.14 (Fowlers_Gap) to 0.96 (Canberra).

4.3 Intercomparison of MIDAS, MERRA-2, and LIVAS DOD products

Following the evaluation of MIDAS DOD against AERONET, the MIDAS, MERRA-2, and LIVAS DOD products are investigated in parallel. For this purpose, the MERRA-2 and MIDAS data have been regridded to 1^∘ × 1^∘ grid cells between 2007 and 2015 to match the spatial resolution and availability of LIVAS. Then, the three data sets have been colocated spatially and temporarily. The intercomparison has been performed only during daytime conditions, and the obtained findings are presented through geographical distributions (Sect. 4.3.1) and averaged monthly time series of regional, hemispherical, and planetary averages (Sect. 4.3.2). Finally, it must be clarified that our focus in this part of the analysis is on the intercomparison among the DOD products and not on interpreting their spatiotemporal features. The latter will be discussed thoroughly in a companion paper analyzing the MIDAS fine-resolution DOD data set.

4.3.1 Geographical distributions

The annual geographical distributions of LIVAS, MERRA-2, and MIDAS DODs are depicted in Fig. 5a, b, and c, respectively, while the corresponding global seasonal maps are provided in Fig. S6. Among the three data sets, both for annual and seasonal geographical distributions, it is apparent that there is a very good agreement in spatial terms, in contrast to the magnitude of the simulated (MERRA-2) and retrieved (MIDAS and LIVAS) DODs. The most evident differences in MERRA-2 (Fig. 5b) and MIDAS (Fig. 5c), with respect to LIVAS (Fig. 5a), are encountered across North Africa, forming clear patterns with positive and negative deviations over the Sahara and the Sahel, respectively. In particular, MERRA-2 DOD positive offsets mostly range from 0.04 to 0.20, while those of MIDAS–LIVAS are lower, placing our DOD product between active remote sensing retrievals and reanalysis data set. Previous studies relying on satellite (Yu et al., 2010; Kittaka et al., 2011; Ma et al., 2013) and ground-based (Schuster et al., 2012; Omar et al., 2013) observations have reported that CALIOP underestimates AOD over the Sahara. Konsta et al. (2018), who utilized higher and more realistic dust lidar ratio (55 sr; adopted also for the region in the current study) compared to the aforementioned works (40 sr), reported similar tendencies in lower magnitude against MODIS. Therefore, additional factors might contribute to the lower lidar-derived DODs over the arid regions in North Africa. For example, it has been observed that CALIOP can misclassify very intense dust layers as clouds, which, on the other hand, can attenuate the emitted laser beam significantly or totally (Yu et al., 2010; Konsta et al., 2018). All these aspects, most likely met over dust sources, act towards reducing DOD (resulting from the vertical integration of the extinction coefficient profiles) and might explain the missing hot spot by LIVAS in and around the Bodélé Depression in contrast to single-view, multi-angle, and geostationary passive satellite sensors (e.g., Banks and Bridley, 2013; Wei et al., 2019a). Across the Sahel, LIVAS provides higher DODs (mainly up to 0.2) against both simulated and satellite products. These differences might be attributed to the misrepresentation of dust sources in MERRA-2 along this zone where vegetation cover has a prominent seasonal cycle (Kergoat et al., 2017). An inaccurate representation of vegetation also impacts the surface reflectance, which, in turn, can introduce critical errors in the MODIS retrieval algorithm. Wei et al. (2019b) showed that MODIS underestimates AOD with respect to AERONET in that region, which is in agreement with the fact that the maximum MIDAS–AERONET negative DOD differences are found at Ilorin and Djougou sites (Fig. 4d).

Figure 5 Long-term (2007–2015) average geographical colocated distributions at 1^∘ × 1^∘ spatial resolution during daytime for (a) LIVAS DOD_532 nm, (b) MERRA-2 DOD_550 nm, and (c) MIDAS DOD_550 nm.

Over the eastern tropical Atlantic Ocean, the difference between LIVAS and MIDAS is negligible, whereas MERRA-2 gives lower DODs by up to 0.08. In the Middle East, MERRA-2 and MIDAS DODs are higher than in LIVAS over the Tigris–Euphrates basin, while an opposite tendency for MERRA-2 is found in the interior parts of Saudi Arabia. Lower DODs are also given by LIVAS over the arid and semi-arid regions eastwards of the Caspian Sea, including also the Aral Sea. This area of the planet is one of the most challenging for spaceborne passive observations due to the terrain complexity prohibiting the accurate characterization of the surface reflectance and type, resulting in unrealistically high MODIS AODs (Klingmüller et al., 2016; see the interactive comment posted by Andrew Sayer) that may also affect MERRA-2 via assimilation. The largest negative MIDAS–LIVAS differences (exceeding 0.2) worldwide are recorded in the Taklamakan Desert, whereas the corresponding results between MERRA-2 and LIVAS are somewhat lower. This might be attributed to an inappropriate selection (overestimation) of the lidar ratio, taking into account that over the region CALIOP mainly underestimates AOD, dust contribution to the total AOD exceeds 70 % (Proestakis et al., 2018) throughout the year and MDF shows robust consistency (Fig. 2). Eastwards of the Asian continent, the situation is reversed, and the LIVAS DODs are lower by up to 0.2 when compared to MERRA-2 and MIDAS, indicating a weaker trans-Pacific transport, predominant during the boreal spring (second row in Fig. S6), being in agreement with the findings of Yu et al. (2010) and Ma et al. (2013). In the Southern Hemisphere, negative MERRA-2–LIVAS and MIDAS–LIVAS differences are computed in Patagonia, attributed to the underperformance of MDF, which are not, however, spatially coherent. On the contrary, in the desert areas of the inland parts of Australia, there is a clear signal of positive MERRA-2–LIVAS deviations, not seen between MIDAS and LIVAS, most likely attributed to the overestimation of aerosol (dust) optical depth by MERRA-2, as recently presented by Mukkavilli et al. (2019). For a global and long-term perspective, based on ∼440 000 colocated data, MERRA-2 agrees slightly better with LIVAS than MIDAS as revealed by the correlation (R=0.74 vs. 0.71) and bias (8.2 % vs. 13.3 %) metrics.

4.3.2 Planetary, hemispherical, and regional intra-annual variability

We compared the monthly variability of the planetary (Fig. 6a) and hemispherical (Fig. 6b, c) averages of LIVAS (black curve), MERRA-2 (red curve), and MIDAS (blue curve) DODs. We note that, for each considered timescale, the averaging has been made following the upper branch shown in Fig. 5 of Levy et al. (2009), where each grid cell is first temporally averaged, and the resulting field is spatially averaged. In the Northern Hemisphere, the annual cycle of DOD is reproduced by the three data sets with maximum levels in June (0.118 for LIVAS–MERRA-2 and ∼0.126 for MIDAS) and minimum ones in November–December (0.034–0.040). Nevertheless, the most evident deviations in terms of magnitude are recorded during the high-dust seasons, with MIDAS giving slightly higher DODs than MERRA-2 and even higher than LIVAS, particularly during the boreal spring. On an annual basis (Table 1), the averaged MERRA-2 and MIDAS DODs for the Northern Hemisphere are equal to 0.056 and 0.060, respectively, and higher than the LIVAS climatological value (0.051). In the Southern Hemisphere (Fig. 6c), DODs range at very low levels (up to ∼0.011), attributed to the low amounts of mineral particles emitted from spatially restricted desert areas and the limited dust transport over oceanic regions. Despite the low annual levels (0.008; Table 1) there is an intra-annual cycle pattern not entirely commonly reproduced by the three data sets. In particular, MIDAS and MERRA-2 DODs are maximized in February, while the highest levels for LIVAS are recorded in September. For all DOD products, the minimum values are found in May, which are slightly lower than those observed during April–July (austral winter). At a global scale (Fig. 6a), the seasonal patterns of DODs are mainly driven by those of the Northern Hemisphere, particularly for LIVAS, and to a lesser degree for MIDAS and MERRA-2. More specifically, there are two peaks (∼0.055; March and June) for MIDAS and flat maximum levels (∼0.05) between March and June for MERRA-2, while there is a primary (∼0.05) and a secondary (∼0.04) maximum in June and March, respectively, in LIVAS. Even though there are month-by-month differences, the LIVAS (0.029), MERRA-2 (0.031), and MIDAS (0.033) global annual DODs are relatively close, indicating a sufficient level of agreement among the three data sets (Table 1). The obtained value for MIDAS is 10 % higher and within the uncertainty estimate of the global DOD average (0.030 ± 0.005) reported by Ridley et al. (2016).

Figure 6 Monthly variability in LIVAS (black curve), MERRA-2 (red curve), and MIDAS (blue curve) DODs, regionally averaged over (a) the whole globe (GLB), (b) the Northern Hemisphere (NHE), and (c) the Southern Hemisphere (SHE). The error bars correspond to the monthly interannual standard deviation computed during the period 2007–2015.

Table 1 Planetary (GLB), hemispherical (NHE and SHE), and regional DOD averages representative for the period 2007–2015, based on colocated LIVAS, MERRA-2, and MIDAS 1^∘ × 1^∘ data. Within the parentheses, the minimum and maximum annual values are given. The regional averages have been calculated following the upper branch (first temporal averaging and then spatial averaging) in Fig. 5 of Levy et al. (2009). The full names of the acronyms for each subregion are given in the caption of Fig. 7.

The consistency among the three DOD data sets, in terms of magnitude and temporal covariation, is highly dependent on the region of interest. Figure 7 shows the defined subdomains considered in this study, Fig. S7 depicts the corresponding intra-annual DOD time series, while Table 1 lists the computed annual averages and their minimum or maximum values between 2007 and 2015. The best agreement among MIDAS, LIVAS, and MERRA-2 is found along the tropical Atlantic Ocean. In the nearby outflow regions (i.e., ETA), considerably high DODs (>0.1) are found between January–August, with a maximum in June, as indicated by the three data sets, with slight underestimations in MERRA-2 (Fig. S7k). Over the western tropical Atlantic Ocean, the sharp increase in DOD from May to June indicates the arrival of considerable amounts of Saharan particles, which are sustained at high levels in summer and diminish during autumn and winter (Fig. S7q). This seasonal fluctuation is almost identically reproduced by the three products. Nevertheless, when the dust activity is well established in the area (i.e., boreal summer), LIVAS shows higher values than MERRA-2 and MIDAS.

Figure 7 Regional domains of the east tropical Atlantic (ETA), west tropical Atlantic (WTA), Mediterranean (MED), Gulf of Guinea (GOG), western Sahara (WSA), sub-Sahel (SSA), Bodélé Depression (BOD), northern Middle East (NME), southern Middle East (SME), central Asia (CAS), Thar Desert (THA), Taklamakan Desert (TAK), Gobi Desert (GOB), East Asia (EAS), western North Pacific (WNP), eastern North Pacific (ENP), and southwestern United States (SUS).

Across North Africa, and particularly in the Bodélé (Fig. S7a) and western Sahara (Fig. S7h), the LIVAS DODs are substantially lower when compared to MIDAS and MERRA-2. In the Bodélé, this is evident for the entire year, and in western Sahara it can be clearly seen between March and June. Similar findings are drawn, either for other source areas, such as central Asia (Fig. S7c), or outflow regions, such as the Mediterranean (Fig. S7m). In southwestern United States (SUS; Fig. S7e), the seasonal variation in DODs is in a very good agreement between MERRA-2 and MIDAS, but a positive offset is seen for the reanalysis data. On the contrary, LIVAS does not reproduce the secondary maximum in July, while it gives very high DODs in November–December, which are not reliable. Over the Taklamakan (Fig. S7f), the LIVAS DODs are higher than the corresponding MIDAS and MERRA-2 regional averages in the high-dust months (i.e., April–May) and in July. On the contrary, in the Gobi Desert, residing eastwards, the LIVAS–MIDAS agreement is very good, while MERRA-2 DODs are less variable, within the course of the year with respect to the observed values (Fig. S7b). Among the three DOD products, a very good temporal agreement is found in the Thar Desert (Fig. S7g), but there are deviations regarding the peak of July, which is higher in LIVAS (0.88) than in MIDAS (0.75) and MERRA-2 (0.48), respectively. Over downwind continental areas of East Asia (Fig. S7i), only a few exceptions break down the consistency between MIDAS and LIVAS, whereas MERRA-2 is able to reproduce the annual cycle, but underestimates the intensity of dust loads. In the southern Middle East (Fig. S7n), the reanalysis and the spaceborne lidar DODs are very well correlated within the course of the year and lower than MIDAS during February–May. In the northern parts (Fig. S7d), MIDAS gives substantially higher DODs against the well-correlated and matched values of LIVAS and MERRA-2. Over the Northern Pacific, Asian dust is transported eastwards during spring, affecting nearby (Fig. S7p) and distant (Fig. S7j) oceanic areas. The signal of this mechanism is clearly evident in MIDAS and MERRA-2 time series, in contrast to LIVAS, which exhibits substantially lower DOD maxima. Moreover, these maxima appear in the western North Pacific Ocean earlier (March) with respect to the other two data sets (April). Based on MERRA-2 and MIDAS in the sub-Sahel (Fig. S7o), a primary maximum and a secondary maximum are recorded in March and October, in agreement with ground-based visibility records (N'Tchayi Mbourou et al., 1997). LIVAS reproduces both peaks but with a weaker intensity in March compared to MIDAS and MERRA-2. However, throughout the year, the maximum LIVAS DOD is observed in June (a local maximum is also recorded in MIDAS), which might be attributed to the strong convection activity favoring the occurrence of haboobs. Saharan dust aerosols, under the impact of the northeasterly harmattan winds, are carried over the Gulf of Guinea (Fig. S7l) during the boreal winter, although DODs among the three data sets reveal a noticeable variability in terms of intensity.

4.4 MIDAS DOD global climatology

The annual and seasonal DOD patterns, representative for the period 2003–2017, are illustrated in Figs. 8 and 9, respectively. Among the desert areas of the planet, the most intense dust loads (DODs up to ∼1.2; Fig. 8a) are observed in the Bodélé Depression located in the northern Lake Chad basin (Washington et al., 2003). Over the region, these high DODs are sustained throughout the year (Fig. 9), while, due to the prevailing meteorological conditions, during MAM (Fig. 9b) and JJA (Fig. 9c), mineral particles are transported westwards, along the Sahel, contributing to the locally emitted anthropogenic dust (Ginoux et al., 2012). Substantial high climatological DODs (up to 0.6; Fig. 8a) are recorded in the western sector of Sahara, in contrast to the eastern parts, and are attributed to the accumulation of dust aerosols, primarily in JJA (Fig. 9c) and secondarily in MAM (Fig. 9b), under the impact of the Saharan heat low (Schepanski et al., 2017). Saharan dust is subjected to short-range transport, frequently affecting the nearby maritime areas of the Gulf of Guinea (Ben-Ami et al., 2009), the Mediterranean Sea (Gkikas et al., 2015), and the Red Sea (Banks et al., 2017). Nevertheless, the strongest signal of Saharan dust transport appears over the tropical Atlantic Ocean, with massive loads of mineral particles confined within the Saharan air layer (SAL; Kanitz et al., 2014), reaching the Caribbean Sea (Prospero, 1999), under the impact of the trade winds. The characteristics of the trans-Atlantic dust transport reveal a remarkable intra-annual variation (Fig. 9) as it concerns the plumes' latitudinal position, longitudinal extension, and intensity, which is at maximum levels during the boreal summer (Fig. 9c).

Figure 8 Annual geographical distributions, at 0.1^∘ × 0.1^∘ spatial resolution, of (a) the climatological DODs, (b) the average of the daily DOD uncertainties, and (c) the percentage availability of MIDAS data with respect to the entire study period spanning from 1 January 2003 to 31 December 2017. Gray coloring represents areas with an absence of data.

Figure 9 As in Fig. 8a but for (a) December–January–February (DJF), (b) March–April–May (MAM), (c) June–July–August (JJA), and (d) September–October–November (SON).

Dust activity over the Middle East is more pronounced in a zone extending from the alluvial plain of the Tigris–Euphrates river system to the southern parts of the Arabian Peninsula (Fig. 8), through the eastern flatlands of Saudi Arabia (Hamidi et al., 2013). Mineral particles emitted from these sources also affect the Persian Gulf (Giannakopoulou and Toumi, 2011) and the Red Sea (Banks et al., 2017); however, the major transport pattern is recorded across the northern Arabian Sea in JJA (Fig. 9c), when dust plumes can reach the western coasts of India (Ramaswamy et al., 2018). In the Asian continent, the Taklamakan Desert (Ge et al., 2014), situated in the Tarim Basin (northwestern China), is one of the strongest dust sources of the planet, yielding DODs up to 1 during spring (Fig. 9b). These intensities are substantially higher than those recorded in the Gobi Desert, located eastwards in the same latitudinal band, due to the different composition of the erodible soils (Sun et al., 2013). Midlatitude cyclones propagating eastwards during springtime (Fig. 9b) mobilize dust emissions from both sources, inducing uplifting and, subsequently, advection of mineral particles towards the continental East Asia (Yu et al., 2019), over the northern Pacific Ocean (Yu et al., 2008), and exceptionally over the United States (Husar et al., 2001). Other hot spots of dust activity in Asia are recorded in the central parts (Li and Sokolik, 2018) and in the Sistan Basin (Alizadeh Choobari et al., 2013). Dust aerosols originating from agricultural activities along the Indus River basin (Ginoux et al., 2012) and natural processes in the Thar Desert (Proestakis et al., 2018) result in the accumulation of mineral particles at the Pakistan–India borders, while, under favorable meteorological conditions, these loads are carried towards the Indo-Gangetic Plain mainly during the premonsoon season (Dey et al., 2004). In North America, dust production becomes more evident in the southwestern United States and northwest Mexico in regional terms and during spring within the course of the year (Fig. 9b). However, DODs are mostly lower than 0.2, with a few local exceedances, indicating relatively weak dust emissions from the natural (Mojave and Sonoran deserts; Hand et al., 2017) and anthropogenic (Chihuahuan Desert and Southern Great Plains; Hand et al., 2016) dust sources of the region. Between the two hemispheres, there is a clear contrast in DODs, which are substantially lower in the Southern Hemisphere, attributed to the weaker processes triggering dust emission from the spatially restricted deserts located in Southern Africa (Bryant et al., 2007), South America (Gassó and Torres, 2019), and in the interior parts of Australia (Prospero et al., 2002).

In addition to the global climatological DOD pattern in Fig. 8a, the average of the daily DOD uncertainties provided within the data set (not to be confused with the uncertainty of the average DOD) and the temporal availability of the MIDAS data set are shown in Fig. 8b and c, respectively. More than 70 % of daily satellite retrievals with respect to the full period are included in the calculation of the mean DODs (Fig. 8a) over the cloud-free desert areas. Over dust-affected downwind regions, the corresponding percentages range from 30 % to 60 % (Fig. 8c). As expected from Eqs. (4) to (7), daily DOD uncertainties (Fig. 8b) scale with DOD and reach up to a 0.4 annual average and 0.5 when averaged over MAM and JJA (Fig. S8) in the regions with strongest DODs.

5 Summary and conclusions

In the current study, we presented the MIDAS (ModIs Dust AeroSol) dust optical depth (DOD) data set, developed via the synergistic implementation of MODIS-Aqua AOD and dust fraction extracted from colocated MERRA-2 reanalysis outputs. The derived fine-resolution (0.1^∘ × 0.1^∘) global data set between 2003 and 2017 provides DOD both over continental and oceanic areas, in contrast to similar available satellite products restricted over land surfaces (Ginoux et al., 2012), thus making a thorough and consistent description of dust loads, not only over the sources but also over downwind regions, feasible. Reanalysis data sets, spanning through decades and available at a high temporal frequency, can fulfill such tasks; however, their coarser spatial resolution imposes a restriction when investigating the mineral loads' features at finer spatial scales. Our developed DOD product aims at complementing existing observational gaps and can be exploited in a variety of studies (e.g., climatology, trends, evaluation of atmospheric–dust models, radiative effects, and data assimilation).

The core concept of the applied methodology relies on the utilization of MODIS AOD and MERRA-2 dust fraction (MDF) for the derivation of DOD on MODIS swaths. The validity of MDF has been justified through its evaluation against reference values obtained by the LIVAS database. Over dust-abundant areas extending across the dust belt, MERRA-2 adequately reproduces the magnitude of dust portion, as indicated by the calculated primary statistics (bias, FB, and FGE), with the maximum underestimations (up to 10 %) being observed in Asian deserts. The agreement between MDF and LIVAS is reduced in the main dust regions of North America and in the Southern Hemisphere. Regarding the temporal covariation of the observed and simulated dust portions, over the period 2007–2015 moderate R values (up to 0.5) are computed above the sources and are attributed to the high spatiotemporal variability of the emission processes. On the contrary, the correlation increases substantially (up to 0.9) over maritime downwind regions (tropical Atlantic Ocean, North Pacific Ocean, Arabian Sea, and Mediterranean Sea) where the main dust transport pathways are recorded. Apart from the geographical dependency of the level of agreement between MDF and LIVAS dust fraction, we also investigated the impact of the spatial representativeness of the CALIOP observations. Through this analysis, we revealed that, for an increasing number of CALIOP L2 profiles (ranging from 1 to 23), that are aggregated for the derivation of the 1^∘ × 1^∘ LIVAS grid cell, the computed metrics converge towards the ideal scores.

Finally, the obtained MIDAS DOD was evaluated against AERONET retrievals and compared with LIVAS and MERRA-2 DODs. AERONET observations were processed to minimize the contribution of other aerosol species, making also the assumption that dust loads mainly consist of coarse particles (their radii are larger than the defined inflection point). Overall, the agreement between ∼10 500 MIDAS–AERONET pairs is very high (R=0.89), whereas the satellite DODs are higher by 2.7 % with respect to the ground-based ones. At station level, the R values are mainly above 0.8 at most sites in the Northern Hemisphere (except western US), while they are mostly lower than 0.5 in the Southern Hemisphere. Moreover, positive MIDAS–AERONET deviations (up to 0.2) are mainly encountered in North Africa and the Middle East, in contrast to negative values (down to −0.14) recorded at the remaining sites. Based on the annual and seasonal global DOD patterns corresponding to the period 2007–2015, the locations with the maximum DODs are in a good agreement among the three data sets. Nevertheless, in many regions (e.g., Bodélé, sub-Sahel, and northern Pacific Ocean) there are deviations in the intensity of dust loads, attributed to the inherent weaknesses of DOD derivation techniques based on different approaches. Despite the regional dependency of deviations among the three data sets, the colocated global long-term averaged DOD is very similar (0.029 for LIVAS, 0.031 for MERRA-2, and 0.033 for MIDAS) and close to that reported (0.030) in Ridley et al. (2016). In the Southern Hemisphere, the corresponding levels are equal to 0.008 for the three data sets, whereas in the Northern Hemisphere, LIVAS DODs (0.051) are lower with respect to MIDAS (0.060) and MERRA-2 (0.056).

As a demonstration of the MIDAS data set, a brief discussion about dust load regime at global scale is made by analyzing the annual and seasonal DOD patterns. The most pronounced dust activity recorded in the Bodélé Depression (DODs up to ∼1.2), across the Sahel (DODs up to 0.8), in western parts of the Sahara (DODs up to 0.6), in the eastern parts of the Arabian Peninsula (DODs up to ∼1), along the Indus River basin (DODs up to 0.8), and in the Taklamakan Desert (DODs up to ∼1). On the contrary, the weaker emission mechanisms triggering dust mobilization over the spatially limited sources of Patagonia, Southern Africa, and the interior arid areas of Australia do not favor the accumulation of mineral particles at large amounts (DODs up to 0.4 at local hot spots), even during high-dust seasons. Over oceans, the main pathways of long-range dust transport are observed along the tropical Atlantic and the northern Pacific, revealing a remarkable variation, within the course of the year, in terms of intensity, latitudinal position, and range. Finally, the Mediterranean and the Arabian seas are affected by advected dust plumes originating from North Africa and the Middle East, respectively. Based on the performed uncertainty analysis, the MIDAS DOD product is highly reliable over dust-rich regions and becomes more uncertain in areas where the existence of dust loads is not frequent.

The exploitation of the MIDAS DOD product will be expanded on in other studies. At present, focus is given to (a) the DOD climatology over dust sources and downwind regions, (b) the implementation of the MIDAS data set in the DA scheme of the MONARCH model (Di Tomaso et al., 2017), (c) the estimation of dust radiative effects and the associated impacts on solar energy production, in North Africa and Middle East, upgrading the work of Kosmopoulos et al. (2018), and (d) the analysis of global and regional trends of dust loads.