15-year variability of desert dust optical depth on global and regional scales

This study aims to investigate global, regional and seasonal temporal dust changes as well as the effect of dust particles on total aerosol loading using the ModIs Dust AeroSol (MIDAS) fine-resolution dataset. MIDAS delivers dust optical depth (DOD) at fine spatial resolution ($\mathrm{0.1}{}^{\circ }\times \mathrm{0.1}{}^{\circ }$) spanning from 2003 to 2017. Within this study period, the dust burden increased across the central Sahara (up to 0.023 yr⁻¹) and Arabian Peninsula (up to 0.024 yr⁻¹). Both regions observed their highest seasonal trends in summer (up to 0.031 yr⁻¹). On the other hand, declining DOD trends are encountered in the western (down to −0.015 yr⁻¹) and eastern (down to −0.023 yr⁻¹) Sahara, the Bodélé Depression (down to −0.021 yr⁻¹), the Thar (down to −0.017 yr⁻¹) and Gobi (down to −0.011 yr⁻¹) deserts, and the Mediterranean Basin (down to −0.009 yr⁻¹). In spring, the most negative seasonal trends are recorded in the Bodélé Depression (down to −0.038 yr⁻¹) and Gobi Desert (down to −0.023 yr⁻¹), whereas they are in the western (down to −0.028 yr⁻¹) and the eastern Sahara (down to −0.020 yr⁻¹) and the Thar Desert (down to −0.047 yr⁻¹) in summer. Over the western and eastern sector of the Mediterranean Basin, the most negative seasonal trends are computed at summer (down to −0.010 yr⁻¹) and spring (down to −0.006 yr⁻¹), respectively. The effect of DOD on the total aerosol optical depth (AOD) change is determined by calculating the DOD-to-AOD trend ratio. Over the Sahara the median ratio values range from 0.83 to 0.95, whereas in other dust-affected areas (Arabian Peninsula, southern Mediterranean, Thar and Gobi deserts) the ratio value is approximately 0.6. In addition, a comprehensive analysis of the factors affecting the sign, the magnitude and the statistical significance of the calculated trends is conducted. Firstly, the implications of the implementation of the geometric mean instead of the arithmetic mean for trend calculations are discussed, revealing that the arithmetic-based trends tend to overestimate compared to the geometric-based trends over both land and ocean. Secondly, an analysis interpreting the differences in trend calculations under different spatial resolutions (fine and coarse) and time intervals is conducted.

1 Introduction

Dust particles emitted from natural or anthropogenic sources are a major contributor to the atmospheric aerosol burden in terms of mass (Zender et al., 2004; Textor et al., 2006; Kok et al., 2017). Among aerosol properties, aerosol optical depth (AOD) adequately describes aerosol load, in optical terms, corresponding to the entire atmospheric column. The proportion of AOD attributed to dust particles consists of dust optical depth (DOD). The spatiotemporal patterns of mineral particles are determined by the dust life cycle components, characterized by a pronounced heterogeneity (Mahowald et al., 2014). The main natural dust sources are located in the Northern Hemisphere (Goudie and Middleton, 2006), with the Sahara region being the most dominant one (Prospero et al., 2002; Goudie and Middleton, 2006; Rajot et al., 2008; Alizadeh-Choobari et al., 2014a). Other active source areas of mineral particles are situated in the Middle East and the region stretching from Mesopotamia to the Oman coasts on the southern Arabian Peninsula (Prospero et al., 2002; Ginoux et al., 2012), in southwestern Asia and the Sistan Basin (Iran–Pakistan–Afghanistan) (Alizadeh-Choobari et al., 2014b; Rashki et al., 2015), in central Asia across the Karakum (Turkmenistan–Uzbekistan) and Kyzylkum deserts (southeast of the Aral Sea in Uzbekistan) (Elguindi et al., 2016), in East Asia with the Taklamakan (Tarim Basin in northwestern China) and Gobi (northern China–southern Mongolia) deserts (Ginoux et al., 2012), and in North America with the Black Rock and Smoke deserts, the Great Salt Lake Desert, and the Chihuahuan and Sononan deserts (Ginoux et al., 2012).

Mineral dust aerosols are uplifted, accumulated into the atmosphere and transported over enormous distances (up to some thousands of kilometers) from their sources (Goudie and Middleton, 2006) driven by the prevailing winds. Schepanski et al. (2018) reported that the transport distance of dust particles is strongly related to their residence time, which is analogous to the dust lifetime, dust layer altitude, atmospheric circulation pattern, buoyancy and gravitational forces. Van der Does et al. (2018) also determined that strong winds, turbulence, electrostatic forces developed by dust particle charging, and thunderstorms or tropical cyclones may potentially enhance the residence time of dust aerosols in the atmosphere. On a seasonal basis, dust particles can be transported from North Africa towards the Atlantic Ocean, reaching the Caribbean, Central America, the southern United States (in boreal summer) and South America (in spring and winter) (Griffin et al., 2002; Prospero and Lamb, 2003; Kalashnikova et al., 2008; Huang et al., 2010; Tsamalis et al., 2013; Prospero and Mayol-Bracero, 2013). Additionally, Saharan dust is advected towards the Mediterranean and Europe (Mona et al., 2006, 2012; Papayannis et al., 2008; Basart et al., 2009; Schepanski et al., 2018; Gkikas et al., 2015, 2016; Logothetis et al., 2020, 2021).

During the last decades, numerous studies have been conducted using observations from various satellite sensors. Prospero et al. (2002) and Ginoux et al. (2012) identified global dust sources by relying on the Total Ozone Mapping Spectrometer (TOMS; Torres et al., 2002) and Moderate Resolution Imaging Spectroradiometer (MODIS; Remer et al., 2008), respectively. More specifically, the studies of Prospero et al. (2002) and Ginoux et al. (2012) were based on the frequency of occurrence (FoO) of the TOMS absorbing aerosol index (AAI) and MODIS-based DOD, respectively, exceeding defined thresholds. In addition, Ginoux et al. (2012) associated the dust frequency with three clusters, such as hydrologic and non-hydrologic natural or anthropogenic, in order to distinguish the dust origin. Similarly, at a regional scale, Schepanski et al. (2012) implemented a comprehensive analysis o the potential differences of Saharan dust active sources within the intercomparison of aerosol property observations derived from MODIS, as well as the Meteosat Second Generation (MSG) and Ozone Monitoring Instrument (OMI). Voss and Evan (2020) presented a global DOD climatology over both land and ocean using MODIS (Aqua and Terra) from 2001 to 2018 and the Advanced Very High Resolution Radiometer (AVHRR) over ocean from 1981 to 2018. Similarly, Clarisse et al. (2019) derived a global seasonal DOD climatology relying on Infrared Atmospheric Sounding Interferometer (IASI) retrievals during 2008–2017. Yu et al. (2019) derived DOD using MODIS, IASI and the Multiangle Imaging Spectroradiometer (MISR); in conjunction with dust vertical profiles from the Cloud–Aerosol Lidar with Orthogonal Polarization (CALIOP) (Shikwambana and Sivakumar, 2018), they investigated the dust deposition and loss frequency across the tropical Atlantic Ocean on a seasonal basis.

The investigation of dust load variation at interannual timescales is quite critical for assessing the associated impacts on climate and the response of these tendencies to environmental factors. Since the majority of remote sensing instruments provide an AOD product, numerous studies on a global scale are focused on the estimation of AOD temporal trends, which are not always representative of DOD, being mixed with other aerosol types (Zhang and Reid, 2010; de Meij et al., 2012; Hsu et al., 2012; Yoon et al., 2014; Pozzer et al., 2015; Alfaro-Contreras et al., 2017; Zhao et al., 2017; Che et al., 2019) and on regional scales (Guo et al., 2011; Li, 2014; Klingmüller et al., 2016; Floutsi et al., 2016; Dahutia et al., 2017; Hu et al., 2018; Zhang et al., 2018). Limited satellite studies are dedicated to the estimation of DOD temporal trends due deficiencies in accurately quantifying the portion of AOD attributed to DOD. Prior studies have investigated the interannual patterns of DOD over the “dusty” regions of the planet. The dust load increased across the Sahara based on a MODIS/Aqua-derived DOD dataset during 2003–2018 (Voss and Evan, 2020). Notaro et al. (2015) detected a regime shift in dust activity between 1998–2005 (inactive dust period) and 2007–2013 (active dust period) across the Arabian Peninsula, which is attributed to the prolonged drought along the Fertile Crescent. Through the synergy of MISR DODs and back trajectories, they revealed that the positive DOD anomalies (increased dust burden) are strongly connected with dust advection from the Fertile Crescent towards the Arabian Peninsula. These findings are consistent with the strong positive AOD (Klingmüller et al., 2016) and DOD (Voss and Evan, 2020) trends reported in the area. Voss and Evan (2020) found a reduction of dust load across the northern African coasts over the period 2001–2018 based on a MODIS/Terra DOD dataset. Declining DOD trends have also been reported in central Asia by Xi and Sokolik (2015), who analyzed MODIS and the Sea-viewing Wide Field-of-view Sensor (SeaWiFS) DODs for a 15-year period (2000–2014). The DOD trend sign has also abruptly changed from positive (1999–2009) to negative (2010–2016) over East Asia and the North Pacific Ocean in springtime based on the Modern-Era Retrospective Analysis for Research and Applications Version 2 (MERRA-2) (Gelaro et al., 2017) measurements (Guo et al., 2019). Across South Asia, a negative shift in DOD interannual variation is recorded during the pre-monsoon season between 2008–2012 and 2013–2017 based on CALIOP observations (Lakshmi et al., 2019). In the southern sector of the Gobi Desert, declining DOD trends are observed from MODIS and CALIOP DOD datasets during 2007–2019 (Song et al., 2021).

This study's main objective is to investigate dust temporal variations at global, regional and seasonal scales using the new MIDAS DOD product from 2003–2017. A few aspects regarding the innovative points of this research are highlighted below to support the scientific contribution to the relevant research field. First, in contrast to the existing studies, this trend analysis relies on fine-spatial-resolution data, making it feasible to depict the spatial patterns of DOD variations in detail. Such information can be critical for interpreting perturbations of the radiation fields, environmental impacts and health effects attributed to dust particles. One more advantage of the high-resolution DOD analysis is the flexibility of the final grid size selection depending on data availability, which is a critical aspect when satellite observations are used. MIDAS data can be easily upscaled at coarser spatial resolutions in order to match spaceborne observations, which have been commonly used in other trend analyses available in the literature (Hsu et al., 2012; Yoon et al., 2014; Notaro et al., 2015; Pozzer et al., 2015; Klingmüller et al., 2016; Alfaro-Contreras et al., 2017; Che et al., 2019; Guo et al., 2019; Voss and Evan, 2020; Song et al., 2021). In addition, fine-spatial-resolution data ensure a more realistic collocation with ground-based measurements for validating the obtained DOD trends. Second, quite a few studies concentrate on pure DOD (Xi and Sokolik, 2015; Guo et al., 2019; Lakshmi et al., 2019; Voss and Evan, 2020; Song et al., 2021) rather than AOD to analyze the trends of mineral particle load. Even though the consideration of the latter parameter is quite reasonable across deserts, its representativeness over downwind areas is questionable due to the co-existence of other aerosol types. Such types can also affect the DOD trend uncertainty. In MIDAS, this issue is addressed by the adjustment of MODIS AOD to DOD in terms of the MERRA-2 dust fraction, while in other studies, aerosol size and natural optical properties, the quality of which above land is downgraded, are used in parallel. Third, taking advantage of MIDAS providing DOD and quality-assured AOD, their trends are discussed jointly for assessing the contribution of dust burden temporal variations to those of the total aerosol load. It should be mentioned that this is the first study assessing the effect of DOD on total AOD trends across the major desert dust areas of the planet, highlighting the potential role of desert dust particles in past, present and future AOD trend studies. Fourth, there is an investigation of the potential impact on trend magnitude, sign and statistical significance when different DOD aggregations (i.e., arithmetic mean vs. geometric mean) are considered among various spatial and temporal scales. Fifth, the DOD interannual variations are discussed for the entire study period on a seasonal basis and sub-periods for detecting alternations in DOD trends within the period of interest.

The current paper is organized as follows. Section 2 describes the (i) MIDAS (Sect. 2.1) and (ii) AERONET (Sect. 2.2) datasets as well as the trend detection methodology (Sect. 2.3). The Results section (Sect. 3) is divided into three sub-sections analyzing (i) the global AOD and DOD tendencies (Sect. 3.1), along with the sensitivity analyses between arithmetic and geometric means (Sect. 3.1.1), fine ($\mathrm{0.1}{}^{\circ }\times \mathrm{0.1}{}^{\circ }$) and coarse ($\mathrm{1}{}^{\circ }\times \mathrm{1}{}^{\circ }$) spatial resolutions (Sect. 3.1.2), MIDAS DOD and AERONET retrievals (Sect. 3.1.3), filtered and non-filtered data trend calculations (Sect. 3.1.4), (ii) global dust temporal trends on a seasonal basis (Sect. 3.2), and (iii) DOD temporal tendencies in specific regions (Sect. 3.3). Finally, a discussion focusing on the main findings of this study is presented in the “Summary and conclusions” section (Sect. 4). The current study represents a practical implementation of the MIDAS dataset and aims to demonstrate its feasibility for the estimation of dust load variation at various temporal and spatial scales.

2 Data and methods

2.1 ModIs Dust AeroSol (MIDAS) dataset

The MIDAS dataset (Gkikas et al., 2021a) provides columnar DOD at 550 nm on a daily basis over the 15-year period spanning from 2003 to 2017 at a global scale and fine spatial resolution ($\mathrm{0.1}{}^{\circ }\times \mathrm{0.1}{}^{\circ }$). Its development has relied on the synergistic implementation of quality-filtered AOD retrievals from MODIS/Aqua (Level 2; Collection 6.1) and MERRA-2 dust fraction (MDF), both reported at 550 nm. More specifically, the multiplication of MODIS/Aqua AOD with MDF provides the MIDAS DOD on the MODIS native grid, which is converted to an equidistant lat–long projection. In order to justify the reliability of MDF, it has been evaluated against the corresponding portion provided by the LIVAS database (Amiridis et al., 2013, 2015). Based on the aforementioned assessment analysis, an adequate representation of the MERRA-2 dust fraction, in optical terms, has been revealed over the main dust sources and the outflow regions, in contrast to areas where dust presence is weak. MIDAS implements the MODIS aerosol products in which AOD is derived using the Dark Target (DT) and Deep Blue (DB) algorithms. Globally, these products have been thoroughly validated against AERONET, providing very good agreement either over oceanic areas for DT (R=0.880; MAE =0.055; RMSE =0.083) or over land for DT (R=0.920; MAE =0.066; RMSE =0.107) and DB (R=0.904; MAE =0.062; RMSE =0.107) (Wei et al., 2019b). Therefore, the combination of highly accurate MODIS AODs (Wei et al., 2019b) and quite reliable MDF results in a trustworthy MIDAS DOD product, as has been justified via its evaluation against AERONET dust-like AOD and its intercomparison versus DOD derived by LIVAS and MERRA-2. For the former analysis, ground-based AODs have been treated appropriately in order to resemble DOD as much as possible, assuming that the contribution of fine mineral particles is negligible and trying to minimize the contribution of non-dust aerosol species to the columnar aerosol load. Under these assumptions, the evaluation metrics at both global and station level reveal quite a high level of agreement between the two datasets. At global scale, there is a high level of agreement between MIDAS and AERONET DODs as indicated by the high correlation (∼0.9) and the low positive bias (0.004 or 2.7 %). Across the “dust belt”, the correlation coefficients can reach up to 0.98 at station level, whereas positive biases (mostly lower than 0.06) are found. Outside this zone, the correlation reduces, and the biases of similar magnitude switch to negative. Likewise, considerable consistency is evident among MIDAS, LIVAS and MERRA-2 DODs at global and hemispherical scales, despite the different approaches applied for the DOD derivation, whereas the intercomparison results are regionally dependent. Summarizing, the reliability of the MIDAS DOD justified in Gkikas et al. (2021a) strengthens its applicability for investigating the temporal trends of dust aerosol burden over long time periods and at various spatial scales.

2.2 AERONET

On a regional basis, the reliability of MIDAS DOD trends has been evaluated using ground-based measurements from AERONET (Holben et al., 1998). AERONET provides information for the spectral columnar AOD using direct spectral solar irradiance. In addition, the physical and optical properties of aerosols can be derived by applying the Dubovik and King (2000) inversion algorithm. DOD ground-based measurements can be approached using AERONET retrievals; however, an aerosol classification scheme is mandatory. Through the MIDAS evaluation procedure, the discrimination of dust load from non-dust aerosol species has been performed in terms of the Ångström wavelength exponent at 440–870 nm (AE_{440−870 nm}) and the difference between single-scattering albedo at 440 nm (SSA_440 nm) and 675 nm (SSA_675 nm). AE_{440−870 nm} has been used as an indicator of particles size (coarse or fine aerosols) in numerous aerosol classification studies (Basart et al., 2009; Mielonen et al., 2009; Lee et al., 2010; Giles et al., 2012; Hamill et al., 2016; Zheng et al., 2017; Che et al., 2018; Logothetis et al., 2020). In addition, the aerosol origin can be acquired using either the difference between SSA at different wavelengths or specific spectral SSA thresholds. In Gkikas et al. (2021a), pure dust conditions were successfully produced using a very restricted threshold limit of AE_440−870 nm (≤0.75) along with a positive difference of SSA_675 nm − SSA_440 nm. Despite the favorable dust conditions retrieved through the SSA retrievals, the data availability strongly decreases, causing difficulties in trend calculations. In the current study, the evaluation of MIDAS DOD trends is performed by using the AERONET coarse-mode AOD at 500 nm, derived from the spectral deconvolution algorithm (SDA) (O'Neill et al., 2001, 2003). The performance of the SDA algorithm, in terms of segregating the fine and coarse aerosol modes, has been justified against ground-based observations (Kaku et al., 2014). Only quality-assured data (cloud-screened), including pre-field and post-field calibrations (Level 2.0, L2) from AERONET version 3 (V3) (Giles et al., 2019), are used. It should be mentioned that cAOD does not profoundly represent the dust load due to the presence of sea-salt particles.

2.3 Temporal trend methodology

The spatiotemporal changes in dust particle burden over the period 2003–2017 are investigated by calculating the annual trends derived by the monthly MIDAS DODs. At each grid cell, the monthly DOD averages are calculated when 20 % (≥6 d) of daily data are available (Hsu et al., 2012). Subsequently, at the grid points with more than 60 months available (5 out of 15 complete years), linear trends are calculated by the implementation of the following equation:

$$\begin{array}{}\text{(1)}& Y_{\mathrm{t}}=\mathit{\mu }+S_{\mathrm{t}}+\mathit{\omega }{X}_{\mathrm{t}}+N_{\mathrm{t}},\end{array}$$

where Y_t represents the monthly averaged values, μ the offset term, S_t the seasonal term (i.e., long-term monthly value), ω the linear trend and N_t the residuals. The seasonality is removed by subtracting S_t from Y_t. The statistical significance of ω is derived according to Weatherhead et al. (1998). N_t follows a first-order autoregressive process (significant lag-1 autocorrelation),

$$\begin{array}{}\text{(2)}& N_{\mathrm{t}}=\mathit{\phi }{N}_{\mathrm{t}-\mathrm{1}}+{\mathit{\epsilon }}_{\mathrm{t}},\end{array}$$

with ε_t as the white noise and φ the lag-1 autocorrelation coefficient. The standard deviation of the trend can be expressed as

$$\begin{array}{}\text{(3)}& {\mathit{\sigma }}_{\mathit{\omega }}\approx {\displaystyle \frac{{\mathit{\sigma }}_N}{n^{\mathrm{3}/\mathrm{2}}}}\sqrt{{\displaystyle \frac{\mathrm{1}+\mathit{\phi }}{\mathrm{1}-\mathit{\phi }}}},\end{array}$$

where σ_N is the standard deviation of N_t and n is the number of complete years depending on the data availability at each grid cell without always considering the entire period (i.e., a constant value of 15 years). When $\mid \mathit{\omega }/{\mathit{\sigma }}_{\mathit{\omega }}\mid \mathit{>}\mathrm{2}$, significant temporal trends are considered at a 95 % confidence level. The methodology of Weatherhead et al. (1998) is commonly applied in numerous studies concerning the detection of temporal trends in AOD (Hsu et al., 2012; Babu et al., 2013; Li et al., 2014; Kumar et al., 2015, 2018; Pozzer et al., 2015; Adesina et al., 2016; Alfaro-Contreras et al., 2017; Zhang et al., 2018; Ningombam et al., 2019). Additionally, for comparison purposes with previous studies, AOD and DOD linear trends are also calculated at 1^∘ spatial resolution. The re-gridding procedure from fine to coarse spatial resolution is implemented following Levy et al. (2009) (upper branch in Fig. 5 of their publication). For the calculation of regional trends (Sect. 3.3), the same approach is adopted. First, daily spatial grids of 0.1^∘ are temporally averaged to create monthly data. Then, monthly grids with 1^∘ spatial resolutions are generated using a weighted aggregation of monthly fine grids. The weighting factors are defined in terms of latitude. More specifically, this weighting scheme considers the fraction between the area covered by each fine grid cell and the total available surface area within the coarse grid cell.

The appropriate selection of the statistical average metric (e.g., arithmetic or geometric mean) is reflected in the background probability distribution, which the raw data resemble. For instance, the vast majority of the studies focusing on AOD statistics have thoroughly considered AOD to follow a Gaussian distribution using the simple arithmetic mean for temporal and spatial aggregations. Nevertheless, the frequency distribution of AOD is more well fitted to the lognormal distribution (O'Neill et al., 2000). Sayer and Knobelspiesse (2019) determined that the calculation of the geometric instead of the arithmetic mean for obtaining temporal and spatial AOD trends may underestimate them compared to those reported in the literature. Here, in order to investigate the potential differences in AOD and DOD temporal trends, a sensitivity analysis using both the geometric and arithmetic mean is established. The MIDAS dataset includes negative DOD values introduced from the applied Dark Target algorithm of MODIS AOD retrievals. Since zero or negative arguments of logarithm cannot be defined, all these values are overwritten to 0.0001 as suggested by Sayer and Knobelspiesse (2019).

3 Results

This section is divided into three main parts. Section 3.1 describes the geographical distribution of AOD and DOD trends at global scale by simultaneously performing a sensitivity analysis of (i) the aggregation metric (geometric vs. arithmetic mean), (ii) the spatial resolution (fine vs. coarse), (iii) MIDAS trend evaluation (DOD vs. cAOD trends) and (iv) the applied temporal criteria (filters vs. no filters). In Sect. 3.2, focus is given to the seasonal DOD trends, whereas in Sect 3.3 emphasis is given to DOD trends in 12 regions of interest.

3.1 Global trends

3.1.1 Geometric vs. arithmetic mean

The overwhelming majority of the published AOD $/$ DOD trend analysis studies have relied on the arithmetic mean in order to produce coarser spatial (e.g., $\mathrm{1}{}^{\circ }\times \mathrm{1}{}^{\circ }$) and temporal (e.g., monthly values) resolutions. In this study, a sensitivity analysis of the aggregation method (i.e., arithmetic vs. geometric mean) has been performed by utilizing the daily coarse-spatial-resolution ($\mathrm{1}{}^{\circ }\times \mathrm{1}{}^{\circ }$) MIDAS DODs. Figure 1 depicts the frequency histogram of MIDAS DODs in log scale using all the available data over the period 2003–2017. It is apparent that the shape of the DOD distribution is close to a lognormal distribution (Fig. 1). The latter considers extreme dust episodes, which force the distribution curve to be right-skewed. The arithmetic mean is about 4 times higher than the geometric mean, highlighting the importance of which metric is more representative of the population and how this can affect the results for specific applications (e.g., radiative forcing). When the geometric standard deviation is considered, 65.43 % of DOD values range between 0.002 and 0.027, while 81.04 % encompasses −0.033 to 0.099 for the arithmetic mean, indicating an overstating of the variability (Sayer and Knobelspiesse, 2019). The averages and the standard deviations are calculated using weighted aggregation expressions with the weighting factors in terms of latitude (see Sect. 2.3). An additional point, revealing that the selection of normal statistics can lead to a misrepresentation of the population for a positive-definite quantity (i.e., dust optical depth), is that the lower bound (−0.033) is negative, which does not have any physical meaning.

Figure 1 Frequency histogram of MIDAS DODs reported at $\mathrm{1}{}^{\circ }\times \mathrm{1}{}^{\circ }$ spatial resolution over the period 2003–2017. The computed geometric mean (green) and arithmetic mean (blue) are given along with ranges corresponding to the geometric mean $/$ geometric standard deviation and geometric mean × geometric standard deviation for the geometric mean and ± standard deviation for the arithmetic mean.

In agreement with Sayer and Knobelspiesse (2019), the calculated trends for AOD and DOD are consistent in terms of sign between fine-spatial-resolution (Fig. 2) and coarse-spatial-resolution (Fig. 3) data. However, in terms of magnitude this is not the case. Figures S1 and S2 depict the frequency histograms of the deviations between the arithmetic and geometric trends. When geometric AOD $/$ DOD averages are considered, the deseasonalized trends are suppressed by up to 91.77 % with respect to the corresponding levels obtained from the arithmetic means, regardless of the underlying surface type. The only exception is found for AOD at 1^∘ spatial resolution and across oceanic territories (Fig. S2c) where geometric trends overestimate arithmetic trends by 71.63 %.

Figure 2 Global maps of temporal trends (statistically significant at 95 % confidence level) at $\mathrm{1}{}^{\circ }\times \mathrm{1}{}^{\circ }$ spatial resolution calculated from the deseasonalized AOD (a) and (c) and DOD (b) as well as (d) monthly values during 2003–2017. The upper panels, (a) and (b), show the arithmetic-based trends, while panels (c) and (d) indicate the geometric-based trends. For the pixels with (1) non-significant trends, (2) the temporal criteria are not met and (3) non-available data are colored white (blank) in each map of the trend analysis. Neutral trends, ranging between −0.002 and 0.002 DOD yr⁻¹, are colored light yellow (or cream) (the central color of the color bar). Figure S3a, b and S4 present the differences between the aforementioned colors at fine and coarse spatial resolution as well as for the seasonal trends, respectively. Based on those figures, the pixels with no available data are colored gray, while those that fail to meet the availability criteria are colored greenish.

Figure 3 Same as Fig. 2 at $\mathrm{1}{}^{\circ }\times \mathrm{1}{}^{\circ }$ spatial resolution.

Across the dust-affected regions of the planet, the calculated AOD and DOD trends have revealed similar patterns. Over areas where the dust contribution to the total aerosol load is negligible, DOD trends are non-significant or neutral. For instance, strong positive AOD trends are depicted across India and the Bay of Bengal (Fig. 2a and c), whereas the recorded annual DOD tendencies are negligible (Fig. 2b and d). Similar findings are evident along the eastern coasts of the US and in the Gulf of Mexico. In the Mediterranean Basin, disparities are recorded between AOD and DOD trends. Negative AOD trends are shown in the entire region, with decreasing DOD trends confined to the southern areas near the North African coast (Fig. 2b and d). In addition, strong positive trends for AOD and DOD are revealed in the central Sahara (up to 0.026 yr⁻¹), across Mauritania–Algeria–Mali–Niger and the Arabian Peninsula. The highest positive tendencies are shown along the Oman–Saudi Arabia borders (up to 0.031 yr⁻¹). In contrast, decreasing AOD $/$ DOD tendencies are observed in the eastern (down to −0.017 yr⁻¹) and western (down to −0.019 yr⁻¹) Sahara, in the Bodélé Depression of the Chad Basin (north of Lake Chad), in the Gobi Desert (northern China–southern Mongolia), and in the Thar Desert (northwestern Indian subcontinent). Among the regions where declining tendencies are evident, the most negative ones are recorded in the Bodélé Depression (down to −0.025 yr⁻¹) and in the Thar Desert (down to −0.029 yr⁻¹). A comprehensive regional analysis including intercomparison with prior findings and the potential trend justification is discussed in Sect. 3.3.

3.1.2 Fine vs. coarse spatial resolution

The second sensitivity analysis aims to highlight differences of AOD $/$ DOD trends when the fine and coarse spatial resolutions of MIDAS data are contrasted. De Meij et al. (2012) have investigated the differences between the fine daily MODIS AOD (L2; at $\mathrm{1}{}^{\circ }\times \mathrm{1}{}^{\circ }$ spatial resolution) and the coarse monthly MODIS AOD (L3; at $\mathrm{1}{}^{\circ }\times \mathrm{1}{}^{\circ }$ spatial resolution) in trend calculations utilizing the Collection 5 (C005) retrievals. According to their study, good agreement was found between the L2 and L3 AOD trends over specific areas (i.e., central Mediterranean, northeastern America and East Asia). At first glance, the trend patterns reproduced by the fine (Fig. 2) and coarse (Fig. 3) MIDAS DODs are spatially consistent. Nevertheless, in terms of magnitude, the absolute values of DODs at coarser spatial resolution are lower in most areas with an evident signal (either positive or negative), such as the southern parts of the Arabian Peninsula (up to 0.014 yr⁻¹), the Bodélé Depression (down to −0.015 yr⁻¹) and the Thar Desert (down to −0.024 yr⁻¹). Coarser grid cells, in contrast to the finer spatial resolution, meet the data availability threshold (≥60 months) (Fig. S3) defined for the calculation of temporal trends because of their more extensive spatial coverage. MIDAS adequately meets the temporal criteria (Sect. 2.2) at both fine (Fig. S3a) and coarse (Fig. S3b) spatial resolutions, providing grid cells of long-term AOD $/$ DOD time series along with significant AOD (Figs. 2a and c, 3a and c) and DOD (Figs. 2b and d, 3b and d) tendencies. Trend analysis for the coarse grids yields a higher number of significant AOD tendencies globally (Fig. 3a and b). For instance, new and significant declining coarser AOD trends are observed in East Asia, particularly across Southeast Asia, the Yellow Sea, the Sea of Japan and the North Pacific Ocean. Similarly, new increasing AOD trends are reported in the southern Arabian Sea and the North Atlantic Ocean near the coast of Venezuela.

Voss and Evan (2020) generated two global DOD datasets using MODIS retrievals, combined with reanalysis data and AERONET inversion retrievals. They estimated the decadal DOD trends (see their Fig. 11a and b) based on MODIS/Terra (2001–2018) and MODIS/Aqua (2003–2018) data projected at an equal lat–long 1^∘ spatial resolution. In order to compare the findings in this study against Voss and Evan (2000), only arithmetic DOD trends are used at the same grid cell spatial resolution (Fig. 3b). In addition, since the MIDAS dataset relies on MODIS/Aqua retrievals (Sect. 2.1), only their DOD MODIS/Aqua dataset is used for comparison. Table 1 displays the ranges of DOD trends between the two datasets over 15 desert and downwind areas of the globe. The selection of the above regions is accomplished based on the areas in which statistically significant DOD trends are revealed in the Voss and Evan (2020) study. Over the Sahara, the Bodélé Depression and the Thar desert, identical significant trends in terms of magnitude and sign are recorded in both studies. Over the Arabian Peninsula, the trends calculated here are common in terms of sign but lower in terms of magnitude. More specifically, the current study reports declining DOD trends along the Mediterranean Basin, while Voss and Evan (2020) did not find any significant trends. On the other hand, they reported strong positive trends over the tropical Atlantic Ocean, the sub-Sahel, the northeastern Middle East, and the northeastern Caspian and Aral Sea, which is not the case here. Trend inconsistencies are revealed because of (a) the derivation algorithm of DOD, (b) the trend detection methodology, (c) the different study periods and (d) the temporal filtering criteria.

Table 1 Comparison of the computed DOD trends between the current study and Voss and Evan (2020) over 15 regions of interest.

3.1.3 Evaluation of MIDAS trends vs. AERONET

The validity of the computed MIDAS DOD trends is thoroughly analyzed in this section. For evaluation purposes, the coarse-mode AOD, hereafter referred as cAOD, retrieved through the AERONET SDA algorithm is applied. The two datasets are collocated by spatially averaging MIDAS DODs included in an area of 3×3 pixels around each AERONET site and also temporally combined using the same daily values. DOD and cAOD trends are calculated using the methodology described in Sect. 2.3. Figure 4a illustrates the geographical distribution of the AERONET stations (in total 41) satisfying the defined temporal criteria (see Sect. 2.3) and finally used in the assessment analysis. According to Fig. 4b, a good linear correlation (R=0.86) of the derived trends is observed accompanied by low bias (MBE $=-\mathrm{0.04}\times {\mathrm{10}}^{-\mathrm{2}}$). In order to quantify the performance of MIDAS in capturing the cAOD trends, the total accuracy (TA) of the correct trends (identical sign) is calculated. More specifically, TA is defined as the percentage of the stations with correct trend signs (the sum of B1 and A2 quadrant numbers in Fig. 4c divided by the total number of stations). MIDAS trends can capture the correct cAOD trend signs with a TA of 80.5 % (Fig. 4c). In addition, the percentages of each quadrant in Fig. 4c show the number of correctly detected stations per sign. The overwhelming majority of the AERONET stations encompassed negative trends (30 out of 41 stations). MIDAS trends capture 26 out of those 30 (86.7 %) stations denoted with negative cAOD trends. Lastly, there are only 11 sites with positive cAOD trends, and MIDAS proved able to detect them at 7 out of 11 (63.6 %) stations.

Figure 4 (a) Geographical location of the available AERONET stations for which the temporal trends of DOD and cAOD are calculated. Green and red dots refer to the stations with similar and different trend signs, respectively. (b) Scatter plot between MIDAS DOD and AERONET cAOD trends. The color bar indicates the total number of months. (c) Confusion matrix of possible trend signs among the two datasets.

Table S1 provides the point-to-point trend results for each collocated MIDAS–AERONET site. It should be mentioned that among the revealed stations, few indicate dust particles as the predominant aerosol type. Based on Logothetis et al. (2020), the station of Solar Village, Arabian Peninsula, is primarily a “dusty” site revealing significantly high percentages of coarse absorbing particles. Across this station, the trends between AERONET cAOD (0.0126 yr⁻¹) and MIDAS DOD (0.0107 yr⁻¹) are in good agreement in terms of magnitude, but both are non-significant at the 95 % confidence level. Similar findings are also observed in the “dusty” site of Tamanrasset, Algeria, showing increasing cAOD and DOD trends of 0.0062 yr⁻¹ and 0.0044 yr⁻¹, respectively. One of the advantages of MIDAS is that it provides dust aerosol burden information in downwind regions of the planet. For instance, most of the stations located in southern Europe receive dust particles transported towards the Mediterranean from the Saharan and Middle Eastern deserts, encompassing identical trend signs (Fig. 4a). Moreover, there are many downwind regions across different regions of the globe with similar trend magnitudes such as Ispra in Italy (cAOD $=-\mathrm{0.0013}$ yr⁻¹, DOD $=-\mathrm{0.0015}$ yr⁻¹), Kanpur in India (cAOD $=-$0.0027 yr⁻¹, DOD $=-$0.0029 yr⁻¹), SERC (cAOD $=-\mathrm{0.0023}$ yr⁻¹, DOD $=-\mathrm{0.0020}$ yr⁻¹) in the US, Shirahama in Japan (cAOD $=-\mathrm{0.0020}$ yr⁻¹, DOD $=-\mathrm{0.0012}$ yr⁻¹) and Xianghe in China (cAOD $=-\mathrm{0.0028}$ yr⁻¹, DOD $=-\mathrm{0.0018}$ yr⁻¹).

3.1.4 Filtering vs. non-filtering trends

The third sensitivity analysis in this section concerns the calculation of temporal trends using filtered (Fig. 2) and non-filtered (Fig. 5) data at fine spatial resolution. Here, the AOD $/$ DOD trends are calculated by applying two consecutive temporal filters (Sect. 2.3). More specifically, the two applied temporal filters include the (1) the calculation of the monthly averages when at least 6 d are available and (2) the trend calculation only in grid cells with more than 60 available months (5 years). The total number of grid points for trend calculation increases when the above two temporal filters are not applied. Overall, months with more than one daily measurement are retained for trend analysis, increasing the total data availability of the entire global grid from 36 % (filtered) to 83 % (non-filtered).

Figure 5 Global maps of unfiltered temporal trends (statistically significant at 95 % confidence level) at $\mathrm{0.1}{}^{\circ }\times \mathrm{0.1}{}^{\circ }$ spatial resolution calculated from the deseasonalized geometric (a) AOD and (b) DOD monthly values during 2003–2017.

According to Eq. (3), σ_ω (Fig. S5) and the statistical significance (Fig. S6) of the trend are controlled by σ_N, φ (Fig. S7) and n (Fig. S3). The differences in data availability between the unfiltered and filtered trend analysis are depicted in Fig. S8. Across the desert areas, the number of filtered months is adequately high and very close to the non-filtered case (Fig. 8). Thus, no trend differences in magnitude and sign are recorded. Over maritime and continental dust-affected areas, new statistically significant AOD and DOD trends are represented when the unfiltered data are employed. Firstly, the number of pixels with statistically significant DOD trends has significantly increased (Fig. 5b), but the majority of the new trends are mainly neutral and located over oceanic territories (yellow pixels). New positive DOD tendencies are observed over the tropical Atlantic and India, while new negative trends are recorded across southeastern China (Fig. 5b). Secondly, the significant AOD trend grid points are also strongly increased. New decreasing AOD trends are observed over the USA, China and the Philippine Sea (Fig. 5a). Especially across the region of southeastern China, the number of significant AOD trends has increased profoundly. In addition, new AOD increasing trends are recorded over the tropical Atlantic, North Pacific (west of Mexico), Arabian Sea, and the oceanic area between 30.0 and 60.0^∘ S. Despite the increase in monthly data availability, trend analysis without temporal filtering may lead to erroneous and non-representative results either for AOD or DOD. In order to investigate the last statement in detail, a sensitivity analysis is performed using a total of 76 AERONET stations (Fig. S9a). Through the evaluation analysis of MIDAS DOD trends (Sect. 3.1.3), the AERONET stations have been selected based on their data availability. More specifically, only the stations with at least 10 years of data have been retained for further analysis. The total number of sites decreased from 76 (here) to 41 due to temporal filtering (see Sect. 3.1.3). Based on Fig. S9b, the unfiltered temporal trends between MIDAS and AERONET document a profoundly lower linear correlation (R=0.51), while the data points are more sparse. When trends are calculated without any temporal filtering, the total accuracy of the correct trends decreases from 80.5 % (Fig. 4c) to 73.7 % (Fig. S9c).

In the following sections and considering all three implemented sensitivity analyses, only the geometric-based DOD and AOD trends (Sect. 3.1.1) at fine spatial resolution (Sect. 3.1.2) including the temporal filtering (Sect. 3.1.4) are discussed.

3.2 Seasonal trends

Dust aerosol burden is subjected to strong intra-annual and interannual variations with different cycles depending on the source or downwind region (Gkikas et al., 2021b). Here, the seasonal DOD tendencies at a global scale (Fig. 6) are calculated based on the methodology proposed by Hsu et al. (2012). The corresponding seasonal AOD trends are depicted in Fig. S10.

Figure 6 Seasonal geographical distributions of DOD temporal trends in (a) December–January–February (DJF), (b) March–April–May (MAM), (c) June–July–August (JJA) and (d) September–October–November (SON).

The detection of the statistical significance of the calculated trends based on Weatherhead et al. (1998) cannot be applied here due to the 9-month gap among the seasons. Therefore, an alternative approach is followed by calculating the seasonal trends using a simple linear regression model on the DOD anomalies and identifying the statistically significant trends based on the two-sided Student's t test. The null hypothesis of the t test assumes a non-significant temporal trend under a defined confidence level (here is 95 %). The total number of months for each season is displayed in Fig. S4. Only the grid points with more than 13 available months (13 from 45 total months) are retained. The performed analysis at global scale and on a seasonal basis highlights many regions with significant DOD trends (Figs. 2 and 6), which are used to define regional domains (Fig. 7).

Figure 7 Regions of interest for seasonal and regional analysis: northern Arabian Peninsula (NAP), southern Arabian Peninsula (SAP), central Sahara (CSA), Gobi Desert (GOB), western Sahara (WSA), Mediterranean (MED), eastern Sahara (ESA), Bodélé Depression (BOD), Thar Desert (THA), eastern tropical Atlantic (ETA), eastern Middle East (EME) and Taklamakan Desert (TAK). Solid rectangles indicate the regions included in the regional analysis (Sect. 3.3), while solid and dashed rectangles are for the seasonal analysis.

The Sahara, hosting some of the most active aeolian natural dust sources of the planet, is a region of great interest for investigating intra-seasonal DOD variations. In its central sector, increasing trends are recorded throughout the year (maximum positive value in JJA) (Table 2). In contrast, negative seasonal DOD trends appear in the majority of regions across North Africa (Fig. 6). The western and eastern parts of the Sahara present strong declining trends maximized during boreal summer (Table 2). From spring to autumn in the Bodélé Depression, substantial decreasing trends are recorded (Table 2). Over the period 2001–2012, dust emissions in the broader area of the Bodélé Depression were decreased in summertime, which was attributed to the increased rainfall caused by the positive trends of the Sahara heat lows (SHL), the warm phase of the Atlantic Multidecadal Oscillation (AMO), and the decreasing trends in terms of occurrence and intensity of nocturnal low-level jets (NLLJs) (Shi et al., 2021). Surface wind speed also affects dust emissions across North Africa (Evan et al., 2016). Surface wind speed and NLLJs are the principal drivers for the interannual variation of dust emissions across the western Sahara, while in summertime dust emissions decreased during 2001–2012 (Shi et al., 2021). The eastern sector of North Africa (northern Libya and Egypt) presents moderate negative trends maximized in winter (down to −0.014 yr⁻¹) and spring (down to −0.011 yr⁻¹). In the western sector (i.e., northern Algeria and Tunisia) strong significant declining trends are observed in summer (down to −0.035 yr⁻¹) (Fig. 6c). The dust sources near the North African coasts are strongly influenced by the surface wind speed, NLLJs, Harmattan surge and the tracks of Mediterranean depressions (Shi et al., 2021).

Table 2 Seasonal trends over 12 regions of interest. Each column includes the range of DOD trend (yr⁻¹) and the median DOD trend. Non-available values are shown in blank cells. The domains of the regions are depicted in Fig. 7.

Numerous studies have highlighted Saharan transatlantic transport (Peyridieu et al., 2013; Alizadeh-Choobari et al., 2014c; Prospero et al., 2014; Gläser et al., 2015; Fréville et al., 2020; Gkikas et al., 2021b). Across the Gulf of Guinea and mid-Atlantic, relatively high DODs are documented in boreal winter, ranging from 0.1 to 0.2 (up to 0.6) (Gkikas et al., 2021b). According to Fig. 6a, strong positive DOD trends are shown over the Gulf of Guinea (up to 0.047 yr⁻¹). In this season, strong northeasterly winds (Harmattan) transport intense loads of Saharan dust towards Nigeria and the Gulf of Guinea (Washington et al., 2006). However, the trend magnitude along the Gulf of Guinea and in the northern regions (from Ghana to Cameroon) seems to be unreliable due to the reduced performance of MDF, surface contamination in the received radiances by MODIS, and the co-existence of aerosols and clouds (Gkikas et al., 2021b). In order to investigate the reliability of MIDAS DOD in the above region, the seasonal trend analysis is implemented at the AERONET station of Ilorin, Nigeria. It is worth mentioning that this station has not been included in the evaluation analysis in Sect. 3.1.3 due to the low number of months (N=48<60). However, 32 out of 48 months are in boreal winter, enabling the application of the methodology during this season. Both datasets represented increasing and statistically significant cAOD (0.011 yr⁻¹, $|\mathit{\omega }/{\mathit{\sigma }}_{\mathit{\omega }}|=\mathrm{2.43}$) and DOD trends (0.021 yr⁻¹, $|\mathit{\omega }/{\mathit{\sigma }}_{\mathit{\omega }}|=\mathrm{3.74}$). As expected, MIDAS revealed higher-magnitude trends than cAOD (approximately 2 times higher). In the case of AOD (Fig. S10a), positive trends identical to DOD are shown in winter. Interestingly, substantial negative DOD trends (down to −0.045 yr⁻¹) are recorded during springtime in an area encompassing parts of Guinea, Sierra Leone and the North Atlantic Ocean (Fig. 6b), while positive tendencies are documented over the tropical Atlantic Ocean in DJF (Table 2). Over the region extending from the North Atlantic Ocean to the eastern Caribbean Sea (lat: 10.0–18.0^∘ N; long: 70.0–45.0^∘ W), moderate positive trends are documented predominantly in MAM (up to 0.010 yr⁻¹) and JJA (up to 0.008 yr⁻¹). Dust transport from North Africa along the Atlantic Ocean reaches the Caribbean Sea, while summer DODs are intertwined with tropical Atlantic cyclone activity. More specifically, Caribbean DOD during summer is negatively correlated with Atlantic accumulated cyclone energy and the Atlantic Meridional Mode index (Xian et al., 2020). In addition, dust aerosol burden originating from the Sahara is directly related to tropical cyclones (TCs) over the tropical Atlantic Ocean, but their relationship is not unambiguously clarified (Evan et al., 2006). Based on model simulations, Pan et al. (2018) investigated the role of dust aerosols for the genesis of TCs over the tropical Atlantic basin by comparison of dust and non-dust model simulations. They revealed that dust amount is directly related to TCs, revealing a bimodal pattern. More specifically, favorable conditions for TC formation are documented by increasing the midlevel moisture and decreasing the vertical wind shear, while unfavorable conditions are documented by decreasing the low-level vorticity and potential intensity.

Increasing tendencies throughout the year appear in the Middle East (Table 2) with few exceptions, particularly in the northern Arabian Peninsula (Fig. 6). More specifically, the southwestern sector of Oman presents the highest increasing trends for all seasons (up to 0.026 yr⁻¹). Similarly, high positive trends in MAM (up to 0.015 yr⁻¹) and JJA (up to 0.018 yr⁻¹) are documented over the western part of Saudi Arabia. Negative tendencies are observed for all seasons at the north of the Oroug Bani M'aradh Wildfire Sanctuary (southern Saudi Arabia) including strongly negative DOD trends in summer (down to −0.020 yr⁻¹) and spring (down to −0.019 yr⁻¹). In the northern part of the Arabian Peninsula, positive trends are detected predominantly in JJA and SON. Dust activity across the Arabian Peninsula is strongly influenced by the intensity of the northwesterly Shamal winds, favored by the low precipitation amounts during summer (Yu et al., 2015). The temporally extended drought (Notaro et al., 2015) along with the cooling of the tropical Indian Ocean and the Mediterranean Sea temperatures (Yu et al., 2015), leading to an enhancement of the Shamal winds, could regulate the summer DOD trends across the Arabian Peninsula. During springtime, La Niña events constitute the principal drivers for dust activity by reducing the rainfall amounts over the Rub' al Khali Desert, which is one of the most active dust sources across the Arabian Peninsula (Yu et al., 2015). Strong positive trends (up to 0.026 yr⁻¹) are encountered in MAM and JJA over Iraq, while significant increasing trends are recorded for all seasons across eastern Iran, with the most evident ones in spring (up to 0.020 yr⁻¹). However, a hotspot of strong declining trends exists in the southeastern area of Iran (34.5^∘ N, 54.5^∘ E), with the most negative values in JJA (down to −0.029 yr⁻¹) and SON (down to −0.025 yr⁻¹). Moderate negative DOD trends are documented during the summertime (down to −0.01 yr⁻¹) across the Alboran Sea (western Mediterranean). The dust aerosol burden is also decreased in the eastern part of the Mediterranean Sea during spring (down to −0.006 yr⁻¹, from the Libya and Egyptian coasts to the Aegean Sea) and autumn (down to −0.005 yr⁻¹, across the Gulf of Sidra).

Statistically significant positive DOD trends are detected across the intersection of Kazakhstan, Uzbekistan and Turkmenistan at the northeastern Caspian Sea shore. In all seasons, the DOD trends exceed 0.011 yr⁻¹, while the maximum trends are recorded in summer (up to 0.035 yr⁻¹) and spring (up to 0.019 yr⁻¹). These findings substantiate the positive decadal DOD trends (∼ 0.18 decade⁻¹) of Voss and Evan (2020) and could be attributed to the amount of drawdown ($\sim -\mathrm{6.72}$ cm yr⁻¹) in the Caspian Sea level during 1996–2015 (Chen et al., 2017). Central, South and East Asia constitute another regions of interest in which robust DOD trends are encountered (Fig. 6). The maximum negative values are depicted over the Thar Desert in JJA and MAM (Table 2). It must be highlighted that the maximum decreasing trends are detected during the high-dust season in the Thar Desert (Proestakis et al., 2018). The reduction of dust load during the pre-monsoon (MAM) could be attributed to the increase in rainfall and soil moisture levels, acting in favor of wet dust deposition and decreasing dust erosion (Pandey et al., 2017; Jin and Wang, 2018; Lakshmi et al., 2019). Moreover, reductions in dust emissions are recorded during summertime, which are strongly linked to soil moisture and wind speed (Shi et al., 2021). In northwestern China (central Asia) lies the Taklamakan Desert, where non-significant annual trends are documented (Fig. 2). However, significant seasonally negative DOD trends are observed (Fig. 6), particularly in summertime (Table 2). Additionally, over the Gobi Desert and East Asia, strongly negative DOD tendencies are documented mainly in spring (Fig. 6b), when the dust activity peaks (Proestakis et al., 2018). The negative DOD trends across the Gobi Desert could be attributed to reduced dust emissions caused by the decrease in surface wind speed recorded between 2010 and 2016 (Guo et al., 2019).

3.3 Regional trends

The regional DOD and AOD trends are calculated for nine specific regions of interest (Table 3 and Fig. 7, solid rectangles) as well as for the whole globe (GLB) and separately over land (GLB-L) and ocean (GLB-O) (Sect. 3.3.1). The full names of each region as well as the calculated regional DOD trends and their uncertainties are included in Table 3. The comparisons among the geometric vs. arithmetic aggregation method and coarse vs. fine spatial resolution are also investigated for the regional DOD (Fig. 8a) and AOD (Fig. 8b) trends. On a regional basis, the temporal trends between the two spatial resolutions are in very good agreement, corroborating de Meij et al. (2012). Since the statistical significance of the trends is strongly influenced by the number of years and the study period, the regional DOD trends are also computed for different time periods considering the systematic change in time period (number of years) and initial year (Fig. 9).

Table 3 Global and regional temporal DOD trends based on the MIDAS dataset. The trends with the ratio $\left|\mathit{\omega }/{\mathit{\sigma }}_{\mathit{\omega }}\right|$ higher than 2.0 are statistically significant at the 95 % confidence level. The star symbol corresponds to statistically significant regions under the 95 % confidence level. The domains of the regions are represented in Fig. 7.

Figure 8 Bar plots indicating the (a) DOD and (b) AOD regional temporal trends. The hatched bars represent regions with significant trends ($\left|\mathit{\omega }/{\mathit{\sigma }}_{\mathit{\omega }}\right|>\mathrm{2.0}$). The error bars denote the uncertainty of the trends based on Eq.(3) (Sect. 2.3).

The explanation of the temporal evolution of calculated trends shown in Fig. 9 is presented in the next section for individual regions.

Figure 9 Regional DOD trends at different time intervals, with at least a 5-year time series. The acronym full names are documented in Table 3. The YY' axis shows the number of trends included in the analysis, while XX' is the starting year. The hatched pixels represent regions with significant trends ($\left|\mathit{\omega }/{\mathit{\sigma }}_{\mathit{\omega }}\right|>\mathrm{2.0}$).

3.3.1 Global land and ocean

Small global DOD trends are recorded during the study period both over land and ocean (Table 3). However, distinguishable DOD trends are detected in specific regions. Prior studies have focused on satellite-based measurements detecting statistically significant AOD trends at global scale. Over oceanic areas, AOD trends based on MODIS Collection 6.0 are reported to be equal to 0.0050 decade⁻¹ and 0.0020 decade⁻¹ during 2000–2009 and 2000–2015, respectively (Alfaro-Contreras et al., 2017). SeaWiFS AOD retrievals recorded higher annual positive significant trends over ocean (0.00080 yr⁻¹) for a 13-year period (1998–2010) (Hsu et al., 2012). Recently, significant positive tendencies have been documented for GLB (0.00066 yr⁻¹) using L3 Collection 6.1 MODIS/Terra measurements spanning from 2001 to 2016 (Che et al., 2019). In this study, significant AOD trends of 0.00052 and 0.00051 yr⁻¹ are revealed over GLB-O and GLB, respectively (Fig. 8b), using geometric means at fine spatial resolution. The differences in trend magnitude among the studies are attributed to the different datasets, aggregation methods and temporal availability.

3.3.2 North Africa

Across North Africa, four sectors have been defined based on the sign of the DOD trends (see Fig. 7). The first one consists of the central Sahara (CSA) where increasing DOD trends are mainly recorded (up to 0.023 yr⁻¹) (Fig. 10c). Voss and Evan (2020) also reported similar DOD trends, in terms of sign and magnitude, based on a MODIS/Aqua dataset over the period 2003–2018. At regional scale, positive significant DOD trends of 0.0021 yr⁻¹ determine those of AOD as expected (regional DOD-to-AOD trend ratio: 0.84) due to the predominance of mineral particles in the area among other aerosol species (Fig. 10d). During different time frames, the sign of DOD trends remains mainly positive (Fig. 9), with intense (from 0.0044 to 0.0095 yr⁻¹) and significant results from 2011 onward. Nevertheless, the time window is very short to infer the amplitude of the calculated trends.

Figure 10 (a) Geographical boundaries, (b) annual DOD, (c) DOD geometric trends and (d) DOD-to-AOD trend ratios for the central Sahara.

In the westernmost section of the Sahara (WSA), in the majority of grid cells (∼ 73 %) decreasing DOD tendencies are recorded (down to −0.015 yr⁻¹), whereas positive trends (up to 0.009 yr⁻¹) are evident at scattered pixels (∼ 27 %) within the domain (Fig. 11c). Overall, the total load (Fig. 8a) and the dust burden (Fig. 8b) decreased during the study period, but the magnitude of this reduction is relatively low and not statistically significant. This behavior is consistent regardless of the spatial resolution or the approach for the calculation of regional values (i.e., arithmetic or geometric mean) (Fig. 8a and b). The same DOD trend pattern is also reflected using different time periods (Fig. 9).

Figure 11 Same as Fig. 10, but for the western Sahara.

In the eastern sector of the Sahara (ESA), statistically significant negative DOD trends (down to −0.023 yr⁻¹) (Fig. 12c) are revealed as are both regional DOD (Fig. 8a) and AOD (Fig. 8b) trends, thus indicating the predominant contribution of mineral particles to the total aerosol load. Reduction of the dust burden has also been recorded from a MODIS/Terra dataset according to Voss and Evan (2020). Dust particles potentially affect total AOD across ESA, indicating a regional DOD-to-AOD trend ratio of 1.06. Over the dust-affected areas of the planet, the DOD-to-AOD trend ratios range from negative (different trend sign) to values higher than unity. Since the dust burden modulates the total AOD over “dusty” regions, the ratio between DOD and AOD is expected to be around unity. Ratios higher than unity are expected as the non-dust AOD signal decreases or increases and the DOD signal shows a reciprocal pattern. According to Fig. 9, two significant outcomes can be extracted in ESA. The magnitude, in absolute terms, of the decreasing DOD trends increases (down to −0.0054 yr⁻¹; 2012–2017) and statistical significance is achieved for shorter time periods for increasing starting years.

Figure 12 Same as Fig. 10, but for the eastern Sahara.

Within the study period, the regional dust load decreased (−0.0055 yr⁻¹) in the Bodélé Depression (BOD), which is the most active aeolian dust source of the planet (Prospero et al., 2002; Washington et al., 2006; Todd et al., 2007; Gkikas et al., 2021b). DOD trends range from −0.021 to −0.003 yr⁻¹ (Fig. 13c), corroborating the findings of Voss and Evan (2020). For starting years between 2003 and 2010 the dust aerosol burden decreases (Fig. 9), with the most negative trend found during 2007–2014 (−0.015 yr⁻¹). As the starting year increases and the time length decreases, the magnitude of DOD trends becomes weaker. More specifically, the sign of DOD trends is shifted using 2011–2013 as starting years, but the results are not statistically significant. The latter indicates that DOD over the most active dust source of the planet has become more intense during the last years. In order to better understand the reasons for the decreasing DOD trends across BOD, we investigated the correlation between the monthly MIDAS DODs and two meteorological variables, such as the wind speed (WS) at 10 m and the accumulated precipitation (STP), along with the volumetric soil water (layer: 0–7 cm) (VSM), retrieved from ERA5 reanalysis products (Hersbach et al., 2018, 2020). Figure S11 displays the trend analysis of the aforementioned parameters over the area of BOD. The solid black box depicted in Fig. S11a refers to the area that is presented in Fig. 13. Based on Fig. S11a, except for the small area where strong negative trends are documented, there are also many scattered pixels of negative trends in the surrounding area (16–19^∘ N–15–20^∘ E). According to Fig. 11b–d, it is evident that DOD values are related to temporal variability of applied parameters. Through time, dust amounts are strongly related to meteorology and geophysical factors. For instance, the increase in STP and VSM levels increases wet dust deposition and decreases dust erosion, causing a reduction of the dust load. On the other hand, the increase in WS acting in favor of dust erosion thus increases the dust amount. Across this region, increasing STP and VSM trends are documented, showing that the reduction of DOD levels is strongly related to the increase in STP and VSM. It should be mentioned that across this area the levels of WS and VSM provide a seasonal pattern, recording relatively low levels through the seasons of the year. During summer, the Intertropical Convergence Zone (ITCZ) location moves to higher latitudes and rapidly affects the amounts of WS and VSM across BOD. The revealed summertime STP and VSM trends are higher in terms of magnitude. In addition, positive trends have also been documented for WS. Generally, the increase in WS acting in favor of dust erosion thus increases the dust amount. However, the level of increased WS is revealed to be relatively small and may not lead to changes in dust amount across BOD. To our knowledge, there are no observed changes in the surface reflectance that could influence the true level of MODIS AOD and provide suspicious tendencies, like those of Klingmüller et al. (2016) study in the area of the Aral Sea, where the substantial trends in terms of magnitude are revealed and proved to be related to the land cover changes caused by the drying of the Aral Sea.

The aerosol burden over the Sahara is modulated by meteorological parameters affecting dust emissions such as precipitation (Pozzer et al., 2015) and wind speed (Che et al., 2019). Based on model simulations, Pozzer et al. (2015) confirmed that positive AOD trends across the Sahara are linked to negative precipitation trends during 2001–2010. Decreasing precipitation levels reduces dust deposition and favors dust transport, leading to an increase in dust amounts. In addition, Che et al. (2019) reported that wind speed is the most critical feature for explaining the variance of AOD (by 29.8 %). During high-dust seasons, the interannual variation of dust emissions is affected by wind speed and NLLJs (in the southern Sahara dust sources) as well as by the Harmattan surge and Mediterranean depressions (in the northern Saharan dust sources) (Shi et al., 2021).

Figure 13 Same as Fig. 10, but for the Bodélé Depression.

3.3.3 Arabian Peninsula

The DOD tendencies over the Arabian Peninsula are presented separately for the northern (NAP) and southern (SAP; including Jordan, Iraq and Syria) sectors. In both regions, positive trends (non-significant, Fig. 8a) are computed, which are stronger in the southern parts of the Arabian Peninsula.

Across the Arabian Peninsula, the sign of the trend provides a bimodal pattern, revealing strong positive trends for time frames up to 2013 followed by negative trends onwards (Fig. 9). Our findings are consistent with the revealed regime shift regarding dust activity on the Arabian Peninsula, as discussed in Notaro et al. (2015). The increased dust activity is also reflected in MIDAS DODs between 2008 and 2012, recording increased and stable annual DOD values (∼0.2) compared to DOD values for the previous 2003–2007 period (∼0.15). Due to these sudden changes in dust amount tendencies (from positive to negative), the regional DOD trends are revealed to be non-significant within the entire study period (Fig. 8a). Klingmüller et al. (2016) documented increasing AOD trends during 2001–2012 followed by declining tendencies onwards, in agreement with our findings (Fig. 9). The predominant role of mineral particles in the total aerosol load across the Arabian Peninsula is also reflected by the moderate to high DOD $/$ AOD trend ratios (up to 0.65; Figs. 14d and 15d). Increases in both DOD and non-dust AOD explain the moderate-magnitude trend ratios.

Figure 14 Same as Fig. 10, but for the northern Arabian Peninsula.

Figure 15 Same as Fig. 10, but for the southern Arabian Peninsula.

In the last few decades, compelling interannual dust activity has been documented for the Arabian Peninsula (Notaro et al., 2015). Numerous studies examined the temporal variability of aerosol loads in the Middle East, showing strong ascending tendencies (de Meij et al., 2012; Hsu et al., 2012; Yoon et al., 2014; Pozzer et al., 2015; Klingmüller et al., 2016; Che et al., 2019; Wei et al., 2019a). Klingmüller et al. (2016) revealed that the positive AOD trends in the Middle East are associated with decreasing trends of the Ångström exponent (AE) and fine-mode fraction (FMF), indicating a shift towards a regime with a stronger presence of coarse particles. These findings are in agreement with the increasing DOD trends documented in this study during the same study period (Fig. 9). Moreover, Klingmüller et al. (2016) implemented a multivariate linear model for annual AOD in order to identify the linkage of AOD trends with critical parameters such as precipitation, surface soil moisture and surface wind speed. Soil moisture is the major controlling parameter in Saudi Arabia and Iraq, whereas precipitation dominates in Iran. For all regions, the addition of surface wind speed as an independent parameter upgraded the model performance. Moreover, Che et al. (2019) also used a multiple linear regression model to investigate the relationship of AOD with specific meteorological parameters. The most appropriate ones, in terms of statistical significance, are chosen to enhance the model performance. Across the Middle East, the major controlling meteorological parameters for AOD variance are sea level pressure (60.9 % of total AOD explained variation) and wind speed, highlighting the large impact of synoptic systems on dust burden over the area.

3.3.4 Mediterranean

The Mediterranean (MED) Basin is a region of great concern due to high interannual variability of aerosol loadings and types (Floutsi et al., 2016). In this study, MED presents significant DOD trends ranging from −0.009 to 0.006 yr⁻¹ (Fig. 16c). The regional analysis documents strong declining significant DOD (Fig. 8a) and AOD (Fig. 8b) tendencies across the MED Basin. Negative DOD trends are also revealed during different sub-periods (Fig. 9). The DOD-to-AOD trend ratio shows a latitudinal reduction moving from northern African coasts to the northern parts of the Mediterranean (Fig. 16d). Higher ratio values are documented in the southern Mediterranean–North African coast region (0.3–1.94, median 0.71) (lat: 30.0–38.0^∘ N; long: 6.0^∘ W–30.0^∘ E) compared to the northern Mediterranean (0.21–0.91, median 0.45) (lat: 38.0–45.0^∘ N; long: 1.0^∘ W–27.0^∘ E). Dust particles originating from North African and Middle Eastern deserts driven by low-pressure systems (cyclones) can be transported towards MED (Gkikas et al., 2015), providing relatively high AOD values in the southern parts. Mineral particles are recorded mainly in summer, spring and winter in the western, central and eastern MED, respectively (Floutsi et al., 2016; Gkikas et al., 2013, 2016, 2021b). Across the north sector of MED, lower AOD values are associated with higher FMF values due to the prevailing anthropogenic fine aerosols (Floutsi et al., 2016). The latter could also be observed from the negligible DOD trends there (Fig. 16c).

Figure 16 Same as Fig. 10, but for the Mediterranean Basin.

Earlier studies investigated the temporal AOD variability in the broader Mediterranean Basin, reporting declining tendencies for the last 2 decades (Papadimas et al., 2008; de Meij et al., 2012; Hsu et al., 2012; Yoon et al., 2014; Pozzer et al., 2015; Floutsi et al., 2016; Che et al., 2019). Across this region, Floutsi et al. (2016) reported significant decreasing AOD trends of −0.0030 yr⁻¹ over the period 2002–2014. Additionally, Nabat et al. (2013) reported decreasing DOD trends (−0.0045 yr⁻¹) across northern Africa. Both studies corroborate the findings of this study: the overwhelming majority of DOD trends are primarily slightly negative (Fig. 16c).

3.3.5 Thar and Gobi deserts

Across the western part of the Indo-Gangetic Plain in the northwest area of the Indian subcontinent, the Thar Desert (THA) is situated. The THA region depicts significant decreasing DOD trends (down to −0.017 yr⁻¹) (Fig. 17c) modulating the corresponding declining AOD tendencies, as indicated by the relatively high DOD $/$ AOD trend ratio (0.67; Fig. 17d). The statistical significance of DOD trends, which are always negative regardless of the sub-period, is strongly affected by the starting year and the length of the time interval in which the tendencies are computed (Fig. 9). More specifically, negative and non-statistically significant DOD trends are observed for all time periods beginning from 2005. During the last 2 decades, strong negative temporal trends have been recorded for the OMI ultraviolet aerosol index (Hammer et al., 2018) along with AOD (Che et al., 2019) and DOD (Voss and Evan, 2020) across THA. The reduction of dust abundance over THA is mainly attributed to the increase in rainfall and soil moisture enhancing wet dust deposition and reducing dust erosion during the pre-monsoon (Pandey et al., 2017; Jin and Wang, 2018).

Figure 17 Same as Fig. 10, but for the Thar Desert.

The Gobi Desert (GOB) resides between the northern part of China and the southern sector of Mongolia (East Asia). In this region, statistically significant DOD trends ranging from −0.011 to 0.004 yr⁻¹ are reported (Fig. 18c), in agreement (both in sign and magnitude) with those reported by Che et al. (2019) and Voss and Evan (2020) for AOD and DOD, respectively. At a regional level, slightly negative DOD (−0.00071 yr⁻¹) and AOD (−0.0010 yr⁻¹) trends statistically significant at the 90 % confidence level ($\mid \mathit{\omega }/{\mathit{\sigma }}_{\mathit{\omega }}\mid \mathit{>}\mathrm{1.65}$) are computed (Fig. 8). Nevertheless, for the time period from 2005 onwards, the DOD trends are statistically significant at the 95 % confidence level (Fig. 9). These temporal DOD trends corroborate the results of Filonchyk et al. (2019), who computed AOD trends equal to −0.004 and −0.002 decade⁻¹ for MODIS/Terra and MISR measurements, respectively, over the period 2000–2017. In addition, GOB AOD trends are strongly influenced by the presence of dust particles, recording a moderate median DOD-to-AOD trend ratio value of 0.62 (Fig. 18d).

Figure 18 Same as Fig. 10, but for the Gobi Desert.

An et al. (2018) conducted a comprehensive analysis to investigate the potential factors driving the reduction of sand and dust storms, taking place mainly in springtime in East Asia between 2006 and 2017. In the aforementioned study, across GOB, the mean surface dust concentration declined by −31.71 µg m⁻³ yr⁻¹ (−12.24 %). These declining trends are regulated by changes in the surface conditions (e.g., vegetation coverage), which are strongly linked to precipitation, soil moisture, ambient temperature and human activities. The increasing trends in vegetation coverage are observed using the normalized difference vegetation index (NDVI) (0.0006 yr⁻¹) from MODIS and are linked with the ascending tendencies of precipitation (0.002 mm d ⁻¹ yr⁻¹) and volumetric soil moisture at 0–0.1 m of depth (0.316 yr⁻¹). All these factors converge towards the reduction of emitted amounts of dust aerosols. Over East Asia, the frequency and intensity of dust outbreaks have decreased during the springtime (An et al., 2018). Dust particles are primarily uplifted into the surrounding air across this region due to massive winds accompanied by cold air from northern latitudes. The decreasing frequency and intensity of dust outbreaks can be explained by the declining tendencies of the northern wind components documented during 2007–2016. The latter is associated with the decline of (1) polar vortex intensity, (2) north-to-south mean surface level pressure gradient and (3) meridional wind component magnitude.

4 Summary and conclusions

Airborne dust affects the global and regional climate via induced perturbation of the solar and terrestrial radiation fields. Therefore, the investigation of DOD temporal variations is crucial to assess its climatic role and the multifarious impacts within the Earth–atmosphere system as well as human health. The present study deals with the investigation of the annual and seasonal trends of AOD and DOD at both global and regional scales over the period 2003–2017, relying on the MIDAS fine-resolution dataset. Taking advantage of the MIDAS strong capabilities, the DOD trends have been analyzed (i) by considering different aggregation approaches (i.e., arithmetic and geometric means), (ii) at fine and coarse spatial resolutions, (iii) applying and not applying temporal filters (iv) at annual, seasonal and sub-period timescales (i.e., sliding window), and (v) in terms of their contribution to the corresponding tendencies of the total aerosol optical depth. Based on this holistic approach, a complete overview of the temporal variability of dust loads is provided, jointly addressing all the factors determining the sign, magnitude and statistical significance of the calculated trends.

The validity of the MIDAS DOD trends has been investigated with respect to AERONET cAOD data. This analysis considered MIDAS tendencies to be capable of following those of AERONET either for increasing (at 26 out of 30 stations, 86.7 %) or decreasing (at 7 out of 11 stations, 63.6 %) trends. Point-to-point investigation of the DOD–cAOD trends highlighted the accurate, in terms of magnitude and statistical significance, MIDAS DOD trends either across “dusty” sites or downwind areas.

The sensitivity analysis between coarse and fine spatial resolution revealed that the magnitude of AOD and DOD trends is lower for the coarse-spatial-resolution data. On a regional basis, the calculated AOD $/$ DOD trends revealed very good agreement in terms of trend magnitude for both spatial scales. In general, coarse-resolution data provide better spatial coverage; however, detailed information, in spatial terms, is provided when fine-resolution data are utilized. In addition, the sensitivity analysis of the different aggregation approaches revealed that the arithmetic-based trends are larger than those of the geometric-based ones (from 52.87 % to 91.77 %) all over the globe. The only exception was found only for AOD at coarse spatial resolution. Lastly, the sensitivity analysis of the temporal filtering in trend calculations indicated new regions (e.g., across oceanic territories) with unphysical unfiltered AOD trends. In addition, unfiltered trends led to poorer performance of MIDAS DOD trends against AERONET retrievals.

Based on our analysis, increasing DOD trends are obtained across the central Sahara and the Arabian Peninsula, whereas opposite tendencies are recorded over the eastern and western Sahara, the Thar and Gobi deserts, in the Bodélé Depression, and in the southern Mediterranean. On a seasonal basis, the maximum positive DOD trends are displayed in the central Sahara and Middle East during summertime. During spring, the strongest reductions of the dust burden have been revealed at the Bodélé Depression, whereas reverse tendencies have been recorded in the western and eastern Sahara and in the Thar and Taklamakan deserts during summer. Similarly, the most negative trends were observed in spring over the Gobi Desert. Positive trends across the area extending from the North Atlantic Ocean to the eastern Caribbean Sea are observed in spring.

At global scale and separately over continents and oceans, negligible and non-significant DOD trends are found. A similar analysis for nine specific regions of the planet has been performed, revealing evident trends that are not all statistically significant. More specifically, strong DOD trends were documented in the central (0.0021 yr⁻¹) and eastern (−0.0018 yr⁻¹) Sahara, the Bodélé Depression (−0.0055 yr⁻¹), the Mediterranean (−0.0011 yr⁻¹), and the Thar Desert (−0.0053 yr⁻¹). In contrast, non-significant regional DOD trends were depicted for the Arabian Peninsula and the Gobi Desert. At a regional scale, the calculated trends vary in terms of sign, magnitude and statistical significance depending on the sub-period, thus indicating a strong interannual variation of dust aerosol load within the period 2003–2017. In the majority of the regions of interest, also including downwind areas, the variations of the total aerosol load are driven by those of mineral particles.

The obtained findings from the current study highlight the applicability of the MIDAS dataset for detecting dust temporal variability across various spatial scales. The high spatiotemporal resolution of MIDAS provides the opportunity to complement and further expand the existing knowledge of this critical aspect that is not yet well covered in the field of dust research. Likewise, by comparing AOD and DOD tendencies, it is feasible to assess the role of mineral particles in variations of the total aerosol load. Our results could be incorporated into chemical models, either for assessing the various impacts of dust and non-dust particles or for evaluating trends based on numerical simulations, to further improve their calibration and forecast performance. Moreover, the findings obtained here can be used for the interpretation of the trends of radiation fluxes at the surface and at the top of the atmosphere obtained by observations or atmosphere–aerosol models (Saidou Chaibou et al., 2020). The recently published Sixth Assessment Report (AR6) of the Intergovernmental Panel on Climate Change (IPCC) states that the increase in the global dust load causes a reduction of the effective radiative forcing (ERF) by −0.25 Wm⁻² (Andrews et al., 2017; Forster et al., 2021). Therefore, fluctuations of the dust atmospheric load provide valuable information in order to understand the associated impacts on past, current and future climate. Finally, the current study focused mainly on the description of dust optical depth trends, also providing a preliminary interpretation of the obtained findings. Nevertheless, an in-depth analysis (a follow-up study) is ongoing to investigate the role of critical meteorological variables (e.g., wind) and of other determinant geophysical factors (e.g., soil moisture, vegetation, land coverage) in the configuration of DOD trends.