A Global Climatology of Dust Aerosols Based on Satellite Data: Spatial, Seasonal and Inter-Annual Patterns over the Period 2005–2019

A satellite-based algorithm is developed and used to determine the presence of dust aerosols on a global scale. The algorithm uses as input aerosol optical properties from the MOderate Resolution Imaging Spectroradiometer (MODIS)-Aqua Collection 6.1 and Ozone Monitoring Instrument (OMI)-Aura version v003 (OMAER-UV) datasets and identifies the existence of dust aerosols in the atmosphere by applying specific thresholds, which ensure the coarse size and the absorptivity of dust aerosols, on the input optical properties. The utilized aerosol optical properties are the multiwavelength aerosol optical depth (AOD), the Aerosol Absorption Index (AI) and the Ångström Exponent (a). The algorithm operates on a daily basis and at 1° × 1° latitude-longitude spatial resolution for the period 2005–2019 and computes the absolute and relative frequency of the occurrence of dust. The monthly and annual mean frequencies are calculated on a pixel level for each year of the study period, enabling the study of the seasonal as well as the inter-annual variation of dust aerosols’ occurrence all over the globe. Temporal averaging is also applied to the annual values in order to estimate the 15-year climatological mean values. Apart from temporal, a spatial averaging is also applied for the entire globe as well as for specific regions of interest, namely great global deserts and areas of desert dust export. According to the algorithm results, the highest frequencies of dust occurrence (up to 160 days/year) are primarily observed over the western part of North Africa (Sahara), and over the broader area of Bodélé, and secondarily over the Asian Taklamakan desert (140 days/year). For most of the study regions, the maximum frequencies appear in boreal spring and/or summer and the minimum ones in winter or autumn. A clear seasonality of global dust is revealed, with the lowest frequencies in November–December and the highest ones in June. Finally, an increasing trend of global dust frequency of occurrence from 2005 to 2019, equal to 56.2%, is also found. Such an increasing trend is observed over all study regions except for North Middle East, where a slight decreasing trend (−2.4%) is found.


Introduction
Dust aerosols (DA) is the most abundant aerosol type on Earth [1,2] contributing almost 30% of the total direct aerosol effect on a global basis [3]. The greater part of DA aerosols is injected into the atmosphere from the great deserts, especially from the Sahara [4,5], while a small part comes either from agricultural soils (<10%, remaining in the are used in the retrieval of atmospheric column ozone, sulfur dioxide, nitrogen dioxide, formaldehyde and other trace gas as well as clouds and aerosols. The OMI aerosol product contains information for the AI, the AOD, and the Single Scattering Albedo (SSA), which have been evaluated with AERONET observations. In the present study, the daily Level 3 AI data of the latest OMI edition OMAERUV [65] are used.
Both MODIS and OMI products used are available on a daily and 1 • × 1 • longitudelatitude resolution. The MODIS Aqua data were preferred to MODIS Terra, because Aqua and Aura are both A-train satellites, being in a polar orbit, crossing the equator within seconds or minutes of each other, thus ensuring the time collocation of the used products in the algorithm. Finally, the study period spans the 15-year period from 1 January 2005 to 31 December 2019, constrained by the simultaneous availability of both MODIS and OMI data.
It should be noted that Level 3 data, as those utilized in the present study, are appropriate for large-scale regional, and especially global-scale analyses [66,67], having the advantage of working on a regular latitude and longitude grid rather than on irregularly spaced pixels along the orbit track (Level 2 data, [68]). The only drawback of Level 3 data is that their resolution is too coarse to accurately detect small scale anthropogenic sources of dust [68], but the contribution of these sources is small compared to the contribution of the extensive physical sources of dust, namely the great desert areas.
Furthermore, it should be noted that time-averaging issues associated with Level 3 gridded products, induced by overlapping orbital tracks, may exist, only towards the poles (poleward of about 77 • N), in contrast to what happens at mid-latitudes, where the time averaging issue practically does not exist (https://modis-images.gsfc.nasa.gov/_docs/L3 _ATBD_C6.pdf). In particular, values from approximately 23 • N to 23 • S correspond to a single overpass, and they are not "daily mean values". However, even up to 60 • N and 60 • S, where the Global Dust Belt zone and the great world deserts are found, the 1-degree MODIS Level-3 products essentially do not involve a mixture of air masses, since they are based on either a single or two orbital swaths with a time lag of about 100 min, which corresponds to a short time for significant changes in dust to occur. It is only in higher latitudes, more poleward than 60 • , that the time averaging issue induced by multiple orbits exists. However, these regions contribute a negligible amount of the global dust load.

Methodology
The satellite algorithm uses as input the MODIS spectral AOD and the OMI AI data described in Section 2.1 and determines the presence or not of DA over a specific place on a specific day. It should be noted that the utilized satellite input data are columnar and that, usually, different aerosol types co-exist in the atmospheric column. In view of this, the presence of DA should be rather taken as a predominance of dust over other aerosol types in the atmospheric column, which is reflected in the satisfaction of specific criteria ensuring the coarse size and absorbing ability of predominant dust aerosol. In a first step, the algorithm uses the spectral AOD information to calculate the α, and subsequently, it applies specific threshold values on both α and AI, as well as on AOD, to identify the presence of DA (Figure 1). These thresholds are as follows: (i) α ≤ 0.4, (ii) AI ≥ 1 and (iii) AOD ≥ 0.1. These thresholds are representative for DA, and they were selected after sensitivity tests (Figures S1-S5) and the existing literature [69][70][71][72][73][74]. These sensitivity tests show that using lower/higher AI thresholds results in larger/smaller frequencies of DA both over land and ocean and an expansion/shrinking of areas covered by DA. On the other hand, although, in general, no substantial differences appear between the results obtained for α ≤ 0.4 and α ≤ 0.7 threshold values, slightly higher frequencies of DA are found above oceanic regions, i.e., Tropical and South Atlantic Ocean, Gulf of Guinea and East Pacific Ocean. The obtained results of the algorithm using the selected thresholds confirm the presence of coarse absorbing aerosols, which is the case for DA. The algorithm identifies the presence of DA on a daily basis and subsequently calculates the frequency of occurrence of DA on a monthly and annual basis. It should be noted that apart from the absolute frequency of occurrence (Nabs) of DA, given in days/month or days/year, the relative frequency of occurrence (Nrel) is also calculated, given as a percentage of days (in a month or year) characterized by the presence of DA with respect to the number of the days (of the month or year) for which the algorithm ran. The algorithm ran for each year of the study period, producing monthly and annual mean Nabs and Nrel values, and subsequently the corresponding average 15-year (climatological) values (monthly and annual) were calculated.
Apart from temporal averages, regional monthly and annual mean absolute and relative frequencies for the globe, but also for specific world areas, namely great deserts, were also produced. The absolute frequencies result from the sum of grids over which DA were detected, while the relative frequencies were computed by dividing the number of grids with DA detected over the area of interest (region or globe) by the total number of grids of the area for which the algorithm operated. Here, it should be noted that the denominator in this calculation can be a large number since it includes a significant number of grids over which DA either are not at all or are rarely detected. In fact, this results in a lower frequency of occurrence than what would it be if the averaging was made only over grids with DA, as done in some previous studies. The selected areas of special interest with relevance to DA are those where dust aerosols are systematically the predominant aerosol type in the atmospheric column. The Sahara is such a region, being the greatest world desert and contributing more than 50% of the global dust emissions [75]. The other selected regions of interest are smaller deserts, namely the Taklamakan [43] and Gobi [76] Asian deserts and the Middle East deserts [4,44], as well as land and oceanic areas undergoing frequent transport of DD. Such areas selected in this study are the Tropical Atlantic Ocean [77][78][79], the Gulf of Guinea [55,80] and the Mediterranean Basin [42,81]. In the present study, ten (10) regions with a high frequency of dust occurrence are discussed. The selected world study regions are displayed in Figure 2. It should be noted that the areas of the Sahara and Tropical Atlantic are divided into East and West Sahara due to their The algorithm identifies the presence of DA on a daily basis and subsequently calculates the frequency of occurrence of DA on a monthly and annual basis. It should be noted that apart from the absolute frequency of occurrence (N abs ) of DA, given in days/month or days/year, the relative frequency of occurrence (N rel ) is also calculated, given as a percentage of days (in a month or year) characterized by the presence of DA with respect to the number of the days (of the month or year) for which the algorithm ran. The algorithm ran for each year of the study period, producing monthly and annual mean N abs and N rel values, and subsequently the corresponding average 15-year (climatological) values (monthly and annual) were calculated.
Apart from temporal averages, regional monthly and annual mean absolute and relative frequencies for the globe, but also for specific world areas, namely great deserts, were also produced. The absolute frequencies result from the sum of grids over which DA were detected, while the relative frequencies were computed by dividing the number of grids with DA detected over the area of interest (region or globe) by the total number of grids of the area for which the algorithm operated. Here, it should be noted that the denominator in this calculation can be a large number since it includes a significant number of grids over which DA either are not at all or are rarely detected. In fact, this results in a lower frequency of occurrence than what would it be if the averaging was made only over grids with DA, as done in some previous studies. The selected areas of special interest with relevance to DA are those where dust aerosols are systematically the predominant aerosol type in the atmospheric column. The Sahara is such a region, being the greatest world desert and contributing more than 50% of the global dust emissions [75]. The other selected regions of interest are smaller deserts, namely the Taklamakan [43] and Gobi [76] Asian deserts and the Middle East deserts [4,44], as well as land and oceanic areas undergoing frequent transport of DD. Such areas selected in this study are the Tropical Atlantic Ocean [77][78][79], the Gulf of Guinea [55,80] and the Mediterranean Basin [42,81]. In the present study, ten (10) regions with a high frequency of dust occurrence are discussed. The selected world study regions are displayed in Figure 2. It should be noted that the areas of the Sahara and Tropical Atlantic are divided into East and West Sahara due to their different magnitude of dust loadings and spatial extent, as is discussed in Section 3. Similarly, Middle East area is divided into North and South Middle East.

Results and Discussion
The algorithm was applied to every day of each year of the study period of 2005-2019 and worked on a daily and 1° × 1° latitude-longitude grid-level basis. Subsequently, the monthly and annual mean global distribution of the absolute (Nabs) and the relative (Nrel) frequency of occurrence of DA were computed. The outputs were then averaged over the entire study period, yielding the 15-year averaged values, which are presented in this section.

Geographical Distribution
The global distribution of the absolute annual mean frequency of occurrence of DA, expressed in number of days/year, averaged for the entire period of 2005-2019 is displayed in Figure 3a. Figure 3b displays the corresponding distribution for the DA relative frequency of occurrence, with respect to the number of days in a year for which the algorithm operated. White areas in the figures correspond to cases (grids) for which the algorithm never operated due to missing input data (non-operating algorithm). The two distributions are quite similar, indicating that the spatial patterns of the absolute frequency are not dominated by data availability. In both distributions, it is clear that the greatest part of the globe has low frequencies of DA (bluish colors) up to 5 days/year or 0.1 (i.e., 10%). Higher frequencies are observed mainly in the "global dust belt". This belt includes the greatest world deserts, namely the Sahara, Middle East, Taklamakan, Thar and Gobi deserts, expanding to the west and including a significant part of the northern Tropical Atlantic Ocean. It should be noted that relatively small frequencies are observed off the coast of Angola, in the southeastern Atlantic Ocean, not exceeding about 10% or 20 days/year. These occurrences, which take place from July to September (as shown and discussed in Figure 4) do not really correspond to the predominant presence of dust over this oceanic area. As shown by our analysis ( Figure S6 and relevant discussion at the end of this subsection), they arise from the coexistence of sea-salt and biomass-burning aerosols exported from Angola. This coexistence of coarse scattering sea salt and fine-absorb-

Results and Discussion
The algorithm was applied to every day of each year of the study period of 2005-2019 and worked on a daily and 1 • × 1 • latitude-longitude grid-level basis. Subsequently, the monthly and annual mean global distribution of the absolute (N abs ) and the relative (N rel ) frequency of occurrence of DA were computed. The outputs were then averaged over the entire study period, yielding the 15-year averaged values, which are presented in this section.

Geographical Distribution
The global distribution of the absolute annual mean frequency of occurrence of DA, expressed in number of days/year, averaged for the entire period of 2005-2019 is displayed in Figure 3a. Figure 3b displays the corresponding distribution for the DA relative frequency of occurrence, with respect to the number of days in a year for which the algorithm operated. White areas in the figures correspond to cases (grids) for which the algorithm never operated due to missing input data (non-operating algorithm). The two distributions are quite similar, indicating that the spatial patterns of the absolute frequency are not dominated by data availability. In both distributions, it is clear that the greatest part of the globe has low frequencies of DA (bluish colors) up to 5 days/year or 0.1 (i.e., 10%). Higher frequencies are observed mainly in the "global dust belt". This belt includes the greatest world deserts, namely the Sahara, Middle East, Taklamakan, Thar and Gobi deserts, expanding to the west and including a significant part of the northern Tropical Atlantic Ocean. It should be noted that relatively small frequencies are observed off the coast of Angola, in the southeastern Atlantic Ocean, not exceeding about 10% or 20 days/year. These occurrences, which take place from July to September (as shown and discussed in Figure 4) do not really correspond to the predominant presence of dust over this oceanic area. As shown by our analysis ( Figure S6 and relevant discussion at the end of this subsection), they arise from the coexistence of sea-salt and biomass-burning aerosols exported from Angola. This coexistence of coarse scattering sea salt and fine-absorbing-biomass-burning aerosols, along with a weak presence of dust over this area, is misinterpreted by the present satellite algorithm, which uses columnar aerosol products as a presence of coarse-absorbing DA.
Remote Sens. 2020, 17, x. https://doi.org/10.3390/rsxxxxx www.mdpi.com/journal/rs ing-biomass-burning aerosols, along with a weak presence of dust over this area, is misinterpreted by the present satellite algorithm, which uses columnar aerosol products as a presence of coarse-absorbing DA.
(a) (b) Figure 3. Global distribution of: (a) absolute frequency of occurrence of dust aerosols (in number of days/year); (b) relative percentage frequency of occurrence of dust aerosols (with respect to the number of days per year with available satellite data and for which the satellite algorithm operated). The results are averaged over the period of 2005-2019. The white-shaded areas are those for which the algorithm did not operate due to missing input data (non-operating algorithm).
The highest frequencies of DΑ are observed over North Africa (Sahara Desert) where the values ranged from about 60 up to 120 days/year, while in some areas (Bodélé), they are as high as 160 days/year. High frequencies are also observed over the Asian deserts Thar, Taklamakan and Gobi, with values ranging from 70-80 (greenish colors) up to 150 days/year (reddish colors), with the highest frequencies mainly observed over the Taklamakan. Quite high DA frequencies (up to 90 days/year) were observed over the Middle East, especially over its southern part, namely the Arabian desert. The frequencies of DA are generally low (up to 20 days/year) over the oceanic areas, while greater frequencies are observed only over the Mediterranean Sea (up to 30 days/year), the Arabian Sea (up to 40 days/year), the Tropical Atlantic Ocean (up to 50 days/year) and the Gulf of Guinea (up to 30 days/year). The results are averaged over the period of 2005-2019. The white-shaded areas are those for which the algorithm did not operate due to missing input data (non-operating algorithm).
The highest frequencies of DA are observed over North Africa (Sahara Desert) where the values ranged from about 60 up to 120 days/year, while in some areas (Bodélé), they are as high as 160 days/year. High frequencies are also observed over the Asian deserts Thar, Taklamakan and Gobi, with values ranging from 70-80 (greenish colors) up to 150 days/year (reddish colors), with the highest frequencies mainly observed over the Taklamakan. Quite high DA frequencies (up to 90 days/year) were observed over the Middle East, especially over its southern part, namely the Arabian desert. The frequencies of DA are generally low (up to 20 days/year) over the oceanic areas, while greater frequencies are observed only over the Mediterranean Sea (up to 30 days/year), the Arabian Sea (up to 40 days/year), the Tropical Atlantic Ocean (up to 50 days/year) and the Gulf of Guinea (up to 30 days/year).
The seasonal variation of dust aerosols is shown in Figure 4, which displays the monthly frequencies of DA averaged over the 15-year study period. Seasonal patterns are apparent only over the world great deserts and their neighboring areas, namely over the ten (10) selected areas of interest in this study (see Section 2.2 and Figure 2). These seasonal patterns, along with the inter-annual variability of DA, are discussed in detail in Section 3.2, where the regional results are presented. Below, only a brief discussion of the The seasonal variation of dust aerosols is shown in Figure 4, which displays the monthly frequencies of DA averaged over the 15-year study period. Seasonal patterns are apparent only over the world great deserts and their neighboring areas, namely over the ten (10) selected areas of interest in this study (see Section 2.2. and Figure 2). These seasonal patterns, along with the inter-annual variability of DA, are discussed in detail in Section 3.2, where the regional results are presented. Below, only a brief discussion of the seasonal variation of globally distributed DA frequency of occurrence is presented, along with the seasonal and inter-annual variation of 15-year average global DA frequencies.     . Although this seasonality prevails globally, it is much more pronounced over the global dust belt zone and its surrounding areas. The global dust belt zone appears more or less clearly, being expanded or shrunk in every month of the year. A latitudinal shift of this zone is notable, which is driven by the corresponding latitudinal shift of the Intertropical Convergence Zone (ITCZ). Indeed, during winter, the dust sources in southern Sahara, Sahel and sub-Sahel regions (<20 • S) are well activated and combined with the blowing northeasterly (NE) Harmattan trade winds, resulting in the export of dust aerosol burden to the Gulf of Guinea [82][83][84]. On the other hand, during summer, the ITCZ and trade winds are shifted northwards, thus intensifying the Saharan dust sources and exporting DA to the northern Tropical Atlantic [4,82,85,86]. At the same time, the increased precipitation in the Sahel region, due to the African monsoons, decreases the airborne dust over this area [86]. During November-December-January, the largest part of the map is deep blue, indicating that the absolute frequency of dust occurrence is lower than 2 days/month. Higher frequencies of DA are observed only in January over some areas in Africa, and more specifically in the sub-Sahel region and the region north of the Gulf of Guinea, where the absolute frequency takes values up to 5-10 days/month, and especially over the Bodélé depression area, where much higher values (up to 15 days/month) are observed. The period from November to March is known as the dry or Harmattan period characterized by northeasterly trade winds over the area by high frequencies of dust occurrence [87] and dust optical depth [55,88]. In particular, Bodélé is one of the most intense dust sources in the world [4], where DA are observed all year round. Indeed, according to the algorithm results, the annual minimum frequency of DA occurrence over Bodélé, observed in November, is equal to 5-10 days/month. The DA frequencies gradually increase towards the end of boreal winter, and thus in February values up to 10 days/month are observed above the world deserts (e.g., in Northern Africa, Arabian Peninsula and Gobi).
During March-April-May, the frequencies of DA are further intensified and become more widespread. More specifically, the worldwide highest frequencies above the Sahara and Middle East deserts reach values up to 5-10 days/month in March, 7-12 days/month in April and 8-17 days/month in May. During the same period (March-April-May) the DA frequencies also increase over the Asian deserts (Thar, Taklamakan and Gobi) up to 7-15 days/month (in May). In fact, these spring frequencies are the annual maximum values over the Asian deserts, higher than the corresponding frequencies in summer, which are slightly weakened (being the maximum equal to 12 days/month). Spring maximum DA frequencies over the Asian deserts similar to those reported here (up to half days of the season) are also reported by , who used vertically resolved data from Cloud-Aerosol Lidar and Infrared Pathfinder Satellite Observations (CALIPSO) lidar measurements. During spring, a dust transport from Asian deserts, namely Taklamakan and Gobi, towards the Pacific Ocean is noted, as indicated by light blue colors, corresponding to about 3 days/month. Such a dust transport, even reaching the western coasts of North America, has also been reported previously in the literature [43,[89][90][91]. Dust loads and the corresponding shortwave Radiative Forcing (RF) over the West American coast were estimated by [92] using satellite data, i.e., MODIS and CERES, for three springs (namely 2004, 2005, 2006). According to their results, the mean DOD over the area is 0.1, causing a −5.51 W/m 2 shortwave RF. On the other hand, [93] provided a 10-year (2002-2011) climatology of dust aerosols over the same area, using ground-based data from the Interagency Monitoring of Protected Visual Environments (IMPROVE) network. According to their results, the climatological mean spring dust concentration was equal to 0.08-0.60 µgm −3 , while there is transport of DA from Asian deserts towards the Pacific Ocean. Ref. [94], using synergetic (ground-based, satellite and modeled) data, investigated two severe dust events that took place over Asian deserts and found that dust was transported to northern China, the East/Japan Sea, through to the North Pacific Ocean.
In June, the DA frequencies are even higher than in spring over North Africa (Sahara) and especially over its western part (values up to 17 days/month). These high frequencies of dust are attributed to the prevailing local pressure systems, the associated atmospheric circulation and, more specifically, the prevailing Sahara Heat Low (SHL) and west Sahara monsoon, which favor the uplift of dust. The higher DA frequencies over the western than the eastern parts of the Sahara can be explained by the easterly trade winds. Further details about this are given in the next section (Section 3.2.1). In summer, high DA frequencies are also observed over the Arabian Peninsula (up to 18 days/month) and East Sahara (up to 16 days/month), as well as over the Indo-Gangetic Plain (up to 11 days/month) and the Arabian Sea (up to 11 days/month). These results are in agreement with those reported in the existing literature [54,[95][96][97], which are based on other datasets, e.g., European Organization for the Exploitation of Meteorological Satellites (EUMETSAT), CALIPSO, KALPANA and Indian National Satellite-3 (INSAT) 3A as well as ground-based data (from NASA's AEROsol RObotic NETwork (AERONET) and Micro-Pulse LiDAR Network (MPLNET)) and refer to shorter study periods. According to the algorithm results, during summer, and especially in June, the export of African dust to the Tropical Atlantic, which takes place since January, but to a smaller extent and with smaller intensity, is also intensified, in agreement with existing literature [77,98]. During the next two summer months, i.e., in July and August, the high DA frequencies are displaced to northwestern Africa, while the dust export to the Tropical Atlantic Ocean is also shifted slightly northwards and weakened. Finally, in September, and gradually through to the end of the year, the DA frequencies significantly and gradually decreased all over the globe, and high values (up to 16 days/month) are observed only over the deserts of North Africa and Asia.
A remarkable feature in Figure 4 is the presence of DA in the southeastern Atlantic Ocean, off the coasts of Angola, in up to 3-16 days/month during the period from July to September, maximized in August (up to 19 days/month). This appearance of DA is observed over the same location, where absorbing carbonaceous aerosols (ACA) are observed and detected by a companion version of the present algorithm [99] using the same satellite data (MODIS and OMI). These ACA are associated with an effective westward transport of African biomass aerosols originating from extensive agricultural burning and wildfires over Angola and Congo [99]. These ACA aerosols are convectively uplifted to the free troposphere and subsequently transported westward (e.g., up to the Ascension Island 2000 km offshore continental Africa) by a well-established anticyclonic circulation induced by the semi-permanent sub-tropical anticyclone of the South Atlantic [100]. The coexistence of these ACA, which are fine and absorbing, above coarse and scattering sea-salt aerosols in the underlying atmospheric layers leads to a misidentification of predominant aerosols in the atmospheric column by the present satellite algorithm. This is verified by an analysis using MERRA-2 (Modern-Era Retrospective analysis for Research and Applications, Version 2; [101]) vertically resolved (at 50 atmospheric layers from the Earth's surface up to 6 mb) extinction coefficients for five specific (5) aerosol types, namely sulfate, dust, sea-salt, black carbon and organic carbon. The obtained long-term (2005-2019) monthly mean MERRA-2 vertical aerosol profiles ( Figure S6) show that indeed, from July to September, ACA, and more specifically organic carbon aerosols (red-colored curves), exist in the free troposphere (1-4 km), especially in its lowest heights (around 1000 m) overlying sea-salt aerosols (bluish curves) in the boundary layer (lower than 1000 and even more than 500 m). At the same time, the presence of DA (yellow curves) over this area and period of year is weak. It should be noted that this misidentification of DA is an inherent limitation not only of the present but of every satellite algorithm using columnar products, in cases of coexistence of sea-salt (coarse strongly scattering) and carbonaceous (fine strongly absorbing) aerosols in the atmospheric column. The only way to overcome this difficulty is to rely not on columnar but on vertically resolved aerosol products, but these are either spatially and temporally limited-for example, those from Cloud-Aerosol Lidar with Orthogonal Polarization, CALIOP onboard the CALIPSO (Cloud Aerosol Lidar and Infrared Pathfinder Satellite Observation) satellite data-or they are model driven-for example, the MERRA-2 data, which are based on the coupled GEOS-5 (Goddard Earth Observing System, [102,103]) and GOCART (Goddard Chemistry Aerosol Radiation and Transport model, [104,105]) models.

Intra-Annual and Inter-Annual Cycles
The intra-annual variability of the frequency of occurrence of dust aerosols on a global basis is depicted in Figure 5, where the total number of grid cells (1 • × 1 • ) over which DA have been observed, in terms of averaged values over the study period 2005-2019, is given. The maximum frequency of DA is noted in June, when dust is observed over almost 23,400 pixels, i.e., over the 3.6% of the planetary surface. Note that 23,400 is the average number of grid cells over which the algorithm detected DA in June (mean value of June's cumulative DA daily frequencies over the period 2005-2019). This number is divided by the average (over the period 2005-2019) number of grid cells of the globe for which the algorithm operated during June, i.e., 633,463 grids, resulting in the low-frequency value of 3.6%. The seasonal cycle of global dust is clear, with minimum frequencies in boreal winter, taking its minimum in November, when dust is observed over only 3181 pixels or less than 1% of the planet. The error bars in Figure 5 correspond to the year-to-year (during the 15-year study period) variability of the frequency of DA on a global basis. The magnitude of variability is larger in summer (equal to about 7000 grid cells) and smaller in other seasons (minimum equal to about 3000 grid cells in November), but the existence of this year to year variability does not affect the strength of the mean annual cycle of global DA presented in Figure 5.
Remote Sens. 2020, 17, x. https://doi.org/10.3390/rsxxxxx www.mdpi.com/journal/rs At the same time, the presence of DA (yellow curves) over this area and period of year is weak. It should be noted that this misidentification of DA is an inherent limitation not only of the present but of every satellite algorithm using columnar products, in cases of coexistence of sea-salt (coarse strongly scattering) and carbonaceous (fine strongly absorbing) aerosols in the atmospheric column. The only way to overcome this difficulty is to rely not on columnar but on vertically resolved aerosol products, but these are either spatially and temporally limited-for example, those from Cloud-Aerosol Lidar with Orthogonal Polarization, CALIOP onboard the CALIPSO (Cloud Aerosol Lidar and Infrared Pathfinder Satellite Observation) satellite data-or they are model driven-for example, the MERRA-2 data, which are based on the coupled GEOS-5 (Goddard Earth Observing System, [102,103]) and GOCART (Goddard Chemistry Aerosol Radiation and Transport model, [104,105]) models.

Intra-Annual and Inter-Annual Cycles
The intra-annual variability of the frequency of occurrence of dust aerosols on a global basis is depicted in Figure 5, where the total number of grid cells (1° × 1°) over which DA have been observed, in terms of averaged values over the study period 2005-2019, is given. The maximum frequency of DA is noted in June, when dust is observed over almost 23,400 pixels, i.e., over the 3.6% of the planetary surface. Note that 23,400 is the average number of grid cells over which the algorithm detected DA in June (mean value of June's cumulative DA daily frequencies over the period 2005-2019). This number is divided by the average (over the period 2005-2019) number of grid cells of the globe for which the algorithm operated during June, i.e., 633,463 grids, resulting in the low-frequency value of 3.6%. The seasonal cycle of global dust is clear, with minimum frequencies in boreal winter, taking its minimum in November, when dust is observed over only 3181 pixels or less than 1% of the planet. The error bars in Figure 5 correspond to the yearto-year (during the 15-year study period) variability of the frequency of DA on a global basis. The magnitude of variability is larger in summer (equal to about 7000 grid cells) and smaller in other seasons (minimum equal to about 3000 grid cells in November), but the existence of this year to year variability does not affect the strength of the mean annual cycle of global DA presented in Figure 5.  The inter-annual variability of the frequency of occurrence of dust on a global basis during the period 2005-2019 is shown in Figure 6, which displays the annual relative frequency of DA, expressed in terms of percentage of the area of the globe for which the algorithm operated, which (according to the algorithm) is covered by DA. A linear regression was applied to the time series of annual values in order to investigate possible changes and trends. It was found that there is a significant year-by-year variability, within the range of 0.0143-0.0290. The overall maximum frequency of dust on a global scale is observed in 2018, and the minimum in 2007 and 2009. The 2018 maximum can be associated with a minimum in global precipitation (not shown here) based on an analysis using Global Precipitation Climatology Project [106] data. In general, the presence of dust on a global basis increased from 2005 to 2019 by about 56%. This increase is statistically significant at the 95% level, according to the applied Mann-Kendal test results, but also as ascertained from the computed standard error of the slope of the applied linear regression compared to the slope value ( Figure 6).
Remote Sens. 2020, 17, x. https://doi.org/10.3390/rsxxxxx www.mdpi.com/journal/rs algorithm operated, which (according to the algorithm) is covered by DA. A linear regression was applied to the time series of annual values in order to investigate possible changes and trends. It was found that there is a significant year-by-year variability, within the range of 0.0143-0.0290. The overall maximum frequency of dust on a global scale is observed in 2018, and the minimum in 2007 and 2009. The 2018 maximum can be associated with a minimum in global precipitation (not shown here) based on an analysis using Global Precipitation Climatology Project [106] data. In general, the presence of dust on a global basis increased from 2005 to 2019 by about 56%. This increase is statistically significant at the 95% level, according to the applied Mann-Kendal test results, but also as ascertained from the computed standard error of the slope of the applied linear regression compared to the slope value ( Figure 6).

Regional Results
In this section, the seasonal cycle and the inter-annual variation and trends of dust aerosols over ten (10) selected world regions are discussed. For each month (intra-annual variation) or year (inter-annual variation) results are given in terms of relative frequency (sum of grid cells over which dust aerosols are observed by the satellite algorithm divided by the total number of grids of the area for which the algorithm operated, see Section 2.2). The results (statistical metrics computed for each studied region) are summarized and compared in Table 1.

Regional Results
In this section, the seasonal cycle and the inter-annual variation and trends of dust aerosols over ten (10) selected world regions are discussed. For each month (intra-annual variation) or year (inter-annual variation) results are given in terms of relative frequency (sum of grid cells over which dust aerosols are observed by the satellite algorithm divided by the total number of grids of the area for which the algorithm operated, see Section 2.2). The results (statistical metrics computed for each studied region) are summarized and compared in Table 1.

East Sahara
The eastern part of Sahara Desert (ESA, Figure 2) includes a part of Nigeria, Libya, Egypt, Sudan and Chad. In this region, the Qattarah and Bodélé depressions are also included. The latter is characterized as the most intense dust source in the world [4], with high frequencies of DA all year round, as is shown in Figures 3 and 4 (up to 170 days/year and 22 days/month). The regional results of Figure 7a reveal a clear seasonal cycle with maximum frequencies of DA in spring. The dust aerosols are observed most frequently over this region in May, covering an area of about 45,950,000 km 2 , corresponding to 4595 grid cells, slightly more frequently than in April and June (4543 and 4373 grid cells, respectively). During these three months, almost 32% of the included area ( Figure 7a) is covered by DA. The weakest presence of dust over East Sahara is observed in winter, and more specifically in November, when DA are observed over only less than 500 grid cells or 2.3% of the region. As indicated by the standard deviations, which are also displayed in Figure 7a and express the year to year variability, as well as the results of Figure S7ii, this seasonal cycle of DA over East Sahara is followed in every year of the study period. This seasonal cycle of DA over East Sahara can be interpreted in relation with the corresponding cycles of temperature and surface pressure. Indeed, over this region, temperature and surface pressure take lower/higher values in winter and higher/lower values in summer. These conditions favor the uplifting of dust during summer, and also the accumulation of dust aerosols during winter months (because of low wet deposition) leading to a strong presence of DA in spring, which in combination with the gradual heating of surface towards late spring and early summer, result in a late spring/early summer maximum of DA, with a drop in DA frequencies afterward. The inter-annual variability of dust aerosols over East Sahara is presented in Figure  7b. The annual mean values of the relative frequency of occurrence of DA vary between 0.114 and 0.236 pixels or between 11.4% and 23.6% of the region's surface area. According to the applied linear regression, a strong increase (56%) in dust presence over East Sahara during the study period (2005-2019) is observed. This increase was found to be statistically significant according to the applied Mann-Kendall test, as well as based on the computed slope error value, while it is also in agreement with published results in literature [107]   The inter-annual variability of dust aerosols over East Sahara is presented in Figure 7b. The annual mean values of the relative frequency of occurrence of DA vary between 0.114 and 0.236 pixels or between 11.4% and 23.6% of the region's surface area. According to the applied linear regression, a strong increase (56%) in dust presence over East Sahara during the study period (2005-2019) is observed. This increase was found to be statistically significant according to the applied Mann-Kendall test, as well as based on the computed slope error value, while it is also in agreement with published results in literature [107] covering the period 2007-2017. Despite the observed overall increasing tendency of DA during the 15-year period, two different sub-periods with opposite trends can be discerned.
More specifically, from 2005 to 2009, the frequency of occurrence of dust decreased, while it subsequently increased from 2009 to 2018. Such changes in the frequency of occurrence of dust are related with the dust transport [108] and the activation of dust sources [109].

West Sahara
The seasonal cycle and the magnitude of the frequency of occurrence of dust aerosols over West Sahara (WSA, Figure 2) displayed in Figure 8a is different to those for East Sahara. More specifically, the frequency of occurrence of DA is significantly higher, with values up to 6237 grid cells or 51% of the total region's surface area, compared to 4543 pixels and 32% for East Sahara. Besides, over West Sahara, high frequencies of DA occur during summer, with a maximum in July, against a dust cycle with spring maximum over East Sahara. The lowest DA frequencies over West Sahara (varying between 0 and 500 grid cells) are observed in November, December and January. As indicated by the magnitude of error bars in Figure 8a and also by the results of Figure S7iii displaying the annual cycle of DA for every year, the month with the lowest frequency of DA varies, depending on the year, within November, December and January. The seasonality of dust aerosols over West Sahara is driven by Harmattan winds, SHL and the west Sahara monsoons, which blow in the area from the Gulf of Guinea, transporting cold and moist air masses. The prevailing atmospheric conditions (monsoons and SHL) favor convection and as a consequence the uplift of dust and its transport to the upper atmosphere. Once dust is introduced into the atmosphere, it is transported to other regions, namely the Tropical Atlantic Ocean, driven by the atmospheric circulation and especially the Harmattan winds.

East Tropical Atlantic Ocean
As described in the previous sections, dust aerosols can be transported far away from their source areas. This is the case of the Tropical Atlantic Ocean (ETA, Figure 2), where high loads of dust arrive from the Sahara desert, driven by the Harmattan winds, and more specifically by the interaction between the Harmattan winds and the moist west African summer monsoon, producing west African thunderstorms [8,77,78,85,86]. This phe-

East Tropical Atlantic Ocean
As described in the previous sections, dust aerosols can be transported far away from their source areas. This is the case of the Tropical Atlantic Ocean (ETA, Figure 2), where high loads of dust arrive from the Sahara desert, driven by the Harmattan winds, and more specifically by the interaction between the Harmattan winds and the moist west African summer monsoon, producing west African thunderstorms [8,77,78,85,86]. This phenomenon takes place every (boreal) summer, typically between April and November, with different intensity and spatial extent [58,77], leading to dust transport as far as the Caribbean basin [110].
The seasonal variability of the frequency of occurrence of dust over the East Tropical Atlantic Ocean (Figure 9a) reveals a quite clear cycle with a primary maximum in June. This seasonal cycle, also reported in the literature [8,58,78,86], is very similar to the corresponding seasonal cycle for West Sahara (Figure 8a). This is reasonable since dust over the East Tropical Atlantic Ocean originates mainly from central and western Sahara [8,79,86]. However, the comparison between Figures 8a and 9a shows that the frequency of DA is lower over the East Tropical Atlantic Ocean (with values ranging from 153 to 2998 grid cells) than the West Sahara (with values between 550 and 6237 grid cells). It should be noted that standard deviations in Figure 9a (East Tropical Atlantic Ocean) indicate a significant year-to-year variability of DA (also depicted in Figure S7iv), which is quite larger than that in Figure 8a (West Sahara). This, given that West Sahara is the source of dust over the East Tropical Atlantic Ocean, suggests that the larger year-by-year variability of dust over the East Tropical Atlantic Ocean is not only due to a similarly stronger inter-annual variability of dust sources and emitting mechanisms but is also associated with phenomena that drive or influence dust transport [77]. In Figure 9a, apart from the summer maximum in June, a secondary, but much smaller, maximum is also observed in spring (March-April). This maximum can be attributed to dust transport from the eastern Sahara, where dust emissions are maximum in these months (Figure 8a). Indeed, it has been reported in the literature [111] that a great part of dust aerosols in the Tropical Atlantic originates from the Bodélé depression, which is located within the East Sahara study region (Figure 2 of dust over the East Tropical Atlantic Ocean, suggests that the larger year-by-year variability of dust over the East Tropical Atlantic Ocean is not only due to a similarly stronger inter-annual variability of dust sources and emitting mechanisms but is also associated with phenomena that drive or influence dust transport [77]. In Figure 9a, apart from the summer maximum in June, a secondary, but much smaller, maximum is also observed in spring (March-April). This maximum can be attributed to dust transport from the eastern Sahara, where dust emissions are maximum in these months (Figure 8a). Indeed, it has been reported in the literature [111] that a great part of dust aerosols in the Tropical Atlantic originates from the Bodélé depression, which is located within the East Sahara study region (Figure 2).  Figure 9b) and smaller than over the Sahara (West and East). The results also indicate that the occurrence of DA increased during the study period by 30%. This trend is not statistically significant at the 95% level according to the applied Mann-Kendall test, but it is statistically significant based on the computed slope error value). The inter-annual variability (similarly to the seasonal one) of DA over the East Tropical Atlantic Ocean is in quite good agreement with the corresponding variability over Sahara (East and West), revealing a decrease from 2005 to 2011, followed by a subsequent increase until 2018. Figures 7b, 8b and 9b show that there  Figure 9b) and smaller than over the Sahara (West and East). The results also indicate that the occurrence of DA increased during the study period by 30%. This trend is not statistically significant at the 95% level according to the applied Mann-Kendall test, but it is statistically significant based on the computed slope error value). The inter-annual variability (similarly to the seasonal one) of DA over the East Tropical Atlantic Ocean is in quite good agreement with the corresponding variability over Sahara (East and West), revealing a decrease from 2005 to 2011, followed by a subsequent increase until 2018. Figures 7b, 8b and 9b show that there has been an increase of DA during the period of 2005-2019, which gets smaller moving westwards, towards the North African coasts and the adjacent East Tropical Atlantic Ocean. This may be attributed to dust removal mechanisms (wet and dry deposition) as well as to changes in atmospheric phenomena, such as the North Atlantic Oscillation (NAO). Indeed, Ref. [77] found a strong correlation between winter Dust Optical Thickness (DOT) over Northern Tropical Atlantic Ocean and NAO Index (NAOI) (correlation coefficients r = 0.67 and r = 0.74, for TOMS and Meteosat DOT, respectively, from 1984 to 1997). This correlation is even stronger than that reported by [111] for summer using Meteosat data for the period 1984-1994. A similar correlation between dust frequency (DA) and NAOI was found in the present study. More specifically, r values equal to 0.59 and 0.32 were found for winter and summer, respectively. It is also noted that results from a correlation analysis between DOD and NAOI (not shown here) yielded r values equal to 0.47 and 0.35 for winter and summer, respectively.

West Tropical Atlantic
The intra-annual variability of dust aerosols over the West Tropical Atlantic Ocean (WTA, Figure 2) shown in Figure 10a exhibits the same annual cycle with that for the East Tropical Atlantic Ocean, as well as for the West and East Sahara areas, which are the main dust source areas for the Tropical Atlantic. However, as expected, because of the westward transport of African dust, the frequencies of occurrence of dust in this region are further decreased, with respect to those in the East Tropical Atlantic Ocean, varying between 19 grid cells in February and 1095 pixels in July or, in terms of regional surface area coverage by dust, from 0.07% to 11%. The weakened presence of DA over the western than the eastern part of the Tropical Atlantic Ocean is attributed to the gradual deposition of dust (e.g., dry deposition including gravitational settling or wet deposition associated with precipitation and cloud scavenging [112]). The weak presence of dust over the West Tropical Atlantic Ocean is not strange, given the long distance from the Sahara (actually, the western limit of the region is the Caribbean Sea, located 7000 km from the dust source [110]). A noticeable characteristic in Figure 10a is the strong year-by-year variability of the DA's frequency (as indicated by the large error bars in Figure 10a and the diversity of annual curves in Figure S7v). This shows that the dust transport, done by the trade winds and maximized in summer, as well as the dust deposition (wet and dry), vary significantly from year to year, resulting in a strong inter-annual variability. This is also shown in Figure 10b, which displays the series of annual mean frequencies of DA over the study region from 2005 to 2019. On an annual basis, at the West Tropical Atlantic, the frequencies vary between 0.0122 (in 2011) and 0.0518 (in 2018). The applied linear regression and the Mann-Kendall test, but also the computed slope error value, reveal a statistically significant increase in DA during the 15-year study period, equal to 91.8%. This increase, as expected, is in line with the already discussed increase in DA over the East Tropical Atlantic Ocean and the West and East Sahara regions. Of course, the inter-annual tendency of DA over the West and East Tropical Atlantic Ocean does not only depend on the West and East Sahara source regions, but also on the DA transport mechanisms. The identified increase in DA over the West Tropical Atlantic is in line with as well as with increasing NAOI during the study period. Ref. [113] connected the AOD anomaly over the Tropical Atlantic Ocean not only with NAO, but also with ENSO. More specifically, they distinguished between summer (June-July-August) and winter (December-January-February) anomalies, and they associated the summer anomalies with ENSO (r values up to about 0.6-0.7 on a localized grid cell level basis) and the winter anomalies with NAO (r values up to about 0.4, again on a grid cell level basis). On the other hand, Ref. [114] correlated modeled anomalies of summer dust concentration data over the 1965-2008 period with Multivariate ENSO Index (ENSOI) in Barbados and reported an r value equal to 0.20. Here, a correlation coefficient equal to 0.17 was found between summer DA anomalies and ENSOI, whereas a correlation coefficient of 0.12 was found between winter DA anomalies and NAOI. The lower r values in our analysis than those reported by [113] are explained by the fact that our correlation is done using regional DA values averaged over a region in which the r values reported by [113] range from slightly negative up to 0.6 values. It should also be noted that in the present study, the correlation with NAOI and ENSOI is done using DA frequencies, whereas in the other two studies, AOD and dust concentration were correlated with NAOI and ENSOI.

Gulf of Guinea
Apart from the Tropical Atlantic Ocean, desert dust from Sahara is also transported to the Gulf of Guinea (GOG, Figure 2) [55,88,115]. More specifically, during (boreal) winter, when the ITCZ is shifted to the equator, trade winds take a strengthened northeastern direction and transport dust aerosols towards the Gulf of Guinea. Indeed, during winter, 60% of the Sahara desert's DA reach the Gulf of Guinea [80]. On the other hand, during boreal summer, the ITCZ shifts northwards, resulting in eastern trade winds over Sahara, leading to dust export to the Tropical Atlantic Ocean and not to the Gulf of Guinea. This seasonal change in the position and direction of trade winds results in a seasonality of dust over the Gulf of Guinea, which is nicely depicted in Figure 11a. There is a clear annual cycle of DA, with maximum frequencies in winter (December to March, up to 861 grid cells) and minimum in late summer and autumn (August to October, down to 53 grid cells). This seasonal variability of dust is associated not only with trade winds but also with precipitation. During winter, the precipitation amount over the Gulf of Guinea is low (less than 0.5 mm/day, [116]), and the soil of the sub-Saharan Africa (Sahel and sub-Sahel) is arid, so the Harmattan winds transport big quantities of dust to the study region [55]. On the other hand, during summer (when southern winds and monsoons prevail in the region) the amount of precipitation is higher (5-19 mm/day, [116]), and thus the regional dust load is decreased. The year-to-year variability of dust over the Gulf of Guinea is significant, especially in winter, as indicated by the size of error bars as well as by the great variability of winter frequencies of DA over the region from one year to another (see Figure S7vi). Besides, the year-to-year variability of annual frequencies of DA (Figure 11b) is considerable (values within the range 0.0520-0.1483), while there has been an increasing trend (39.2%) during the 15-year study period, which is however not statistically significant according to the applied Mann-Kendal test, but it is so according to the slope error value.

Gulf of Guinea
Apart from the Tropical Atlantic Ocean, desert dust from Sahara is also transported to the Gulf of Guinea (GOG, Figure 2) [55,88,115]. More specifically, during (boreal) winter, when the ITCZ is shifted to the equator, trade winds take a strengthened northeastern direction and transport dust aerosols towards the Gulf of Guinea. Indeed, during winter, 60% of the Sahara desert's DA reach the Gulf of Guinea [80]. On the other hand, during boreal summer, the ITCZ shifts northwards, resulting in eastern trade winds over Sahara, leading to dust export to the Tropical Atlantic Ocean and not to the Gulf of Guinea. This seasonal change in the position and direction of trade winds results in a seasonality of dust over the Gulf of Guinea, which is nicely depicted in Figure 11a. There is a clear annual cycle of DA, with maximum frequencies in winter (December to March, up to 861 grid cells) and minimum in late summer and autumn (August to October, down to 53 grid cells). This seasonal variability of dust is associated not only with trade winds but also with precipitation. During winter, the precipitation amount over the Gulf of Guinea is low (less than 0.5 mm/day, [116]), and the soil of the sub-Saharan Africa (Sahel and sub-Sahel) is arid, so the Harmattan winds transport big quantities of dust to the study region [55]. On the other hand, during summer (when southern winds and monsoons prevail in the region) the amount of precipitation is higher (5-19 mm/day, [116]), and thus the regional dust load is decreased. The year-to-year variability of dust over the Gulf of Guinea is significant, especially in winter, as indicated by the size of error bars as well as by the great variability of winter frequencies of DA over the region from one year to another (see Figure S7vi). Besides, the year-to-year variability of annual frequencies of DA (Figure 11b) is considerable (values within the range 0.0520-0.1483), while there has been an increasing trend (39.2%) during the 15-year study period, which is however not statistically significant according to the applied Mann-Kendal test, but it is so according to the slope error value.

Mediterranean Basin
The Mediterranean Basin (MB) (MED, Figure 2) located near the greatest world deserts, is one of the most interesting areas with relevance to dust aerosols. DA are transported to the region mainly from Sahara Desert in North Africa, but also from Middle East, under favorable atmospheric conditions [117]. The algorithm results (Figure 12a) indicate that the presence of dust over the MB follows an annual cycle with two maxima. The primary maximum occurs in summer, and more specifically in July, when dust covers 1664 grid cells/month or 7.8% of the surface area of the MB. A secondary maximum occurs in spring (April), corresponding to 1151 grid cells/month or 7% of the regional surface area. According to the existing literature [74] but also the algorithm results, these two maxima arise from different parts of the Mediterranean Basin. The primary summer maximum is attributed to dust transport from the central-western Sahara to the Western Mediterranean, driven by the coexistence of a low-pressure system over the central-western Sahara and a high-pressure system over the eastern Sahara [58]. On the other hand, the secondary spring maximum is related to African dust transport mainly to the central, but also to the eastern Mediterranean, caused by moving low pressures over the Sahara (spring Sharav cyclones, [118]). The time series of mean annual values of DA's frequency over the Mediterranean Basin (Figure 12b) indicate that dust covers from 2.02% up to 7.07% of the study area. The observed inter-annual variability is attributed to changes in the produced dust (from sources that are mainly located in North Africa) as well as changes in dust transport [117] and deposition (e.g., precipitation [119]) mechanisms. For example, the high frequencies of dust observed in 2008 are associated with a high number of active dust sources and strong cyclonic circulation over North Africa in this year [109]. On the other hand, the very high DA frequency in 2018 (highest over the 15-year study period) is related to reduced precipitation and remarkable south-southwestern wind directions over the eastern Mediterranean Basin (cyclonic circulation) and strong winds over N. Africa [120]. The North African dust transport to the Mediterranean is associated with the North Atlantic Oscillation (NAO) and more specifically with the positive phase of NAO [111,121]. Such conditions cause drier winters in the MB [119,122] and, as a consequence, they reduce the dust wet deposition [111,123]. The relationship between NAO and dust over the MB is verified by the algorithm results, which show an increasing trend (statistically significant according to the applied Man Kendall test and according to the computed slope error value) of DA's occurrence (equal to 62%) during the study period,

Mediterranean Basin
The Mediterranean Basin (MB) (MED, Figure 2) located near the greatest world deserts, is one of the most interesting areas with relevance to dust aerosols. DA are transported to the region mainly from Sahara Desert in North Africa, but also from Middle East, under favorable atmospheric conditions [117]. The algorithm results (Figure 12a) indicate that the presence of dust over the MB follows an annual cycle with two maxima. The primary maximum occurs in summer, and more specifically in July, when dust covers 1664 grid cells/month or 7.8% of the surface area of the MB. A secondary maximum occurs in spring (April), corresponding to 1151 grid cells/month or 7% of the regional surface area. According to the existing literature [74] but also the algorithm results, these two maxima arise from different parts of the Mediterranean Basin. The primary summer maximum is attributed to dust transport from the central-western Sahara to the Western Mediterranean, driven by the coexistence of a low-pressure system over the central-western Sahara and a high-pressure system over the eastern Sahara [58]. On the other hand, the secondary spring maximum is related to African dust transport mainly to the central, but also to the eastern Mediterranean, caused by moving low pressures over the Sahara (spring Sharav cyclones, [118]). The time series of mean annual values of DA's frequency over the Mediterranean Basin (Figure 12b) indicate that dust covers from 2.02% up to 7.07% of the study area. The observed inter-annual variability is attributed to changes in the produced dust (from sources that are mainly located in North Africa) as well as changes in dust transport [117] and deposition (e.g., precipitation [119]) mechanisms. For example, the high frequencies of dust observed in 2008 are associated with a high number of active dust sources and strong cyclonic circulation over North Africa in this year [109]. On the other hand, the very high DA frequency in 2018 (highest over the 15-year study period) is related to reduced precipitation and remarkable south-southwestern wind directions over the eastern Mediterranean Basin (cyclonic circulation) and strong winds over N. Africa [120]. The North African dust transport to the Mediterranean is associated with the North Atlantic Oscillation (NAO) and more specifically with the positive phase of NAO [111,121]. Such conditions cause drier winters in the MB [119,122] and, as a consequence, they reduce the dust wet deposition [111,123]. The relationship between NAO and dust over the MB is verified by the algorithm results, which show an increasing trend (statistically significant according to the applied Man Kendall test and according to the computed slope error value) of DA's occurrence (equal to 62%) during the study period, which coincides with a corresponding increasing trend in North Atlantic Oscillation Index (NAOI).

Taklamakan Desert
The Taklamakan desert (TAK, Figure 2) is an important dust source, which provides the atmosphere with approximately 200 Tg of dust per year [23]. According to the existing literature, dust emissions in the area are maximum during spring months, after the dry season (minimum precipitation in winter, [96]) and during the strongest uplifting [124]. On the other hand, dust emissions become minimum during winter (November-December-January) after the rainy monsoon period and when downward vertical winds prevail [23,59,124,125]. These seasonal patterns of dust emission are in line with the average intraannual variability of dust over the Taklamakan desert, estimated by the algorithm (2005-2019). Indeed, the results displayed in Figure 13a show a distinct maximum in spring (April) equal to 791 grid cells/month corresponding to a 56% coverage of the Taklamakan study area by DA. Apart from the primary spring maximum of DA, a secondary one in summer is also observed (447 grid cells/month and 34% coverage) which may be attributed to the maximum upward vertical wind velocities prevailing during this season [124]. The inter-annual variability of the frequency of occurrence of DA over Taklamakan (Figures 13b and S7vii)

Taklamakan Desert
The Taklamakan desert (TAK, Figure 2) is an important dust source, which provides the atmosphere with approximately 200 Tg of dust per year [23]. According to the existing literature, dust emissions in the area are maximum during spring months, after the dry season (minimum precipitation in winter, [96]) and during the strongest uplifting [124]. On the other hand, dust emissions become minimum during winter (November-December-January) after the rainy monsoon period and when downward vertical winds prevail [23,59,124,125]. These seasonal patterns of dust emission are in line with the average intra-annual variability of dust over the Taklamakan desert, estimated by the algorithm (2005-2019). Indeed, the results displayed in Figure 13a show a distinct maximum in spring (April) equal to 791 grid cells/month corresponding to a 56% coverage of the Taklamakan study area by DA. Apart from the primary spring maximum of DA, a secondary one in summer is also observed (447 grid cells/month and 34% coverage) which may be attributed to the maximum upward vertical wind velocities prevailing during this season [124]. The inter-annual variability of the frequency of occurrence of DA over Taklamakan (Figure 13b and Figure S7vii

Gobi Desert
The Gobi Desert (GOB, Figure 2) is located in northern China and southern Mongolia, near the previously discussed Taklamakan desert (Figure 2). Although these two great Asian deserts have a quite similar annual cycle of DA, there are slight differences between them. The Gobi has just a single spring maximum, in April (470 grid cells/month and 15% areal coverage, Figure 14a) and not a secondary one in summer, as the Taklamakan does. The same seasonality of DA occurrence over the Gobi was also found by [23], who investigated the spatial and temporal characteristics of dust in East Asia, using vertically resolved CALIPSO data for the period of January 2007-December 2011. The lower dust frequencies over the Gobi than the Taklamakan can be attributed to the higher precipitation in the Gobi compared to the Taklamakan [116]. Concerning the inter-annual variability of dust over the Gobi (Figure 14b), it was found, again, that it is very similar to that of the nearby Taklamakan desert. For example, peaks in 2006, 2010, 2014 and 2018 are observed for both deserts, exhibiting a notable 4-year periodicity. The range of variability of the annual DA frequency is smaller in absolute terms over the Gobi than the Taklamakan (0.0327-0.0729 against 0.2430-0.3809, respectively), but it is larger in relative percent terms (79% against 43%). Such a stronger variability of dust emission and transport over the Gobi than the Taklamakan has also been reported by [59], who used MISR satellite observations and trajectory analysis for the period of 2001-2011, and it is attributed to the complicated meteorological conditions (and especially the frontal activity) over the Gobi [126]. The applied linear regression indicates an increasing (21.1%) but not significant (according to the Mann-Kendall test and the computed error slope value) trend of DA during the period 2005-2019, which is smaller than the corresponding estimated increase for the Taklamakan desert. It should be noted that the estimated increase of the frequency of occurrence of dust over the Gobi is in line with the results published in the literature [127][128][129].

Gobi Desert
The Gobi Desert (GOB, Figure 2) is located in northern China and southern Mongolia, near the previously discussed Taklamakan desert (Figure 2). Although these two great Asian deserts have a quite similar annual cycle of DA, there are slight differences between them. The Gobi has just a single spring maximum, in April (470 grid cells/month and 15% areal coverage, Figure 14a) and not a secondary one in summer, as the Taklamakan does. The same seasonality of DA occurrence over the Gobi was also found by [23], who investigated the spatial and temporal characteristics of dust in East Asia, using vertically resolved CALIPSO data for the period of January 2007-December 2011. The lower dust frequencies over the Gobi than the Taklamakan can be attributed to the higher precipitation in the Gobi compared to the Taklamakan [116]. Concerning the inter-annual variability of dust over the Gobi (Figure 14b), it was found, again, that it is very similar to that of the nearby Taklamakan desert. For example, peaks in 2006, 2010, 2014 and 2018 are observed for both deserts, exhibiting a notable 4-year periodicity. The range of variability of the annual DA frequency is smaller in absolute terms over the Gobi than the Taklamakan (0.0327-0.0729 against 0.2430-0.3809, respectively), but it is larger in relative percent terms (79% against 43%). Such a stronger variability of dust emission and transport over the Gobi than the Taklamakan has also been reported by [59], who used MISR satellite observations and trajectory analysis for the period of 2001-2011, and it is attributed to the complicated meteorological conditions (and especially the frontal activity) over the Gobi [126]. The applied linear regression indicates an increasing (21.1%) but not significant (according to the Mann-Kendall test and the computed error slope value) trend of DA during the period 2005-2019, which is smaller than the corresponding estimated increase for the Taklamakan desert. It should be noted that the estimated increase of the frequency of occurrence of dust over the Gobi is in line with the results published in the literature [127][128][129].

Thar Desert
The Thar Desert is located in the northwestern part of the Indian subcontinent ( Figure  2) and is one of the major Asian dust sources, characterized by extremely dry conditions (rainfall < 0.5 mm/year, Adler et al. 2003) and the absence of vegetation. The intra-annual variability of DA frequency of occurrence in this region is influenced by the South Asian monsoon (June to September, [130]). As observed in Figure 15a, the frequency of DA is maximum in late spring and early summer (dry period), when dust dominates over more than 30% of the region's surface area. On the other hand, dust emissions become minimum over the area after the wet period (October to January), when, according to the algorithm, the DA frequency is almost zero. Regarding the inter-annual variability of the DA frequency, the results of the algorithm indicate a slight increase during the study period (2005-2019). However, this increase is not statistically significant neither according to the applied Mann-Kendal test nor according to the slope error value. A decrease in dust activity in the region, connected with an increase in precipitation, has been reported in the literature [131,132], but it should be noted that these previous studies refer to DOD and not to the frequency of dust occurrence, which is the case here. In addition, and this might be more important, these studies do not include 2018, which is a dusty year that strongly influences the estimated trend. Indeed, when the last two years, 2018 and 2019, are excluded, the trend turns to be slightly decreasing (−9.5%).

Thar Desert
The Thar Desert is located in the northwestern part of the Indian subcontinent ( Figure 2) and is one of the major Asian dust sources, characterized by extremely dry conditions (rainfall < 0.5 mm/year, Adler et al. 2003) and the absence of vegetation. The intra-annual variability of DA frequency of occurrence in this region is influenced by the South Asian monsoon (June to September, [130]). As observed in Figure 15a, the frequency of DA is maximum in late spring and early summer (dry period), when dust dominates over more than 30% of the region's surface area. On the other hand, dust emissions become minimum over the area after the wet period (October to January), when, according to the algorithm, the DA frequency is almost zero. Regarding the inter-annual variability of the DA frequency, the results of the algorithm indicate a slight increase during the study period (2005-2019). However, this increase is not statistically significant neither according to the applied Mann-Kendal test nor according to the slope error value. A decrease in dust activity in the region, connected with an increase in precipitation, has been reported in the literature [131,132], but it should be noted that these previous studies refer to DOD and not to the frequency of dust occurrence, which is the case here. In addition, and this might be more important, these studies do not include 2018, which is a dusty year that strongly influences the estimated trend. Indeed, when the last two years, 2018 and 2019, are excluded, the trend turns to be slightly decreasing (−9.5%).

Thar Desert
The Thar Desert is located in the northwestern part of the Indian subcontinent ( Figure  2) and is one of the major Asian dust sources, characterized by extremely dry conditions (rainfall < 0.5 mm/year, Adler et al. 2003) and the absence of vegetation. The intra-annual variability of DA frequency of occurrence in this region is influenced by the South Asian monsoon (June to September, [130]). As observed in Figure 15a, the frequency of DA is maximum in late spring and early summer (dry period), when dust dominates over more than 30% of the region's surface area. On the other hand, dust emissions become minimum over the area after the wet period (October to January), when, according to the algorithm, the DA frequency is almost zero. Regarding the inter-annual variability of the DA frequency, the results of the algorithm indicate a slight increase during the study period (2005-2019). However, this increase is not statistically significant neither according to the applied Mann-Kendal test nor according to the slope error value. A decrease in dust activity in the region, connected with an increase in precipitation, has been reported in the literature [131,132], but it should be noted that these previous studies refer to DOD and not to the frequency of dust occurrence, which is the case here. In addition, and this might be more important, these studies do not include 2018, which is a dusty year that strongly influences the estimated trend. Indeed, when the last two years, 2018 and 2019, are excluded, the trend turns to be slightly decreasing (−9.5%).

North Middle East
According to the algorithm results, dust aerosols are most frequently observed over the North Middle East desert (NME, Figure 2) from March to July, with a peak in May (Figure 16a, 729 grid cells/month), when they cover 21% of the total surface area of the region. On the contrary, during winter, their frequency is reduced and becomes minimum in November (92 grid cells/month or 2% of the regional surface area). Nevertheless, the large year-to-year variability (large error bars in Figure 15a and diversity of annual curves of DA frequency in Figure S7x), indicating a great year-by-year variability, partly reduces the certitude of the aforementioned annual cycle of DA. This seasonal variability of dust over North Middle East agrees with previous studies, such [45], based on groundbased and satellite observations as well as model simulations, covering the period 2004-2008. The algorithm results show that dust is observed throughout the year. Indeed, it is known that low altitude, soil erosion (due to cultivation and deforestation) and the Tigris and Euphrates rivers, which bring sediments in the region, are some of the factors making North Middle East an important dust source [4,44]. In addition, dust over North Middle East originates not only from regional sources (e.g., Iraq) but also from northern Africa [130] through cyclonic circulation occurring above the northern Middle East and the eastern Mediterranean Basin [133].

North Middle East
According to the algorithm results, dust aerosols are most frequently observed over the North Middle East desert (NME, Figure 2) from March to July, with a peak in May (Figure 16a, 729 grid cells/month), when they cover 21% of the total surface area of the region. On the contrary, during winter, their frequency is reduced and becomes minimum in November (92 grid cells/month or 2% of the regional surface area). Nevertheless, the large year-to-year variability (large error bars in Figure 15a and diversity of annual curves of DA frequency in Figure S7x), indicating a great year-by-year variability, partly reduces the certitude of the aforementioned annual cycle of DA. This seasonal variability of dust over North Middle East agrees with previous studies, such [45], based on ground-based and satellite observations as well as model simulations, covering the period 2004-2008. The algorithm results show that dust is observed throughout the year. Indeed, it is known that low altitude, soil erosion (due to cultivation and deforestation) and the Tigris and Euphrates rivers, which bring sediments in the region, are some of the factors making North Middle East an important dust source [4,44]. In addition, dust over North Middle East originates not only from regional sources (e.g., Iraq) but also from northern Africa [130] through cyclonic circulation occurring above the northern Middle East and the eastern Mediterranean Basin [133]. Figure 16b reveals that there is great year-by-year variability of dust, with annual mean frequencies ranging from 0.0499 to 0.1638. A strong dust activity is observed during the period of 2008-2012, which is suppressed in the subsequent years. The applied linear regression results in a weak decrease (−2.4%), which is not statistically significant according to either the Mann-Kendall test or the slope error value. This variability of dust over North Middle East has been associated with phenomena such as the Southern Oscillation (SO). Thus, Ref. [97]

South Middle East-Arabia
The maximum frequencies of occurrence of dust over the South Middle East desert (SME, Figure 2) are observed a bit later than over the North Middle East desert, namely in June (Figure 17a, 1903 grid cells/month, corresponding to 41% of the total surface area of the region). Such a seasonal variation has been reported by [4] but referring to the period 1980-1992 and was claimed to be connected to the prevailing sea-level pressure and updraft motions as well as to precipitation and soil moisture [134]. Indeed, during winter, the region is under the prevalence of the Siberian high-pressure system, while in summer it is affected by the Pakistan heat low favoring convection. On the other hand, precipitation is minimum in summer and higher in winter (maximum in spring, [116]). The inter-annual variability of the frequency of occurrence of dust (Figure 16b)

South Middle East-Arabia
Τhe maximum frequencies of occurrence of dust over the South Middle East desert (SME, Figure 2) are observed a bit later than over the North Middle East desert, namely in June (Figure 17a, 1903 grid cells/month, corresponding to 41% of the total surface area of the region). Such a seasonal variation has been reported by [4] but referring to the period 1980-1992 and was claimed to be connected to the prevailing sea-level pressure and updraft motions as well as to precipitation and soil moisture [134]. Indeed, during winter, the region is under the prevalence of the Siberian high-pressure system, while in summer it is affected by the Pakistan heat low favoring convection. On the other hand, precipitation is minimum in summer and higher in winter (maximum in spring, [116]. The interannual variability of the frequency of occurrence of dust (Figure 16b) reveals higher annual mean values and stronger year-by-year variability (0.1093-0.2463 versus 0.0497-0.1638) compared to the North Middle East desert. Dust activity in the region has been reported to be connected to ENSO. More specifically, Ref. [97] found that during El Nino (negative phase) the dust optical depth is reduced, while during La Nina it is increased. The applied linear regression to the time series of annual mean frequencies reveals a not statistically significant at the 95% level increasing trend (29.4%) according to the Mann-Kendall test, but marginally significant according to the computed slope and slope error values.

Conclusions
A satellite algorithm was used for detecting dust aerosols (DA) and building a climatology of the frequency of occurrence of DA for the 15-year period 2005-2019. The algorithm uses daily information of MODIS-Aqua based aerosol optical depth (AOD) and Ångström exponent (a) and OMI-Aura based Aerosol Index (AI). Applying specific thresholds to these aerosol properties, the algorithm detects coarse and absorbing dust aerosols over grid cells at 1° × 1° latitude-longitude resolution all over the globe. After its

Conclusions
A satellite algorithm was used for detecting dust aerosols (DA) and building a climatology of the frequency of occurrence of DA for the 15-year period 2005-2019. The algorithm uses daily information of MODIS-Aqua based aerosol optical depth (AOD) and Ångström exponent (a) and OMI-Aura based Aerosol Index (AI). Applying specific thresholds to these aerosol properties, the algorithm detects coarse and absorbing dust aerosols over grid cells at 1 • × 1 • latitude-longitude resolution all over the globe. After its operation on a daily basis over the period of 2005-2019, the algorithm computed the frequency of occurrence of DA on a mean monthly and annual basis. In addition, 15-year average frequencies of DA as well as regional averages for specific world areas of interest (great deserts and downwind regions) were calculated. Finally, the inter-annual variability of DA over the globe and the selected study regions was studied, emphasizing the observed possible trends.

•
The greatest frequencies of occurrence of dust (40-165 days/year) were observed over the North Hemisphere, across the Global Dust Belt (GDB) extending from the western coast of North Africa to China, including the regions of Middle East and Central and Southern Asia. GDB includes the greatest world deserts, namely the Sahara, Middle East, Arabian, Taklamakan and Gobi, from which dust aerosols are transported to remote regions such as the Tropical Atlantic Ocean and the Mediterranean Basin. The maximum frequencies (165 days/year) were observed over the western part of the Sahara and the broader Bodélé area, whereas very high frequencies were observed over the Taklamakan (140 days/year) and Gobi (110 days/year) deserts in Asia. High frequencies of DA (80 days/year) were also observed over the northern part of Middle East as well as over the eastern part of Sahara. Over all these areas, dust aerosols were encountered in more than 50% of the days for which the satellite algorithm operated. • A notable seasonal and inter-annual variability of dust aerosols was revealed by the results of the satellite algorithm. In general, the strongest frequency of occurrence on a global basis was observed in boreal summer, and more specifically in June, when dust aerosols were detected in more than 10 days/month over extended desert areas, namely the Sahara, Middle East and Asian deserts. During June, dust aerosols cover about 3.6% of the global areas for which the algorithm operated (i.e., 633,463 grid cells). On the other hand, the minimum frequencies of DA were observed in boreal winter, and more specifically in November, when dust was observed over desert areas only on a few days (mainly up to 2 or 3 days/month). This seasonal variability of DA on a global basis slightly differs from one region to another. Thus, DA frequencies are maximum in spring (mainly late) over East Sahara, Taklamakan, Gobi, Thar and North Middle East deserts, and in (mainly early) summer over West Sahara, Tropical Atlantic Ocean and the South Middle East desert. Double maximum frequencies, primarily in summer and secondary in spring were are observed over the Mediterranean basin.  56.2%, which is not statistically significant according to the applied Mann-Kendal test but marginally significant according to the slope error. Increasing trends of the frequency of occurrence of dust were observed for all the selected study regions, except for the North Middle East, where decreasing trend was noted. The strongest increase, equal to 92% (statistically significant according to the applied Mann-Kendal test) was observed over the West Tropical Atlantic Ocean. According to the applied Mann-Kendall test, the statistically significant trends are those calculated for the Globe, East Sahara, West Sahara, Mediterranean, Taklamakan and West Tropical Atlantic Ocean. However, when the statistical significance is ascertained from the computed standard error of the slope of the applied linear regression, three (3) other regions were found to have undergone significant trends, namely East Tropical Atlantic Ocean, South Middle East and Gulf of Guinea. Such inter-annual trends of the frequency of occurrence of global or regional dust may have significant implications in various aspects ranging from climate to health and ocean biogeochemistry.

•
The seasonal and inter-annual variability of the presence of dust is connected to dust sources' activation (e.g., desertification) and/or to dust removal and deposition mechanisms (e.g., changes in precipitation). Years with low precipitation are associated with strong dust presence in the atmosphere due to the reduced wet deposition [124]. Thus, for example, the strong increase in DA over the West Tropical Atlantic Ocean, compared to the slighter increase over the East Tropical Atlantic Ocean, given that both regions undergo transport of dust from Sahara, can be attributed to (weakened) deposition.

•
The present algorithm has limitations, since it uses columnar aerosol products. Thus, in cases where coarse strongly scattering sea-salt aerosols in the oceanic boundary layer coexist with overlying (in the free troposphere) fine and strongly absorbing carbonaceous aerosols, the algorithm misidentifies the predominant aerosol type in the atmospheric column as dust (coarse moderately absorbing). This is the case in southeast Atlantic Ocean, off the coasts of Angola, where biomass burning aerosols from the African mainland are exported westwards, into the ocean, from July to September and overlie the marine boundary layer sea-salt aerosols.
In the future, the results of the present study will be compared to similar results obtained with the same algorithm using different satellite data, such as MISR (Multi-angle Imaging SpectroRadiometer). In addition, the current algorithm, using columnar satellite data, can be improved by using contemporary, but spatially and temporally limited as of now, vertically resolved satellite products (for example CALIOP-CALIPSO data) or climatological 3-D aerosol climatologies (MERRA-2), thus providing more accurate results as well as information about the altitude and layers of dust aerosols.