15-year analysis of direct effects of total and dust aerosols in solar radiation/energy over the Mediterranean Basin

: The direct radiative effects of atmospheric aerosols are essential for climate, as well as for other societal areas, like the energy sector. The goal of the present study is to exploit the newly developed ModIs Dust AeroSol (MIDAS) dataset for quantifying the direct effects on the downwelling surface solar irradiance (DSSI), induced by the total and dust aerosols amount, under clear-sky conditions and the associated impacts on solar energy for the broader Mediterranean basin, over the period 2003 – 2017. Aerosol optical depth (AOD) and dust optical depth (DOD) derived by the MIDAS dataset, along with additional aerosol and dust optical properties and atmospheric variables were used as inputs to radiative transfer modeling to simulate DSSI components. A 15-year climatology of AOD, DOD and of clear-sky Global Horizontal Irradiation (GHI) and Direct Normal Irradiation (DNI) was derived. The spatial and temporal variability of the aerosol and dust effects on the different DSSI components was assessed. Aerosol attenuation of annual GHI and DNI range from 1-13% and 5-47%, respectively. Significant attenuation by dust 2-10% and 9-37% was found over North Africa and the Middle East, contributing to 45-90% of the total aerosol effects. The mean GHI and DNI attenuation during extreme dust episodes reached values up to 12% and 44%, respectively, for different areas. After 2008 a decline of aerosol effects on DSSI was found, attributed mainly to the dust component. Sensitivity analysis using different AOD/DOD inputs from Copernicus Atmosphere Monitoring Service (CAMS) reanalysis dataset, revealed CAMS underestimation of aerosol and dust radiative effects compared to MIDAS, due to slight AOD and stronger DOD underestimation, respectively.


Introduction
Aerosols modulate the radiation field of the earth-atmosphere system with several implications for life on earth. The physical mechanisms through which they influence the radiation budget are manifold. Aerosols interact directly (direct effects) with the shortwave (SW) and longwave (LW) radiation through scattering and absorption. Moreover, through their interactions with clouds, they have semi-direct and indirect effects, by altering the atmospheric conditions related to cloud formation/dissipation due to absorbing aerosols and by acting as cloud condensation and ice nuclei (altering the microphysical and hence the optical properties of clouds) (e.g. [1]). The direct radiative effects of aerosols refer to changes in radiative fluxes due to the direct aerosol-radiation interaction (REari -Radiative Effects due to aerosol-radiation interactions as was renamed in the IPCC AR5 [2]). Cloud free aerosol direct radiative effects depend on aerosol optical properties the spectrally resolved aerosol extinction coefficient or the integrated Aerosol Optical Depth (AOD), the Single Scattering Albedo (SSA) and scattering phase function or integrated asymmetry parameter and other environmental parameters like surface reflectance and atmospheric trace gases [3]. Their magnitude corresponds to the perturbation of the radiation fields induced by aerosols, within the Earth-Atmosphere system, with respect to an atmospheric state without their presence. Based on their sign, it is expressed the gain (loss) of energy, for the SW and LW radiation, at the surface and at the top of at the atmosphere when positive (negative). Due to the aerosol induced changes on the Earth's energy budget, they are among the main climate change drivers. Aerosol radiative forcing estimates constitute the quantification of the aforementioned effects due to the aerosol perturbations induced by anthropogenic activity, useful in the framework of climate change policy , as it can be related to changes in global mean surface temperature [4]. Their estimates to the total anthropogenic radiative forcing are related with the greatest uncertainties [5]. Several reasons contribute to these not yet well-constrained estimates, such as the heterogeneity of the processes governing aerosols' production and removal, which in turn regulate optical and microphysical properties, both determining the associated aerosol-radiation interactions. Focusing on earth's surface, apart from the importance of aerosols on climate, studies quantifying the impact of total aerosols (natural and anthropogenic) on incoming solar irradiance are also essential for other societal benefit areas such energy sector (e.g. [6,7]).
Dust particles constitute a major component of atmospheric aerosol load [8]. Their emission is primarily wind driven over the arid or semi-arid regions of the planet. The Sahara Desert in North Africa, the Arabian and the Asian deserts are the major dust sources on Earth, with the highest contribution to global dust load (more than 50%) to be emitted from North Africa [9][10][11]. Even though the dust sources are localized, the spatial distribution of dust over the globe is extensive due to dust mobilization under favorable meteorological conditions. The Mediterranean basin located to the crossroads of air masses from all over the globe [12], experiences high aerosol concentrations of both natural and anthropogenic origin and the aerosol spatiotemporal variability over the area has been investigated in several studies e.g. [13][14][15][16]. In addition, due to its proximity to the Earth's most active dust sources located across the Sahara Desert and Middle East deserts, the Mediterranean basin often experiences high dust aerosol loads (dust intrusions). Several studies have investigated the Saharan dust transport towards the Mediterranean Basin by exploiting, either solely or combined, in-situ measurements, remote sensing (ground based or satellite) retrievals and atmospheric-dust numerical products (e.g. [17,18]). Across the Mediterranean Sea, dust concentrations fade down from south to north attributed to the activation of physical mechanisms responsible for the removal of mineral particles from the atmosphere, either due to dry or to wet deposition [19]. It should be noted that long-range transport of Saharan mineral particles has also been recorded over northern Europe [20,21]. Additionally, there is a distinct seasonal cycle of the longitudinal spatial distribution of dust with higher concentrations in the eastern and central parts of the Mediterranean in spring that are shifted to the western parts in summer [19,22,23]. Those seasonal patterns have been related to the prevailing meteorological conditions over the area [19,24]. This intricate aerosol regime constitutes Mediterranean basin on of the most interesting areas for investigating the aerosol and dust radiative effects.
Towards the reduction of greenhouse gas emissions and climate change mitigation, the deployment of renewable energy technologies is one of the main pillars [25], while at the same time they provide energy in a sustainable manner. Among renewables, the share of solar technologies is growing fast due to the combination of the substantial reduction in their cost and the fact that solar energy is the most abundant among other renewable resources. For electricity production, the systems in use are the photovoltaic (PV) cells and -Investigate the advantages and limitations of existing model and satellitebased aerosol time series for solar energy related applications. After the description of the datasets and the methodology in Section 2, the AOD and DOD climatology is discussed in Section 3.1 along with the intercomparison between MIDAS and CAMS datasets. Section 3.2 presents the aerosol and dust effects on the surface GHI and DNI based on both datasets. In Section 3.3, focus is given on the radiative effects for each DSSI component under extreme dust outbreaks. The interannual variability of the effects of aerosols and dust on GHI and DNI components is presented in Section 3.4 and the clear-sky climatology of GHI and DNI is given in Section 3.5. Finally, our concluded remarks are provided in Section 4.

Data and Methodology
In the present study, the aerosol and especially dust effects on DSSI in terms of GHI and DNI were investigated over the Mediterranean basin. Initially, two different AOD and DOD datasets were explored, the newly developed satellite-based MIDAS and the model derived CAMS reanalysis datasets. As a second step, the quantification of the longterm effects of total aerosols and dust on different DSSI components with respect to aerosol free conditions were derived, using both the aforementioned AOD and DOD datasets as inputs to RTM simulations. For the simulations, additional aerosol optical properties and atmospheric parameters were used as inputs as well. An intercomparison of the results from the two datasets was performed in order to address the last two scientific questions listed in Section 1. The description of the utilized datasets and the RTM simulations are provided in Sections 2.1 and 2.2 respectively.

Data
We have used various aerosol and atmospheric related parameters as inputs in the RTM. Table 1 presents an overview of those datasets and a more analytical description is provided in the corresponding subsections.  15-year period (2003-2017). In brief, MIDAS DOD product has been developed through the synergy of quality assured MODIS-Aqua Level 2 AOD and the dust fraction (MDF) to the total aerosol load, in optical terms, acquired from the MERRA-2 reanalysis. Along with DODs, it is provided also the associated grid-cell uncertainty estimated using reference AERONET retrievals [42] and LIVAS [43] products for AOD and MDF, respectively. A comprehensive evaluation of the MIDAS DOD versus AERONET DOD-like and an intercomparison against MERRA-2 and LIVAS DODs justified its reliability and validity as well as its caveats which should be taken into account. For the estimation of the radiative effects attributed to the total aerosol load, we are using the MODIS-Aqua AOD stored in the MIDAS files. Actually, the MIDAS AOD is the raw MODIS-Aqua AOD (Collection 6.1; [44]), by merging Dark Target and Deep Blue retrievals according to Sayer et al. [45], in which quality filters (see Section 2.1 in [31]) have been applied and it has been reprojected on an equal latitude-longitude grid as those of DOD.

AOD, DOD and Total Column Water Vapor (TCWV) -Copernicus Atmospheric Monitoring Service (CAMS) reanalysis dataset
The CAMS reanalysis, available from 2003 onwards, is the ECMWF's (European Centre for Medium-Range Weather Forecasts) global reanalysis dataset of atmospheric composition, consisting of three-dimensional time-consistent atmospheric composition fields, including aerosols and chemical species [40]. It is based on ECMWF's Integrated Forecast System (IFS), including an aerosol module described in Morcrette et al. [46]. Five species of tropospheric aerosols are included in the CAMS aerosol model, including dust. For dust sources, the parameterization of Ginoux et al. [47] is implemented. The satellite derived aerosol products that were assimilated in the CAMS reanalysis were the MODIS-Aqua and MODIS-Terra AOD retrievals [37] and, in addition, the retrievals from the Advanced Along-Track Scanning Radiometer (AATSR) onboard Envisat from 2003 to March 2012. More details regarding the updates in the meteorological part of IFS and in the aerosol and chemical modules, the data assimilation process and the emission datasets are given in Innes et al. [40] and the references therein. CAMS reanalysis products are available from the Copernicus Atmosphere Data Store (ADS, https://ads.atmosphere.copernicus.eu/#!/home) on 3-hourly basis. Using the CDS API service, AOD and DOD at 550nm and TCWV were obtained for the same period with MIDAS dataset (2003-2017) on a 0.4 o x 0.4 o lat/lon grid.

Additional aerosol and dust optical properties (SSA and AE)
For the additional UV-visible-near IR optical properties for all aerosols and dust climatological values from the second version of the Max-Planck Aerosol Climatology (MACv2; [48]) were utilized, which are given at a global scale with a 1° x 1° spatial resolution. The interannual variability of total aerosols optical properties is provided over the time period 2001-2016 as monthly means. For five aerosol species, including dust, monthly climatological values (intra-annual variability) are given, representative of the aforementioned period.
In the current analysis, SSA at 550 nm was used for total aerosols and dust (DU SSA). AODs at 470 nm and 850 nm, for total aerosols, were obtained from MACv2 and the Ångström exponent (AE) was calculated between the two selected wavelengths (AE 470-850nm). For dust, a fixed climatological value of 0.4 for AE 440-675 nm was used as proposed for the region by Taylor et al. [49].

Total Ozone Column (TOC)
To obtain TOC data for the entire period, a new dataset was constructed by combining data from Ozone Monitoring Instrument (OMI) onboard NASA's Aura satellite from 1/10/2004 until 31/12/2017 and from Total Ozone Mapping Spectrometer (TOMS) onboard the Earth Probe (EP) satellite from 1/1/2003 to 30/9/2004. Τhe satellite-based TOC retrievals were collected from the daily global OMI TOMS-Like TOC Level 3 (OMTO3d) gridded on a 1° x 1° grid product [50] and from the EP TOMS Level 3 (TOMSEPL3) version 8 product [51], which provides daily data on a global grid of 1° x 1.25°.

Library for Radiative transfer (libRadtran) simulations
The simulations of DSSI components for cloud-free conditions were performed using the uvspec model from the libRadtran package [52]. Using the radiative transfer solver sdisort [53], pseudospectral simulations were performed every 1nm for the spectrum range of 280-3000nm, using for the molecular absorption the parameterization of LOWTRAN band model [54], as adopted from the SBDART code [55]. The Kurucz 1.0 nm [56] extraterrestrial solar spectrum and the standard US atmospheric profile [57] were utilized and the surface albedo was set to 0.2.

Database (DB) for Radiative properties
For quantifying the impact of total aerosols and dust on the DSSI components (GHI and DNI), we have performed RTM simulations (see section 2.2.3) using as inputs satellite (MIDAS) retrievals and reanalysis (CAMS) products of AOD and DOD, complemented by additional aerosol optical properties and atmospheric parameters acquired from the MACv2 climatology, spaceborne observations (OMI, TOMS) and reanalysis products (CAMS), which are described in detailed in the previous section.
Those datasets differ in spatial and temporal resolution. First, the spatial and temporal homogenization of datasets with missing values was performed and then the geolocation and synchronization among all datasets in order to generate a complete database (DB) of all the input parameters needed for the RTM simulations (SZA, AOD or DOD, SSA, AE, WV, TOC) on a 0.4° x 0.4° lon/lat grid, on an hourly basis, which was selected to be the frequency of RTM simulations, in order to account for the sun elevation. Figure 1 provides a schematic overview of this process.
We selected to aggregate the MIDAS dataset to the coarser spatial resolution of 0.4°x0.4°. This in order to avoid datasets with larger number of missing values and to improve the daily and monthly statistics of AOD and DOD. The median value was selected as the aggregation method, as a non-parametric measure of central tendency based on the findings of Sayer and Knobelspiesse [58]. In order to homogenize the MIDAS dataset in time and space, the missing values were filled by monthly means. Seasonal means were utilized when the monthly data availability was low (<20%). In cases with low seasonal availability the spatial gaps were filled using bilinear interpolation. For the 1h RTM simulations, the daily MIDAS values assumed constant within day. We took into account the diurnal variability of CAMS datasets (AOD, DOD and TCWV) and the 3h values were assumed constant within the 3-hour time interval. For TOC products, a temporal fitting was performed filling the days missing with monthly mean values and a spatial bilinear fit was performed in order to fill the gaps in space. The original spatial resolution of TOC datasets was 1° x 1° and 1° x 1.25° for OMI and TOMS, respectively, but they were both regridded to 0.4° x 0.4°, using bilinear interpolation. The daily values were kept constant for the 1h time interval of the RTM simulations. Using bilinear interpolation, the 1° x 1° fields of MACv2 were re-gridded to the 0.4° x 0.4° equal lat-lon grid. The monthly values of SSA and AE were assumed constant within every month.
Those assumptions and especially for AE and SSA are simplifications, but they are secondary compared to AOD and DOD changes, which variations are also minor compared to daytime differences due to SZA.

RTM methodology
The hourly values of clear-sky DSSI in terms of global and direct components were obtained relying on pre-calculated Look-Up Tables (LUTs), for reducing the computational cost, similarly done in Kosmopoulos et al. [59]. Spectral LUTs were constructed containing simulated surface spectral irradiances (global and direct), for a wide range of SZAs and atmospheric factors causing DSSI attenuation under cloudless conditions, which were used as input parameters in uvspec (Section 2.2.1) and presented in panel (a) of Figure 2. All possible combinations of the input parameters of Figure 2(a) resulted to 74,520 LibRadtran simulations of the spectral LUTs. The output spectral irradiances were integrated over SW to get the total irradiances. In order to discretize further, we applied interpolation on the spectrally integrated values, and finer LUTs (panel (b) of Figure 2) were derived, as they could be resulted from over 200 million hypothetical RTM runs. Using the fine LUTs and inputs from DB described in section 2.2.2, instantaneous values of total irradiances (global and direct) every 1h during daylight were extracted for every grid cell (0.4° x 0.4°) for the study period 2003-2017. This procedure was repeated five times, as different experiments, described in Table 2, in order to quantify the total aerosol as well as the dust effect on different DSSI components, for the two different datasets (MIDAS, CAMS).  Integrating over sunlight hours the 1h instantaneous values of total irradiances, daily Global Horizontal Irradiations (GHI) and Direct Normal Irradiations (DNI) components (in MJ/m 2 ) were calculated and those values were post -corrected for the Earth-Sun distance, and for the surface elevation following the methodology described in Fountoulakis et al. [30].

Satellite derived and modeled AODs and DODs -climatology and intercomparison
One of the aims of this study is to investigate the benefits and the drawbacks of choosing between satellite and model derived AOD and DOD datasets for estimating their radiative effects. To this end, as described in Section 2, MIDAS and CAMS datasets were explored. The geographical distribution of the long-term annual averaged values of AOD (Figures 4a and 4b) and DOD (Figures 5a and 5b) were derived for both datasets. The corresponding season results are given in the Supplement (Figures S1-2, S4-5). It has to pointed out here that our analysis expands further CAMS DOD product evaluation [60] in terms of its spatial and temporal variability performance. The frequency histogram of the CAMS-MODIS biases and their mean annual geographical distribution are presented For this comparison, the aggregated MIDAS datasets were used, before filling the missing values (see section 2.2.2). CAMS datasets which have a diurnal variation (3h time resolution) were synchronized with MIDAS datasets (MODIS Aqua overpass time) in order to achieve an exact collocation. In Figure 3, the MIDAS cloudless sky data availability (expressed in percentage) is illustrated. Three regions can be distinguished. North Africa, which due to its scarce cloudiness, has the highest data availability, with more than 70% of daily satellite retrievals with respect to the whole period. Over the Mediterranean Sea, data availability decreases down to 60%. Over south Europe, the MIDAS data amount further decreases and it is minimized (~20%) in mountainous regions (i.e., Alps). For the temporal aggregation, only grid points with at least 20% data availability on annual and seasonal basis were used, in order to ensure the representativeness of the results.

Aerosol optical depth
In general, the AOD spatial features are similarly reproduced for both datasets (Figure 4 a,b) and the regional averages are almost equal (0.19±0.06 for MODIS and 0.18±0.06 for CAMS, Table 3). Our findings are also in a good agreement with those of previous studies focusing on the same region of interest (ROI) [13,16]. Between the two datasets, differences were found on the maximum AOD levels. The annual mean AOD values for each individual pixel, range from 0.05 to 0.48 for MIDAS while those of CAMS yield values from 0.05 to 0.37. For both datasets, the maximum and minimum AOD values were found over the same areas. Over North Africa and parts of Middle East maximum AOD values were derived mainly attributed to desert dust. Large AOD values related to anthropogenic activity [61,62] were found over the megacity of Cairo, Egypt and the Po Valley, Italy. Low AOD values were found over the greater part of Iberia Peninsula and South France.
Overall, the CAMS simulated AODs are slightly underestimated (mean bias -0.006) with respect to MIDAS AOD (Figure 4c), which was found to be statistically significant on 95% confidence level (p-value<0.05, t-test for the differences). This result agrees with previous CAMS AOD comparison with satellite products [60], where lower CAMS AODs relative to MODIS were reported over NAMEE domain. The geographical distribution of annual mean bias (Figure 4 d) revealed areas with annual mean bias that differ a lot from the regional averaged value. The greatest CAMS AOD underestimation (with annual mean bias up to -0.15) was found for the most southeastern part of Spain, an confined area strongly influenced by aerosols, that was also indicated in previous studies with MODIS data [13]. There are also areas with CAMS AOD overestimation, with the greatest mean annual value (up to 0.1) over Northwest Africa, which may be explained by the CAMS model overestimation of the organic matter over that area [60]. This CAMS AOD overestimation is more pronounced in summer ( Figure S3c). There is a clear seasonal cycle (Figures S1-2) for AOD over the Mediterranean basin. CAMS AOD reproduces MODIS AOD seasonal variability with the same regional patterns, but again there were differences in the magnitude of maximum seasonal averaged AODs between those datasets. In general, seasonal averaged AODs are higher compared to annual means. In summer, AODs are maximized (0.66 for MODIS and not exceeding 0.56 for CAMS) over North Africa, particularly in its western parts. Over south Europe, large AOD values were found over its southeastern part for spring and summer, which are mainly due to emissions of anthropogenic aerosols like sulfates [16], with peak in spring over Po Valley with mean values of 0.37 for MODIS and 0.28 for CAMS. In Table  3, the regional averaged seasonal mean AOD values are summarized. The distinct seasonal cycle revealed in both datasets with maximum values during summer and minimum values during winter, is the same with the seasonal cycle reported by Papadimas et al. [13]. The seasonal variations of AOD have been linked with the atmospheric circulation and the meteorological conditions over the study area that are affecting the aerosol emission, removal, and transport processes [63,64].

Dust optical depth
There is a clear latitudinal gradient of DOD (Figure 5a,b). For both datasets, the larger DOD values were found over North Africa and parts of Middle East, where major dust sources (Sahara and Arabian peninsula deserts) are located [9,11]. But a CAMS DOD deficiency is profound. For the annual averaged MIDAS DOD a maximum of 0.35 was found over a persistently dust active region of salt lakes (local "chotts") and dry lakes in the boarders of Tunisia and northeast Algeria. Large values (0.32) were also found over the desert of central Algeria. Over the eastern Libyan Desert and Egypt, for most individual  (Table 3), there is a slight difference (0.08±0.07 for MIDAS DOD and 0.06±0.06 for CAMS DOD) in absolute values. A relatively high underestimation (mean bias almost -0.03) of CAMS DODs against MIDAS (Figure 5c) was derived, which was found to be statistically significant on 95% confidence level (p-value<0.05, t-test for the differences). This is almost 40% lower values compared to MIDAS, which is in agreement with the underestimation of CAMS DOD (up to 46%) with respect to AERONET observations reported by Bennouna et al. [60] over the same area. According to the same study, the higher CAMS DOD underestimation was found during wintertime, which was attributed to overestimations in biomass burning organic matter (OM) and in secondary organics over heavily populated areas in summer. The geographical distribution of annual mean bias (Figure 5 d) revealed that almost everywhere MIDAS DOD is larger than CAMS DOD. The largest values of CAMS DOD underestimation (up to -0.15) were found over the aforementioned dust sources of the Saharan and Middle East deserts, which is also found for all seasons ( Figure S6). The geographical distribution DOD over the Mediterranean basin follows a seasonal cycle ( Figures S4-5) and according to the regional averaged seasonal DOD (Table 3) it has maximum in spring and summer, in agreement with previous studies over the same area [16,19,22]. Again, the main characteristic comparing the results from the two datasets is the deficiency of seasonal averaged CAMS DOD compared to MIDAS. The dust activity weakens in winter and it is observed mainly in the northeastern Africa, with DODs up to 0.25 and 0.13 for MIDAS and CAMS, respectively. The high dust activity begins in spring, with large DOD values (maximum values up to 0.42 for MIDAS and 0.27 for CAMS) were found over an extended area, the central and eastern parts of North Africa and the part of Middle East. In summer, the high dust activity is limited to northwestern Africa and the highest mean seasonal DOD values were found for this season, with the differences between the two datasets be minimized (0.45 for MIDAS and 0.44 for CAMS). This seasonal cycle of dust activity and transport over the Mediterranean is in agreement with previous studies, which also investigate the atmospheric circulation patterns favoring this cycle [24]. The differences that were found between the two different datasets (MIDAS, CAMS) and especially at the maximum AOD/DOD levels, were investigated further. It was found that a great portion of MIDAS high AOD and DOD values is missing from the CAMS dataset. Table 4 summarizes the amount of data that are higher than 1, 1.5, 2 and 3 in terms of AOD and DOD for both datasets. It also shows the percentage of the missing high values from CAMS datasets compared to MIDAS. For MIDAS AOD dataset the 0.05% are values greater than 2, while this percentage is only 0.0015% for CAMS dataset, resulted to a high segment (97%) of high MIDAS AODs that do not exist in CAMS dataset. For DOD, over 90% of missing high values has a lower threshold of DOD 1.5, due to strong CAMS DOD underestimation. It is clear from the results that for very high aerosol burdens there are significant differences between the explored datasets, especially when considering the dust component, differences that possibly can lead to uncertainties of the aerosol and dust radiative effects.
Finally, based on MIDAS dataset results of this section, the long-term dust contribution to total aerosols in optical terms ranges from 40% to 90%, over North Africa and Middle East, making dust the most important aerosol component over those areas.

Aerosol and Dust effects on DSSI
In this section, the quantification of total aerosol and dust radiative effects on GHI and DNI over the Mediterranean basin is presented. Using the simulated daily irradiations described in section 2.2.3, mean annual and seasonal integrals (INT) of GHI and DNI, for cloudless conditions, were calculated for the five different experiments described in Table 2, and using the following formula: where i stands for aerosols and dust, the relative change (expressed in percentages) on DSSI due to aerosols and dust presence was calculated with respect to an aerosol free atmosphere. Using the same formula, the daily relative changes were calculated as well.
Those changes were normalized with the aerosol free atmosphere and expressed as relative changes, in order to give a common measure of the aerosol/dust effect on DSSI for areas with different latitude and throughout a year. Moreover, in this section we investigated the effects on DSSI when different datasets (MIDAS or CAMS) of AOD/DOD are used. For this intercomparison between MIDAS and CAMS datasets, in total ~23Μ data points were compared. The average data points (days) compared for each or the ~8000 pixels of the Mediterranean basin, were ~3000 per pixel or ~200 per pixel per year. As mentioned, missing data are related with MIDAS gaps due to cloudy pixel scenes.
The change in the mean annual integral of GHI due to the presence of total aerosols and dust is presented in Figures 6 and 7 respectively, estimated using MIDAS (panel a) and CAMS (panel b) datasets. The same set up, but for DNI is given in Figures 8 and 9. The corresponding seasonal results are given in the supplement (Figures S7-14). In all cases, the patterns of GHI and DNI changes are consistent with those of AOD and DOD. The higher the AOD and DOD values are, the higher their radiative effect. Due to the interactions of the incoming solar radiation by the overlying aerosol (dust) layers, the GHI and DNI amounts reaching at the ground are reduced with respect to an aerosol-free atmosphere, thus explaining the existence of negative values throughout the domain. The day-to-day variations of those effects are presented in Figure 10. Hereon, AOD GHI attenuation stands for the GHI reduction by total aerosols and DOD GHI attenuation means reduction of the GHI by dust component. The same nomenclature is used for the DNI component.

Aerosol effects on GHI
There is a qualitative agreement in the geographical distribution of the annual AOD GHI attenuation between the two datasets ( Figure 6) and their regional averages using the MIDAS dataset (5.2%, Table 5) is almost the same with CAMS (5.1%). The differences in the magnitude of annual mean AODs and especially for the maximum values were also inherited to their radiative effects. The long-term GHI reduction due to aerosols were found to range from 1 to 13% for MIDAS and from 2 to 10% for CAMS. The ~23% relative difference in the maximum levels of annual AOD between the datasets (section 3.1) resulted to this difference in the maximum effects on GHI, which in relative terms is of the same magnitude.
In general, three sub-domains (D) are highlighted for the annual AOD GHI attenuation, based on aerosol load spatial patterns. The highest effects were found over North Africa and the Middle East (D1), where the annual AOD GHI attenuation varies from 4% to 13% based on MIDAS dataset (4% to 10% for CAMS). Lower values were found for central and southeastern Europe and Anatolian Peninsula (D2) ranging from 3% to 8% (3% to 7%), with the largest values over the Po Valley. Over the Iberian Peninsula and south France (D3), the lowest annual AOD GHI attenuation was found ranging from 1% to 6% (2% to 5%), with exemption over southeastern Spain, with values up to 8% (MIDAS dataset). There are clear seasonal variations (Figures S7-8) of the AOD GHI attenuation geographical distribution, with seasonal averaged values of GHI reduction reaching higher values than annual ones. The maximum MIDAS AOD GHI attenuation up to 14% was found in spring over North Africa, while the corresponding CAMS values were around 10%. For CAMS, the maximum reduction was 13% and was found for summer over the northwestern Africa, where similar attenuation was estimated from the MODIS dataset. In Table 5, the regional averages of the seasonal AOD GHI attenuation are summarized. The seasonal cycle of AOD GHI attenuation for both datasets differs from the seasonal cycle of the corresponding AOD values. The most notable difference is that the peak of AOD GHI attenuation was derived in winter for MIDAS dataset instead of summer, when the peak of MIDAS AOD was found. The unexpected peak of the GHI attenuation in winter was mainly due to significant attenuation of the GHI over Egypt and Eastern Libya, which was subsequently attributed to minimum seasonal SSA values in winter over the area ( Figure S11a). For CAMS AOD dataset, a small shift of maximum GHI attenuation to spring instead of summer was found.

Dust effects on GHI
Under the absence of non-dust aerosol species, the spatial patterns of GHI attenuation, based on MIDAS and CAMS DODs (Figures 7, S9-10), show a clear south-north gradient regulated by the reduction of dust load amount from sources to distant downwind regions. The CAMS DOD underestimation is also depicted in the corresponding lower attenuation of the GHI for the CAMS DOD. The maximum of annual DOD GHI attenuation values was found over North Africa and parts of the Middle East, ranging from 2 to 10% for MIDAS and from 2 to 8% for CAMS. This attenuation of GHI by dust accounts for ~45%-90% of the overall attenuation by aerosols over this area on annual basis, and this percentage takes higher values (up to 95%) on a seasonal basis over those highly dust affected areas. In summer, the seasonal mean reached values up to 11% for MIDAS and 10.5% for CAMS, over the northwestern Africa. Except for summer, the CAMS DOD GHI attenuation is significantly lower than the MIDAS DOD GHI attenuation for the rest of the year. Regarding the regional averaged values (Table 5), the annual GHI attenuation due to MIDAS DOD (2.4%) is almost 30% larger than the attenuation from CAMS DOD (1.7%). The seasonal cycle of GHI attenuation attributed to dust is the same with the seasonal cycle of DOD with maximum in spring and summer, but for MIDAS dataset the spring peak is higher than summer (21%) which could not be explained only by the corresponding DOD differences between the two seasons (9%). The sharp MIDAS peak in spring may also be explained by the small DU SSA over North Africa and parts of Middle East during this season ( Figure S12b).  Table 5. Change (in %) of the regional averaged mean annual and seasonal integrals of GHI due to total aerosols (AOD) and dust (DOD) from both datasets of MIDAS and CAMS. The maximum of seasonal values is denoted as bold, emphasizing the peak of the seasonal cycle.

Aerosol effects on DNI
In general, the attenuation of DNI by total aerosol and dust is much larger than the attenuation of GHI. Aerosols can remove significant fraction of the solar radiation from the direction of the solar beam (i.e., DNI) mainly through scattering. Most of the scattered radiation interacting with aerosols (including the scattered radiation which has been removed from the direct solar beam) is redistributed in the atmosphere through multiple scattering. Thus, aerosols attenuate less effectively the scattered (and correspondingly the total) relative to the direct solar radiation. This is not the case for highly absorbing aerosol mixtures, which however are not common and exist mainly over polluted urban environments.
The average of AOD DNI attenuation (Figure 8) ranged from 5% to 47% for MIDAS and from 10% to 39% for CAMS. The AOD differences between the CAMS and MIDAS datasets were amplified in terms of DNI attenuation. For D1 the maximum of DNI attenuation was found for both data sets, with values ranging from 15% to 47% for MIDAS and to 39% for CAMS. In this domain, areas like Morocco, North Algeria, North Tunisia and the areas around Red Sea, had annual attenuations less than 20% which make these areas favorable for CSP (DNI related) installations. For D2, an average DNI reduction between 15 and 25% was derived for both datasets, except for Po Valley where values up to 32% and 26% were found for MIDAS and CAMS respectively. The lowest values were found for D3 ranging from 5% to 25% (except for southeastern Spain where it was 35% for MIDAS) and from 10% to 20% for CAMS. The seasonal AOD DNI attenuation values (Figure S11-12) reached higher values up to 53% for MIDAS and 49% for CAMS, which were found over northwest Africa in summer. The resulted seasonal cycle (Table 6) of aerosol effects on DNI was the same with the seasonal cycle of corresponding AOD values, for both datasets, with maximum during summer and minimum during winter.

Dust effects on DNI
The peak of DNI attenuation due to dust was found over North Africa and Middle East, with values ranging from 9 to 37% for MIDAS and from 9 to 28% for CAMS ( Figure  9). In summer, over northwestern Africa those reductions reached values up to 40% and 38% for MIDAS and CAMS ( Figures S13-14), respectively. The contribution of dust to the overall DNI attenuation by aerosols varies between ~45-90%. For the regional averaged values (Table 6) a larger DNI attenuation due to MIDAS DOD (10.7%) was found than due to CAMS DOD (7.5%), a difference that is attributed to the strong underestimation of CAMS DOD.  Table 6. Change (in %) of the regional averaged mean annual and seasonal integrals of DNI due to total aerosols (AOD) and dust (DOD) from both datasets of MIDAS and CAMS. The maximum of seasonal value is denoted as bold, emphasizing the peak of the seasonal cycle.

Daily variability
The day-to-day variability of GHI attenuation due to total aerosols for all pixels (Figure 10a), is much greater than the annual and seasonal values. The underestimation of CAMS AOD values is reflected in the constant lower values of their effects on daily GHI. There is no value of GHI reduction due to CAMS AOD above 50%, which is related to the lack of CAMS AOD values greater than 3 (see section 3.1). It is noteworthy that there are days that the aerosol effects on GHI reached values up to ~75% using the MIDAS AOD dataset. The daily values of GHI reduction due to dust using MIDAS DOD dataset are constantly greater that those when CAMS DOD is used (Figure 10b), with only exception the lower bin around zero. There are days for MIDAS dataset where the DOD GHI attenuation exceed 60%, while an upper limit of ~45% was resulted for CAMS. The strong impact of the aerosol particles on the direct component of solar radiation reaching the earths' surface is depicted in the distributions of the DNI attenuation due to both total aerosols ( Figure 10c) and dust (Figure 10d) with values up to 100% for MIDAS datasets, while for CAMS lower values were resulted. This intercomparison between MIDAS and CAMS AOD and DOD effects on DSSI showed how the AOD and DOD differences between the two datasets were expressed in differences in their radiative effects. Due to CAMS underestimation for high AODs and the strong DOD underestimation, the MIDAS dataset will be used and discussed in the subsequent analysis. The resulted radiative effects on surface indicated the important role of total aerosols and especially dust over MENA region, where DSSI attenuation by clouds is comparable or even lower than aerosols.

Extreme Dust events
The broader Mediterranean basin is an area frequently affected by dust outbreaks [36,65], resulting in extremely high concentrations of mineral particles and maximized DODs. Thus, we aimed to quantify the impact of those extreme DOD values on GHI and DNI. To this end, the methodology proposed by Gkikas et al. [63] was applied to MIDAS DOD values in order to define the extreme Dust Episode Day (eDED), by adapting also the objective and dynamic algorithm of Gkikas et al. [66] to MIDAS DODs. First, for every pixel the mean DOD ( ) and the associated standard deviation ( ) values were calculated, using the daily values of DOD over the time period 2003-2017. An extreme dust episode occurs on a specific day and location (pixel), when DOD values are higher than a critical value (threshold): This algorithm is characterized as dynamic, since the DOD threshold values are not constant for each pixel. Finally, in order to define a day as an eDED, at least 300 pixels should undergo an extreme dust episode, providing that the data availability for this day is more than 50%. From our analysis, 67 eDEDs were found for the whole study period 2003-2017, or on average 4.5 eDED yr -1 . Figure 11a presents the geographical distribution of mean MIDAS DOD from the resulted 67 eDEDs. The highest value up to 0.5 of mean DOD from the extremes were found over the northwestern Africa and lower mean values up to 0.43 and 0.35 were found over Egypt and the Middle East respectively. The highest values of the associated impacts on GHI and DNI up to 12% and 44% for GHI and DNI, respectively (Figures 11 b and 11 c), were found over the northwestern Africa, specifically over south Tunisia and the central Algeria. Large values of eDEDs mean attenuation were also found over Libya, Egypt and Middle East with values ranging from 4% to 9% for GHI and 17% to 35% for DNI for the bigger part of this areas. Cyprus is the Mediterranean island that was affected the most by the resulted extreme dust events with mean values of GHI and DNI attenuation up to 6.5% and 24%, respectively. For the south European countries, lower values of mean eDEDs attenuation were derived, up to 4% for GHI and 19% for DNI, with exception a maximum over southeastern Spain with values up to 5% and 23%, respectively.
It should be emphasized here that those results correspond to the long-term (2003-2017) average of the eDEDs and their radiative effects over the ROI. Individual dust events are associated with extremely high dust concentrations, resulted to GHI and DNI attenuations up to 50% and 90%, respectively [32].

Interannual variability
There is a global interest for the surface solar radiation changes due to their implications on climate and life on earth and it well documented a twenty year decrease worldwide until 1980 (dimming period) followed by an increase (brightening period) and those changes can be attributed primarily to aerosol's variability [67,68]. Using the simulated clear-sky DSSI values of this study based on highly quality aerosol inputs and for an aerosol free atmosphere (clean-sky), we aimed to investigate the clear-sky DSSI attenuation interannual variability attributed to total aerosols, but also to distinguish into the contribution of dust and non-dust component.
The analysis of the 15-year time series of annual GHI and DNI attenuation attributed to total aerosols and dust using MIDAS AOD and DOD datasets is presented in this section. Also, by subtracting DOD from AOD, the non-Dust Optical Depth (nDOD) was derived which is attributed to all other aerosol components besides dust. The main assumption here in order to derive this new nDOD is that the dust particles are externally mixed with the rest aerosol chemical species. For the RTM simulations, the same additional optical properties for the nDOD runs were assumed with those of total aerosols and using the RTM results without aerosols the corresponding annual GHI and DNI attenuations attributed to all other aerosol components except dust were derived (nDOD GHI and DNI attenuation hereafter). The interannual variability of the AOD (red line), DOD (blue line) and nDOD (green line) GHI and DNI attenuation is presented (Figure 12 b and c) for three different domain averages. The selection of the domains for the spatial averaging was based on the south to north gradient of dust and the geographical limits of those domains are illustrated in Figure 12a.
The year-to-year variability of GHI attenuation by the different aerosol components has lower values (0.5 to 1%) compared to the DNI attanuation (2-4.2%). For the DSSI attenuation (both GHI and DNI) by total aerosols a succesive decline was found after 2008, which is more prominent for D3. This reduction in the AOD DSSI attenuation is in line with the brightening effect over the Mediterranean [69,70]. The resulted decline in the DSSI attenuation bz aerosol is attributed mainly to the dust component for D1 and D2, where the variability of annual DOD DSSI attenuation is also great and highly correlated with annual attenuation values by total aerosols (correlation coefficients -cc ranging from 0.85 to 0.92). For D3, which has the sharpest decrease of DSSI attenuation by total aerosols, this is attributed to both dust (cc=0.84) and the other components (cc=0.87).

GHI and DNI clear-sky climatology
The availability of solar resources at the earth's surface is essential information for the different phases of a plant's deployment and operation. Average annual solar irradiation is a primary site selection criterion [27] and a low seasonal variability is preferable in order to match the power demand. According to the results of previous sections, there are areas that the GHI and DNI attenuation attributed to aerosols can reach values up to 13% and 50%, respectively, and those are areas with high solar energy potential (e.g., North Africa). So, the clear-sky GHI and DNI mean annual integrals based on high quality AOD retrievals are for great importance for such aerosol affected areas with scares cloudiness.
The clear-sky climatology of GHI and DNI was derived relying on MODIS AOD as input in the RTM. Using the daily irradiations (see section 2.2.3) annual and seasonal integrals of GHI and DNI were derived for every year and their mean values were calculated for the entire period (2003-2017) and presented in Figures 13 and S15-16. Given the fact that the cloud effect is missing from those figures, the description of the results is focused on the south part of the domain, over North Africa and Middle East, areas with high solar energy potential, scarce cloudiness and high aerosol loads.
The patterns of both GHI and DNI spatial variability are mainly driven by those of MODIS AOD (Section 3.1) and the surface altitude, and for GHI a latitudinal gradient (south-to-north) is evident as well. Over North Africa and Middle East the cumulative annual GHI and DNI range from 7500 to 8800 MJ/m 2 and from 7000 to 12000 MJ/m 2 , respectively ( Figure 13). Maximum values are observed in the Atlas Mountains (Northwest Africa), in the western parts of Libya and the southeastern parts of ROI. Regarding the spatial variability of the seasonal integrals, we will focus on spring and summer, when the effect of clouds is minimal for those latitudes. This time of the year, the distribution of aerosols expands to the western parts of North Africa and the spatial variability of DSSI components follows that pattern. While the maximum levels of DNI over the high-altitude areas are ~3300 MJ/m 2 in spring and ~3100 MJ/m 2 in summer, the minimum levels were found ~1800 MJ/m 2 for aerosol affected areas. Those differences are less pronounced for GHI (400 MJ/m 2 difference in both seasons).  Table 7. Change (in %) of the regional averaged mean annual and seasonal integrals of DNI due to total aerosols (AOD) and dust (DOD) from both datasets of MIDAS and CAMS. The maximum of seasonal value is denoted as bold, emphasizing the peak of the seasonal cycle.

Summary and conclusions
The broader Mediterranean basin hosts and receives a variety of aerosol types, which are quite variable in spatial and temporal terms. The overarching goal of the present study is to provide an insight of the perturbation of the surface solar radiation, and the subsequent impacts on solar energy production, attributed to the presence of all aerosols whereas focus is also given especially on dust. To realize, AOD/DOD from two different datasets (MIDAS and CAMS) were used as inputs to the LibRadtran RTM (in terms of precalculated LUTs) along with other critical parameters under clear-sky conditions such as additional aerosol optical properties of SSA and AE (MACv2 climatology) and key atmospheric parameters such as TCWV (CAMS reanalysis) and TOC (OMI-TOMS satellite retrievals). Model outputs were the GHI and DNI, which are of particular interest for different solar power systems (PV and CSP respectively). Our study domain encloses the broader Mediterranean basin and the study period spans from 2003 to 2017 (15 years).
Using the high-quality satellite derived MIDAS AOD/DOD datasets a 15-year climatology of total aerosols and dust was established for the broader Mediterranean basin. Based on our results, the largest AODs were found over dust sources or areas affected by dust transport, with maximum long-term averaged AOD up to 0.48 over Northwest Africa (up to 0.66 for summer season). Over the same area, the peak of MIDAS DOD values was derived as well, with mean annual value up to 0.35 (up to 0.45 for summer). Dust was found to contribute to total aerosol loads, in optical terms, from 40% to 90% over North Africa and Middle East, making dust the most important aerosol component over those areas.
Aerosols attenuate GHI by 1-13% and DNI by 5-47%. The larger values of their effects (4-13% for GHI and 15-47% for DNI) were found over North Africa and Middle East. Over the same area, the GHI and DNI reduction by dust ranges from 2-10% and 9-37%, respectively, contributing from 45-90% to total aerosol effects on DSSI. During the dry seasons of the year, when the cloud effects over those areas are comparable or even lower than aerosols, the maximum of aerosol and dust attenuation on GHI (up to 14% and 11%, respectively) and DNI (up to 53% and 40%, respectively) was found, with the contribution of dust to AOD DSSI attenuation reaching values up to 95%, making dust the most important regulator of DSSI levels for those regions. On a daily basis, the GHI reduction due to total and dust aerosol reached substantially higher values up to ~75% and ~60% respectively. There were days that the DNI component was totally blocked (-100%) under high aerosol and dust loads.
The investigation of the intra annual variability of the effects of aerosols and dust on GHI revealed, apart from their seasonal variations, the significant role of SSA on calculating aerosols radiative effects. The combination of low SSA values with considerable AODs/DODs resulted to the winter peak of AOD GHI attenuation which is reversed compared to the seasonal cycle of MIDAS AOD (peak in summer and low in winter). The same reasons explain the sharp peak in spring of DOD GHI attenuation.
The interannual variability of the aerosol and dust effects on DSSI attenuation was assessed for three domains. The radiative effects of all the other aerosol components except dust (nDOD) were also estimated for this analysis. After 2008, a successive decline of aerosol effects on DSSI was found for all domains, which was attributed mainly to the dust component for the south part of the domain.
Since it is well documented that the Mediterranean basin is frequently affected by dust intrusions, an assessment of the GHI and DNI attenuation was conducted for extreme dust events over the area. By adopting the methodology proposed by [63,66] to MIDAS DOD dataset, 67 eDEDs (4.5 eDEDs yr -1 ) were found over the study area for the time period 2003-2017. During those events the mean DOD reached values up to 0.50 (over North Africa) and the corresponding GHI and DNI attenuation were 12% and 44%, respectively. South Europe was also found to be affected by eDEDs, with the largest GHI and DNI attenuation values were resulted in southeastern Spain reaching 5% and 25%, respectively.
Taking advantage of the 15 years of high-quality daily satellite retrievals of AOD from MIDAS dataset, a clear-sky GHI and DNI climatology for the broader Mediterranean basin was derived. An added value of this new clear-sky climatology is that the DSSI values were simulated using, apart from satellite derived AOD, a climatology of additional aerosol optical properties (SSA, AE), model and satellite products for key atmospheric factors (TCWV and TOC). For the south part of the domain, in summer when the levels of DSSI are maximum, the main attenuator of GHI and DNI are aerosols (e.g. [30]). So, we focused on North Africa and Middle East and according to our findings for summer the spatial variability of GHI and DNI is 14% and 42% respectively.
The intercomparison between MIDAS and CAMS datasets revealed that CAMS slightly underestimates AOD and this is mainly evident over areas hosting major aerosol sources and strongly underestimates DOD up to 40% (-0.03) with respect to MIDAS DOD, which is in agreement with those reported by Bennouna et al. [60] versus ground-based retrievals. The CAMS underestimation of high AODs resulted weaker GHI and DNI attenuations on average by 1-4% and 4-11%, respectively. Likewise, due to the pronounced CAMS DOD underestimation, weaker attenuations were found (by 0.5-4% for GHI and 1-15% for DNI). This sensitivity analysis between the two datasets emphasizes on the importance of using reliable aerosol and dust optical properties to accurately simulate DSSI.
Concluding, this study aims to contribute towards the better understanding of the role of aerosols and especially of dust on surface solar radiation in terms of GHI and DNI over the Mediterranean basin. The results of this analysis, apart from their importance from the perspective of climate science, provides valuable information in terms of management and future planning of PV and CSP installations.
Supplementary Materials: The following supporting information can be downloaded at: www.mdpi.com/xxx/s1, Figure S1: Geographical distribution of long term average of seasonal mean AOD at 550nm from MODIS. Blank grid points are those that didn't fulfill the criterion of at least 20% data availability on annual basis.; Figure S2: Geographical distribution of long term average of seasonal mean AOD at 550nm from CAMS. Blank grid points are those that didn't fulfill the criterion of at least 20% data availability on annual basis. Figure S3: Geographical distribution of annual mean CAMS-MIDAS AOD biases. Blank grid points are those that didn't fulfill the criterion of at least 20% data availability on annual basis. Figure S4: Geographical distribution of long term average of seasonal mean DOD at 550nm from MIDAS. Blank grid points are those that didn't fulfill the criterion of at least 20% data availability on annual basis. Figure S5: Geographical distribution of long term average of seasonal mean DOD at 550nm from CAMS. Blank grid points are those that didn't fulfill the criterion of at least 20% data availability on annual basis. Figure S6: Geographical distribution of annual mean CAMS-MIDAS DOD biases. Blank grid points are those that didn't fulfill the criterion of at least 20% data availability on annual basis. Figure S7: Change (in %) of the mean seasonal integral of GHI due to the presence of aerosols under MODIS AOD. Blank grid points are those that didn't fulfill the criterion of at least 20% data availability on annual basis. Figure S8: Change (in %) of the mean seasonal integral of GHI due to the presence of aerosols under CAMS AOD. Blank grid points are those that didn't fulfill the criterion of at least 20% data availability on annual basis. Figure S9: Change (in %) of the mean annual integral of GHI due to the presence of dust under MIDAS DOD. Blank grid points are those that didn't fulfill the criterion of at least 20% data availability on annual basis. Figure S10: Change (in %) of the mean annual integral of GHI due to the presence of dust under CAMS DOD. Blank grid points are those that didn't fulfill the criterion of at least 20% data availability on annual basis. Figure S11: Geographical distribution of seasonal mean SSA (MACv2, [48]). Figure S12: Geographical distribution of seasonal mean DU SSA (MACv2, [48]). Figure S13: Change (in %) of the mean annual integral of DNI due to the presence of aerosols under MODIS AOD. Blank grid points are those that didn't fulfill the criterion of at least 20% data availability on annual basis. Figure S14: Change (in %) of the mean annual integral of DNI due to the presence of aerosols under CAMS AOD. Blank grid points are those that didn't fulfill the criterion of at least 20% data availability on annual basis. Figure S15: Change (in %) of the mean annual integral of DNI due to the presence of dust under MIDAS DOD. Blank grid points are those that didn't fulfill the criterion of at least 20% data availability on annual basis. Figure S16: Change (in %) of the mean annual integral of DNI due to the presence of dust under CAMS DOD. Blank grid points are those that didn't fulfill the criterion of at least 20% data availability on annual basis. Figure S17: Mean seasonal integrals for clear -sky GHI using MODIS AOD. Figure S18: Mean seasonal integrals for clear -sky DNI using MODIS AOD.
Funding: This study has been funded by European Commission project EuroGEO e-shape (grant agreement no. 820852).