Effects of Aerosols and Clouds on the Levels of Surface Solar Radiation and Solar Energy in Cyprus

: Cyprus plans to drastically increase the share of renewable energy sources from 13.9% in 2020 to 22.9% in 2030. Solar energy can play a key role in the effort to fulﬁl this goal. The potential for production of solar energy over the island is much higher than most of European territory because of the low latitude of the island and the nearly cloudless summers. In this study, high quality and ﬁne resolution satellite retrievals of aerosols and dust, from the newly developed MIDAS climatology, and information for clouds from CM SAF are used in order to quantify the effects of aerosols, dust, and clouds on the levels of surface solar radiation for 2004–2017 and the corresponding ﬁnancial loss for different types of installations for the production of solar energy. Surface solar radiation climatology has also been developed based on the above information. Ground-based measurements were also incorporated to study the contribution of different species to the aerosol mixture and the effects of day-to-day variability of aerosols on SSR. Aerosols attenuate 5–10% of the annual global horizontal irradiation and 15–35% of the annual direct normal irradiation, while clouds attenuate 25–30% and 35–50% respectively. Dust is responsible for 30–50% of the overall attenuation by aerosols and is the main regulator of the variability of total aerosol. All-sky annual global horizontal irradiation increased signiﬁcantly in the period of study by 2%, which was mainly attributed to changes in cloudiness.


Introduction
Mitigation of climate change is one of the main challenges of humanity for the 21st century [1,2]. It is currently accepted from the vast majority of the scientific community that since the mid 20th century global average temperature increases quickly mainly due to increasing anthropogenic greenhouse gas (GHG) concentrations [3]. The energy sector and optical properties are generally available from only few, sparsely distributed monitoring stations (globally, and in Cyprus), only satellite retrievals can provide such information on a wide spatial scale. Aerosol and cloud optical properties are currently available in high spatial and temporal resolution from, e.g., the moderate resolution imaging spectrometer (MODIS) [48] and the satellite application facilities on climate monitoring (CM SAF).
The present study aims to highlight the role of aerosols, with special focus on dust aerosols, in the production of solar energy in Cyprus. The recently developed high quality and high resolution ModIs Dust AeroSol (MIDAS) dust climatology [49] was used for this purpose. The spatial and the temporal variability of GHI (of interest mainly for PV) and direct normal irradiance (DNI) (of interest mainly for CSP) were studied with respect to the variability of aerosols and clouds. Using high quality satellite aerosol and cloud products 14-year climatologies of the DNI and GHI were created. The aerosol, cloudiness and SSR trends over Cyprus were studied for the same period. The advantage of the new climatology compared to existing climatologies (e.g., CM SAF) is that SSR has been simulated using satellite high resolution retrievals instead of climatological or assimilated aerosol optical depth (AOD) (e.g., [50]).
The paper is divided as follows. The methods and the datasets used for the analyses are discussed in Section 2. Results are presented in Section 3, and the main conclusions of the study are summarized in Section 4. Section 3 is divided in five subsections. In Section 3.1 the aerosol optical properties and the composition of aerosol mixture over Cyprus are discussed. In Section 3.2 the effects of aerosols, dust, and clouds on GHI and DNI are discussed. GHI and DNI climatologies are presented in Section 3.3. In Section 3.4 an example of the effect of a strong dust event on the levels of GHI and DNI is discussed. In Section 3.5 the economic impact of the attenuation by clouds, dust and total aerosols on different types of PV and CSP installations is discussed.

Ground Based and Satellite Datasets
Aerosol optical depth (AOD) and dust optical depth (DOD) from the MIDAS dataset [49] were used in the present study. MIDAS [49] AOD and DOD at 550 nm are available at a global scale and fine spatial resolution (0.1 • × 0.1 • ), over a 15-year period (2003-2017). For the former parameter, the MODIS onboard Aqua satellite (MODIS-Aqua) (moderate resolution imaging spectroradiometer) [51] AOD Level 2 (L2) retrievals from Collection 6.1 (C061; [48]) were processed by applying quality control screening. The MIDAS DOD constitutes a product developed by multiplying MODIS-Aqua AOD with the DOD/AOD ratio, derived by the Modern-Era Retrospective Analysis for Research and Applications version 2 (MERRA-2; [52]).
The expected errors (EEs) in the MODIS AOD product [53] used for the development of the MIDAS climatology are about ±0.15 × AOD + 0.05 over land and ±0.05 × AOD +0.03 over ocean [54]. The reliability of the MIDAS DOD was justified via its comparison against the aerosol robotic network (AERONET, [55]) AODs. Moreover, an intercomparison study of the MIDAS, MERRA-2 and LIVAS [56] DODs has been performed revealing a moderateto-high consistency among the three datasets depending on the region of interest. Overall, both assessment analyses have confirmed the quality of the MIDAS DOD thus making meaningful its exploitation in radiation studies. The same MODIS AOD product, which was used for the development of MIDAS climatology, was also used in the present study. Hereon it is referred as MODIS AOD.
Attenuation of the monthly integrals of SSR by clouds was calculated from CM SAF Surface Radiation Data Set-Heliosat (SARAH)-Edition 2.1 [57] retrievals of the SSR (hereon this product is referred as CM SAF-SARAH2.1). Reflection measurements from instruments on-board the Meteosat satellite missions [48] were used as inputs in order to simulate the SSR data record under all-sky conditions. More details for this simulation approach can be found in [58]. Climatological water vapor and ozone from reanalysis data have been incorporated in the algorithm and assimilated climatological aerosol optical properties [59]. Uncertainties in the retrieved monthly CMF were estimated to 3% as discussed in Section 2.3.
The lidar observations in Limassol (34.7 • N, 33.0 • E; 10 m above sea level, a.s.l.), Cyprus, were conducted in the framework of the Cyprus, Cloud, Aerosol and Rain Experiment (Cy-CARE) campaign lasting from October 2016 to the end of March 2018. The Polly instrument (portable Lidar system) [60] is the key instrument of LACROS (Leipzig Aerosol and Cloud Remote Observation System; http://lacros.rsd.tropos.de/, last access: 10 April 2021) [61,62] for aerosol profiling, and it is installed inside a container. This lidar has 13 channels and continuously measures elastic and Raman backscatter signals at the laser wavelengths of 355, 532 and 1064 nm and respective Raman backscattering wavelengths of 387 and 607 nm (nitrogen Raman scattering) and 407 nm (water vapor Raman scattering) [63,64].
In addition, a CIMEL sun photometer is operating at Cyprus University of Technology (CUT-TEPAK) in Limassol, since 2010, as part of AERONET [55]. The photometer performs direct and sky radiance measurements at nine band passes. Direct measurements are performed every 10-15 min (under cloudless skies) and AOD is retrieved at the corresponding wavelengths. Sky radiance measurements are performed less frequently, at specific SZA, and are fed in an inversion radiative transfer algorithm in order to retrieve other optical properties of aerosols [55]. In this study version 3 AERONET retrievals [65] were used. AERONET data are classified at three quality levels: level 1.0 (unscreened), level 1.5 (automatic cloud screening) and level 2.0 (quality assured with post field calibration). For CUT-TEPAK, currently level 2.0 data are available up to 2017. Timeseries of level 2.0 for 2010-2017 and level 1.5 for 2017-2020 were analyzed. Although the uncertainties in level 1.5 were larger than the level 2.0 uncertainties, the climatological analysis performed in this study and comparisons with other data sources (MODIS) did not result in significant differences between the two datasets.

Data Processing
In Figure 1, histograms of differences between MODIS AOD at 550 nm with collocated AERONET retrievals are shown. AOD at 550 nm from AERONET was calculated using the 440-870 nm Ångström Exponent (AE) and AOD at 500 nm. Daily averages were calculated when at least 16 measurements were available in the day and were compared with MODIS AOD (measured on a daily basis, at the time of MODIS-Aqua overpasses). The AERONET-MIDAS comparison resulted to a mean difference of 0.02 with standard deviation (σ) = 0.09, which is within the uncertainty of the AERONET product [66], and of similar magnitude as the average agreement between AOD measurements from different colocated ground based instruments [67]. The distribution shown in Figure 1 is slightly skewed towards MODIS overestimation. Comparison between AOD from overpasses (MODIS) and daily averages (AERONET) is possibly responsible for part of the differences. Overestimation of AOD by MODIS may be related with the uncertainties of MODIS retrievals in coastal areas, where pixels include both land and ocean [48,68]. Nevertheless, the agreement between satellite (MODIS) and ground based (AERONET) data sets is satisfactory since it is within the combined uncertainties of the two datasets. Since long-term continuous AERONET retrievals in Cyprus are available only for Limassol, AOD and DOD at 550 nm (hereon AOD and DOD refer to AOD and DOD at 550 nm unless something else is specified) from MIDAS were used for the study of aerosol variability on a wider spatial scale. AOD and DOD corresponding to the five locations shown in Figure 2 were calculated as averages of the nearest four satellite pixels for each location.  Monthly averages of AOD and DOD for each of the five locations were calculated for months for which overpasses were available for more than 6 days.
The calculations of cloudless and aerosol-free sky (clean-sky) SSR in terms of global (I g ), diffuse (I s ), and direct beam (I d ) components were performed using precalculated look-up tables (LUT), for improving the computational time similar to [69]. All quantities discussed in the study were then calculated from direct and diffuse components with the assumption that diffuse irradiance is homogeneous. The following quantities were calculated: Solar irradiance (direct + diffuse) at an inclined horizontal surface with inclination angle equal to the latitude of the location; • Solar irradiance (direct + diffuse) at an inclined horizontal surface with inclination angle equal to the latitude of the location, which follows solar azimuth; • Solar irradiance (direct + diffuse) at a surface that is constantly perpendicular to the solar beam.
LUT contain SSR simulated using the LibRadtran radiative transfer model (RTM) [70] for different SZAs and a wide range of atmospheric parameters, which attenuate SSR for cloudless conditions,. Using the radiative transfer solver sdisort [71] pseudospectral simulations were performed with a spectral resolution of 1 nm in the range of 280-3000 nm, using the parameterization of LOWTRAN band model [72] as adopted from the SBDART algorithm [73] to describe molecular absorption.
The Kurucz [74] extraterrestrial solar spectrum (with a resolution of 1.0 nm) and the US standard atmospheric profile [75] were used for the simulations. Surface albedo was set to 0.2. Cloudless-sky SSR was finally simulated for all possible combinations of the input parameters presented in Table 1. The output spectra were integrated with respect to wavelength, and the short wave (SW) irradiances were derived.
The simulated SSR was subtracted from the LUT for the desired conditions using linear interpolation at multiple dimensions. In order to quantify the effects of total aerosol and dust on SSR, the irradiance was extracted from the LUT for AOD and DOD from MIDAS, and for zero AOD and DOD (clean-sky) with a temporal resolution of 1 h, for the period 2004-2017. Average annual climatological TOC (= 310 DU) and TCWV (= 1.8 cm) for Cyprus were used for the retrieval of SSR. Climatological TOC was derived from daily TOC for 2004-2017 from the ozone monitoring instrument (OMI). Annual average of the TCWV was calculated from 3 h values from the Copernicus Atmospheric Monitoring Service (CAMS) [76] for the same period. Climatological values of 0.92 and 1.3 were used for total aerosol SSA at 550 nm and AE 440-675 nm respectively, as they were estimated by Taylor et al. [77], for the area of study. When SSR was subtracted for DOD the corresponding values for SSA and AE were 0.93 and 0.6 respectively. In a similar way, SSR was subtracted from the LUT for Limassol using aerosol optical properties (AOD, SSA, and AE) from the CIMEL.
In case of missing AOD or DOD values, the corresponding averages for the particular month were used for the simulations, unless retrievals were available for less than 20% of the month days. In the latter case seasonal average values were used. If, however, retrievals were available for less of 20% of the month days for all three months of the season, then bilinear spatial interpolation between the nearest pixels for which data were available was performed.
For the study of the effects of dust and total aerosols the ratios between the monthly integrals of SSR for cloudless-skies (the MIDAS AOD or DOD at 550 nm was used for the simulations) and clean-skies were compared to each other.
Attenuation of the monthly integrals of SSR by clouds was calculated from CM SAF-SARAH2.1 [57] retrievals for 2004-2017. The ratio between all-sky and cloudless-sky irradiances, commonly called the cloud modification factor (CMF), was calculated for the diffuse and the direct components of the SSR. Then, the monthly CMF was calculated for each quantity (GHI, DNI, etc.) depending on the contribution of each component. Monthly integrals of the cloudless-sky irradiance were multiplied with the corresponding monthly CMF for the retrieval of the all-sky irradiance. Monthly climatologies of the all-sky irradiances (averages of all months for 2004-2017) were finally calculated for Cyprus by averaging over the whole 14-year period.
Analysis using RTM simulations showed that for usual aerosol conditions the GHI and DNI increase with altitude by 1-3%/km, which is in agreement with the results of similar analyses in other studies [78]. Thus, an average increase of 2%/km was assumed and the whole dataset of the monthly all-sky irradiances was post-corrected for the effect of altitude. Active remote sensing instruments such as PollyXT lidars are a key technique for characterizing aerosols as they can provide vertically resolved information on extensive (e.g., aerosol backscatter coefficient, aerosol extinction coefficient and volume depolarization ratio) and intensive properties (e.g., Ångström exponent, lidar ratio and particle depolarization ratio) of different aerosol types. The extensive properties depend on the aerosol concentration, whilst intensive ones are type-sensitive and provide separate classification for each detected layer (e.g., [79][80][81][82][83][84][85]).
Moreover, a careful investigation of the air mass origin and dust transport path was performed by means of backward trajectory analysis. This analysis was carried out using the hybrid single-particle Lagrangian integrated trajectory (HYSPLIT) model (https://ready.arl.noaa.gov/HYSPLIT_traj.php, last access: 1 April 2021; [86]) together with the global data analysis system (GDAS) meteorological files (spatial resolution of 1 • × 1 • , every 3 h) as data input. The kinematic back-trajectories were calculated using the vertical velocity component given by the meteorological model with a 96-120 h pathway (4-5 d back).
In the context of this study, radiometric measurements performed at MORDOR were also analyzed. Daily integrals of GHI and DNI were calculated and compared with the corresponding simulated quantities for Limassol.

Uncertainty in the Simulation of GHI and DNI
The uncertainty in radiative transfer simulations depends mainly on the uncertainty of the input parameters to the RTM. Uncertainties in the CM SAF-SARAH2.1 product are mainly determined by uncertainties in the determination of the clear-sky reflection and the self-calibration method [59]. The former uncertainties are significant mainly over areas with long-lasting cloudiness, which is not the case for Cyprus. The latter are of the order of 3% (standard uncertainty) for the monthly average values. As discussed in Section 2.1 the errors in the AOD from MODIS are generally estimated to be of the order of 0.05-0.06 for usual AOD levels at Cyprus. The corresponding error in the simulated GHI is estimated to be less than 2% (standard uncertainty of 1%) for GHI and 8% (standard uncertainty of 5%) for DNI. A standard surface albedo of 0.2 was used as input to the RTM. Assuming that surface albedo may vary spatially and temporally by up to about ± 0.15 (e.g., [87]) in Cyprus (maybe with the exception of very high altitude areas in winter where snowfall, and thus albedo above 0.35, is possible), the corresponding standard uncertainty in the simulated GHI may be of the order of 2%. Annual average climatological values of the TCWV and TOC were also assumed. On a seasonal scale, TCWV ranges from 1.3 cm (in winter) to 2.2 cm (in summer), while the annual average (which was used for the simulations) is 1.8 cm. These differences in TCWV however have a small effect on GHI and DNI (of 1%). Finally, the seasonal variability in TOC is of the order of ± 30 DU around the annual average, which again induces negligible uncertainties, below 1%. Finally, it has to be noted that the simulations are for open horizon. Thus, obstacles blocking the direct solar beam, or even part of the diffuse radiation may induce differences between the simulated and the corresponding real radiometric quantities.

Aerosols in Cyprus
A climatology of direct sun retrievals is presented in Table 2. AOD was highest in May-August and lowest during January and December, which is similar to the pattern reported for other Eastern Mediterranean sites as well [26,27,88,89]. AE follows a more complicated pattern with highest values in August and lowest in April and May. AE is associated with particles' size, with higher values indicating more small particles in the mixture. This seasonality suggests that mixtures containing large particles were more frequent in spring. The seasonal mean size distribution of columnar aerosol mixtures in the area was calculated using the inversion product of AERONET and is presented in Figure 3. It is clear, that large particles were dominant in aerosol mixtures during spring, which denotes significant contribution from dust aerosols. In winter, aerosol mixtures generally contained fewer large particles relative to other seasons. Small particles dominated the mixtures during December-February, which were usually linked with anthropogenic activities and combustion particularly. We should also notice that the period when anthropogenic aerosols played the greatest role was when AOD had the lowest values. Table 2. Climatology of AOD and AE as derived from level 2.0 and level 1.5 version 3 AERONET direct sun product for CUT-TEPAK for the 2010-2020 period. In the last three columns there are the number of measurements (# Data), days for which measurements were available (# Days), and months for which measurements were available (# Months).  In Figure 4, the commonly used aerosol classification following the method of Dubovik et al. [90] is presented. This method was structured using the climatology of AOD and Ångström exponent of AERONET stations where the dominant aerosol type is well known, to set threshold values for classifying the mixtures. The drawback of this approach is that these threshold values could differ in cases of complicated aerosol mixtures. In urban areas, it is very rare to have a mixture where only a single type of aerosols could be identified. Hence, the classification should be considered as an identification of the dominant type of particles that have optical properties that are closest to the measured ones. Utilizing this scheme, we found that 25% of the recordings were dominated by dust particles and also a considerable part is characterized as mixed (co-occurring big and polluted particles). Dust particles in the area could be originating either from the Saharan Desert or Arabian Peninsula. Marine and continental aerosols were found in 13% and 16% of the cases, which were associated with very low AOD values. Biomass burning type was rare in the area (2% of the cases) but was present when AOD was high, which translated to a relative significant radiative effect. Aerosol properties were available from CIMEL (on a regular basis and on a long-term period) only for Limassol. Thus, the annual evolution of AOD at 550 nm, at five different locations in Cyprus ( Figure 2) was studied using the MIDAS product. The results are presented in Figure 5. Letters in the horizontal axis represent the twelve months of the year from January to December.

Month AOD
The results obtained from MIDAS were not directly comparable with those discussed in Table 2 because they referred to different periods. Analysis of the series for the five stations revealed that the AOD changed in a statistically significant way (at the 95% confidence level) in 2004-2017 for some of the stations (see Table 3). For example, AOD at Limassol was decreasing with an average rate of 0.011/year in spring (overall decrease of 0.15 in the 14-year period). It is interesting that in all three cases for which statistically significant trends in AOD were found, the same trends were also found for DOD, which denotes that changes in AOD are fully attributable to changes in DOD.  In all cases presented in Figure 5, the contribution of dust to the aerosol mixture was maximum in spring (which agreed with the results of Figure 3), and there was also a secondary maximum in autumn. Average DOD in autumn (0.05-0.1 depending on the station) was however about half of the corresponding average DOD in spring (0.15-0.2). The AOD was maximum in spring and summer with values around 0.3 and 0.25 respectively. Maximal differences between the AOD at the five locations were found for autumn, when AOD at Nicosia was generally around 0.1 while AOD at Limassol was 0.15-0.25.

Effect of Aerosols and Clouds on SSR
In Figure 6, the overall effects of aerosols, dust, and clouds on annual GHI and DNI are presented. The effects of aerosols and dust were quantified for MODIS AOD and MIDAS DOD, while the effects of clouds were quantified for CMF from CM SAF-SARAH2.1 (see Section 2). Aerosols were found to attenuate 5-10% of GHI, with about 30-40% of the overall attenuation by aerosols to be due to dust. Attenuation of DNI by aerosols was stronger, between 15 and 35%, and dust attenuated 10-17% of DNI (30-50% of the overall attenuation by aerosols). Attenuation of GHI and DNI by aerosols was minimum at the north of Troodos Mountains, with the differences between the north and the south of the mountains being more pronounced for DNI. The differences between the north and the south of the mountains were possibly due to the effect of the mountains to the weather patterns, and subsequently to aerosol over the area. Further analysis, which was out of the scope of the present study, is necessary in order to determine the exact reasons that are responsible for the weaker attenuation by aerosols at the north of the Mountains. Clouds attenuate 12-14% of GHI and 22-25% of DNI. The spatial variability of their effect was small with maxima over the Troodos Mountains and locally at the north shores of the island.
The annual variability of the effects of aerosols and clouds on GHI and DNI was investigated for the same five locations as in Section 3.1 and the results are presented in Figure 7. Shaded areas in the graphs represent the standard deviation of the climatological averages for the 14-year period. As expected, in Omodos, which is at high altitude at the south side of Troodos the effects of clouds are generally the strongest among the five locations. What is interesting is that attenuation by aerosols is also strong at Troodos despite the high altitude. For all five locations there is a clear annual cycle of the attenuation by clouds. In winter, clouds decrease GHI by 25-30% and DNI by 35-50%. In spring and autumn attenuation by clouds becomes weaker. In summer, the role of clouds in the attenuation of GHI was practically the same with that of aerosols (also see Table 4), while the attenuation of DNI by clouds was 7-8 times weaker compared to the attenuation by aerosols. Dust blocked on average 3 times more DNI than clouds in summer.
The median, the minimum, and the maximum attenuation of the annual GHI and DNI for the 14-year period are presented in Table 4. The same quantities for summer are provided in parentheses. For the four lower altitude locations the effect of clouds on annual and summer GHI and DNI was found to be relatively stable in the 14-year period, changing by less than ±1% and ±2% respectively in the years. At Omodos, the variability was slightly larger, and, in specific years, attenuation of annual GHI and DNI by clouds was 3% larger than the corresponding medians.
Attenuation of GHI by aerosols (and dust) was constant within 3-4% for the annual and the summer integrals for all five sites. Although attenuation of GHI by aerosols did not change significantly from year to year, this was not the case for DNI. The largest variability in the attenuation of annual DNI by dust was found for Nicosia, where the difference between minimum and maximum attenuation was 8%. The corresponding difference for total aerosol was similar (9%). For the summer integrals the corresponding differences were even larger, 9% and 13% respectively. The results presented in Figure 7 and Table 4 show that the variability in the levels of GHI and DNI in Cyprus was strongly correlated with the variability of aerosols, and especially dust aerosols. It is noteworthy that low DOD coincided with low AOD, and high DOD coincided with high AOD for all five stations, which shows that dust was the most significant regulator of the levels of GHI and DNI.
Analysis of the 14-year series of the monthly and annual attenuation by aerosols and dust did not yield any statistically significant trends (at the 95% confidence level) for any of the five sites (despite the significant decrease of DOD and AOD in some cases). In the case of clouds statistically significant trends were found for some of the stations for specific months. The most interesting finding of this analysis was that in February, a statistically significant negative trend of the attenuation by clouds was found for all sites, for GHI and DNI (see Table A1 in the Appendix A). Depending on the quantity (GHI or DNI) and the site the decrease was 0.7-1.5% per year, which means an overall decrease of the attenuation of SSR by clouds of 10-20% in the 14-year period. Negative trends in the attenuation of annual GHI by clouds were found for all stations but were statistically significant only for Paphos. It is interesting that the corresponding trends in the attenuation of annual DNI were in all cases positive and not statistically significant. More detailed information regarding the changes in attenuation by clouds can be found in Table A1 in the Appendix A.
Combined effects of aerosols and clouds led to a statistically significant average increase (i.e., brightening effect) of the all-sky annual GHI (average for the five stations) of 0.13% ± 0.06% per year in 2004-2017. No statistically significant change of DNI was found for the same period. The latter results are in agreement with the results of recent studies which report solar brightening (i.e., increase of the GHI) over different Mediterranean sites in the last two decades mainly as a result of decreasing attenuation by clouds and/or aerosols [91][92][93][94].

Climatology
Climatology of the GHI and DNI was developed using SSR simulations based on MODIS AOD and CMF from CM SAF-SARAH2.1. The seasonal integrals (i.e., cumulative GHI and DNI) for the four seasons of the year are presented in Figure 8. Despite the small geographical extent of Cyprus, Troodos Mountains seem to have drastic effect on the spatial variability of solar potential, especially for DNI.
In winter there were more clouds at the top and the south of Troodos relative to the north. Furthermore, there were less aerosols at the north. These differences are depicted very clearly in the levels of DNI for which minimum levels (1100 MJ/m 2 ) were found at the south slopes of the mountains and maximum levels (1300 MJ/m 2 ) were found at the north, near the zone that separates the Greek Cypriot and Turkish Cypriot communities. Differences were less pronounced for GHI (minimum was 800 MJ/m 2 and maximum was 900 MJ/m 2 ). The pattern of spatial variability in autumn was similar to the winter pattern. GHI differed by up to 100 MJ/m 2 between the south and the north of the mountains, while the corresponding difference for DNI was 300 MJ/m 2 .
In spring and summer, the effect of clouds was less pronounced, and the spatial variability depended mainly on the distribution of aerosols over the island. Maximum levels of GHI (2100 MJ/m 2 in spring and 2700 MJ/m 2 in summer) were found at the north side of the mountains while minimum levels were found at the south (2000 MJ/m 2 in spring and 2600 MJ/m 2 in summer). Again, the differences were more pronounced for DNI. In both seasons the difference between the DNI north and south of the mountains may be up to 300 MJ/m 2 .
The maps of the annual integrals of GHI and DNI are shown in Figure 9. The cumulative GHI is 6800-7200 MJ/m 2 with maximum values at the north slopes of Troodos Mountains. The cumulative DNI was 7500-8500 MJ/m 2 , with maximum levels slightly shifted to the north relative to GHI.
At this point it must be noted that the standard deviation in the annual GHI was 120-160 MJ/m 2 , with maximum values over areas where cloudiness played a more significant role (i.e., the south slopes and the higher altitude regions of Troodos Mountains). For DNI the standard deviation was 250-300 MJ/m 2 . In both cases the standard deviation was less than half of the maximum differences shown in Figure 9. Furthermore, analysis of annual GHI and DNI for each year of the 14-year period revealed that although the average levels of irradiance differed from year to year, the pattern of the spatial distribution of both quantities was similar to that shown in Figure 9. Thus, the results presented in Figure 9 were considered reliable.  Although detailed comparison between SSR simulated using MODIS aerosols and SSR from CMSAF was out of the scope of the present study, the GHI and DNI from the two climatologies were compared for the five locations shown in Figure 2 in order to provide an estimate of the magnitude of the differences between the two datasets. Comparison between the two datasets revealed a relatively good agreement at lower altitude areas and larger differences at higher altitudes. Examples of the differences for five locations are shown in Table 5. At the four lower altitude sites the differences were generally small, and within the uncertainties of the simulations. Larger differences between SSR from the two datasets were found for the higher altitude site of Omodos ( 800 m a.s.l.). Since the effect of clouds was the same for the two datasets, differences were mainly due to differences in the AOD used for simulations, and the way that altitude was taken into account in the two datasets. Interpolation from different grid points over complex terrain/fast changing altitude sites might also lead to large differences, which could have an effect in the case of Omodos. A study where both datasets will be validated against ground-based measurements is planned for the future and will give more answers regarding the observed differences. In a recent study, comparison between ground based measurements of the GHI and the corresponding product from CM SAF-SARAH2 for the Iberian Peninsula yielded a relatively good agreement (R 0.83) [95].

Intense Aerosol Events
As discussed in previous sections aerosols play a key role in regulating the levels of SSR in Cyprus. Aerosols can vary significantly from day to day and extreme events (e.g., dust storms) can strongly reduce SSR relative to its climatological levels. In Figure 10 there is an example of an event that lasted 6 days at the station of CUT-TEPAK at Limassol. In particular, the vertical distribution of aerosol attenuated backscatter coefficient and of the particle linear depolarization ratio with respect to time of the day, as recorded by the lidar system, is presented. The upper panels (a and b) are for 10 September 2017 (DOY 152) when aerosol load was relatively low (average AOD at 500 nm is about 0.11). The middle panels (c and d) are for 15 September 2017 (DOY 157). The thick dust layer originated by the Middle East regions can be noticed between 1500 and 6000 m. The lower panels are for 21 September 2017 (DOY 163) when a thick layer of dust originated by the Northern Africa can be seen at the lower 1 km of the atmosphere up to 6 km. The dust particles as shown by the particle depolarization ratio (Figure 11f) reached the surface after 14:00 UTC.
The origin and the evolution of the events presented in Figure 10 were identified by processing measurements from the lidar [36,[96][97][98] and careful investigation of the air mass origin and long-range aerosol transport by means of backward trajectory analysis from the HYSPLIT (hybrid single particle Lagrangian integrated trajectory) model [86,99,100]. Access is provided via the NOAA-ARL Website (http://www.arl.noaa.gov/HYSPLIT.php (accessed on 11 June 2021)). HYSPLIT backward trajectory analysis together with the profiles of the particle depolarization ratio was used to identify the reference case and the cases with long-range transport of desert dust from the Middle East deserts and from Northern Africa (Sahara region).
Measurements of aerosol optical properties (AOD at 500 nm, 440-870 nm AE, and 440 nm SSA) from the CIMEL sun photometer were used as libRadtran inputs in order to simulate SSR for cloudless skies, instead of MODIS AOD and climatological values for SSA and AE. CIMEL measurements were chosen instead of MODIS retrievals because they were more representative for the station and they were available in a finer temporal resolution, although (as discussed in Section 2) the agreement between CIMEL daily averages and MODIS AOD is good.  (a) Daily modeled GHI and DNI for clean (aerosol-and cloud-free)-skies (GHI clean and DNI clean ), for cloudless-skies (GHI cimel and DNI cimel ), and measured GHI and DNI (GHI mordor and DNI mordor ), (b) ratios between the GHI and DNI for cloudless-and all-skies (modeled and measured respectively) with the corresponding aerosol-free cases, (c) differences between the modeled cloudfree daily GHI and DNI from the corresponding clean-sky quantities, and differences between the measured daily GHI and DNI and the corresponding clean-sky quantities, and (d) daily average AOD at 500 nm from CIMEL (blue line) and daily average AOD multiplied with the average daily air mass (purple line). The graphs are for DOY 153-305 (June-October). Figure 11 shows that during the dust event, which took place in September (between DOY 157 and 163), daily cloudless-sky GHI (GHI cimel ) was 10-12% below clean-sky GHI (GHI clean ), and cloudless-sky DNI (DNI cimel ) was 40-50% below clean-sky DNI (DNI clean ), corresponding to energy losses of 2.5-3 MJ/m 2 for GHI and 16-18 MJ/m 2 for DNI. Ac-cording to the results presented in Section 3.2, average reduction of the GHI and DNI at Limassol in September due to the presence of aerosols is 7% and 25% respectively.
Measurements of the GHI and the DNI, which were also performed in the CUT-TEPAK facilities during the same period, were added to the graph in order to show the effect of clouds and additionally to verify the validity of the simulations. Comparison of the measurements with the simulations (which include the effect of aerosols) shows that there was an offset between them, even for completely cloudless days. The offset was maximum in July ( 6% for GHI and 4% for DNI) and then gradually decreased to less than 1% in October. The reason is that there was an obstacle at the west of the sensor blocking sunlight in the evening. In June-July all DNI was blocked at SZA larger than 65 • . Gradually the SZA above which the DNI was blocked became larger, exceeding 80 • in October. Thus, in October the differences between the measured and the modeled daily GHI and DNI in cloud-free days were smaller than 1%. From the comparison of the measured and the modeled irradiances it is quite clear again that the role of clouds was secondary relative to the role of aerosols, especially for DNI.

Economical Impact of the Attenuation of SSR by Clouds, Aerosols, and Dust
In this section, financial analysis of the effects of clouds, aerosols, and dust on different PV installations was performed. The methodology for the analysis is the same as in Kosmopoulos et al. [27]. Analyses were performed for five different types of installations: (a) Horizontal PV panels (e.g., on the terrace of a building); (b) Tilted PV panels with tilt angle equal to the latitude of the site; (c) One axis solar tracking PV system (following solar azimuth); (d) Two axis solar tracking PV system (following solar zenith and azimuth); (e) CSP installation.
In all cases/simulations the nominal power of the installations was assumed to be 500 KW followed by a representative feed-in tariff of 0.0741 EUR/kWh [101]. The solar energy received from each of these installations was calculated as described in Section 2. The economic and energy impact were quantified in terms of monthly means, total energy loss (EL), total financial loss (FL), and solar energy potential. For the PV calculations (horizontal, tilted, and with 1-axis and 2-axis trackers), we assumed the crystalline silicon as the panels' material followed by combined losses of 29%, an efficiency value of 12% and a spatial coverage of 80% [24]. For the CSP simulations we took into account the effect of the heat transfer and the cumulative losses of shading, peak optical efficiency, and incident angle modifier [102]. Calculations were performed for the five locations discussed at Section 2 (Larnaca, Paphos, Limassol, Nicosia, and Omodos) in a similar way as in Kosmopoulos et al. [27]. The quantities presented in Figure 12 are averages of the results for the five locations and are indicative of the FL due to the effects of different parameters over the island.
As the direct component of the solar irradiance became more significant for the production of energy, the FL due to aerosols became larger relative to the FL due to clouds. Thus, for the CSP, the 1-axis tracking system, and the 2-axis tracking system installations, the annual FL due to aerosols were 18,085 euro (EUR), 24,199 EUR, and 31,741 EUR respectively, while the corresponding FL due to clouds were 8444 EUR, 15,437 EUR, and 17,327 EUR. This pattern is directly related to the results of Figure 7 indicating the major role of aerosol as attenuator during summertime, where the higher energy potential was followed by the higher aerosol loads and the subsequent effect on the EL and FL. For CSP and 2-axis tracking system installations, even the economic impact of dust was larger than the impact of clouds (FL of 9464 EUR and 17,373 EUR respectively). For horizontal PV installations the economic impact of clouds was more than double of the economic impact of aerosols (15,936 EUR relative to 6246 EUR). For tilted PV the overall FL due to aerosols and clouds were similar (11,287 EUR and 12,675 EUR respectively). The FL solely due to dust in this last case was 5476 EUR. The annual revenue is, as expected, maximum for the 2-axis solar tracking PV system (103,409 EUR) relative to other types of installations. It is worth mentioning that the financial gain of setting PV panels to a tilt equal to the latitude of the site was only 3% relative to having then at horizontal position (77,178 EUR instead of 74,847 EUR), mainly because of the very strong effect of aerosols on the direct solar beam. Due to the effects of aerosols and clouds the tilt angle that ensures optimal performance of the PV is not usually equal with the latitude of the installation location [103,104]. According to Jacobson and Jadhav [104] the optimal tilt for Larnaca is 30 • instead of 34.9 • (latitude of Larnaca). Analysis for Limassol using the methodology proposed by Raptis et al. [103] also resulted in an optimal tilt angle of 30 • instead of 35.2 • , which is the latitude of Limassol. As shown in Figure 13 changing the tilt angle from 35.2 to 30 • resulted in a gain in daily average energy of only 8 kWh/m 2 (gain of 10 MJ/m 2 in annual energy). The annual financial benefit from this energy gain (for a 500 kW system) would be 300 EUR.

Summary and Conclusions
In this study fine resolution satellite AOD and DOD products from the newly developed MIDAS climatology [49] were used for the quantification of the attenuation of SSR by aerosols and dust in Cyprus. CMF from the CM SAF-SARAH2.1 [105] product was used for the quantification of the effect of clouds. The spatial and temporal variability of the effects of aerosols and clouds on GHI and DNI, and the corresponding economic consequences on the production of solar energy were also investigated. Special emphasis was given to the effects of dust aerosols, which were found to constitute a large fraction of the aerosol mixture in Cyprus. Using the aforementioned information for clouds and aerosols climatological maps of the GHI and DNI were created.
In winter, clouds reduce GHI by 20-30% and DNI by 30-50% constituting the main attenuator of SSR in Cyprus. In spring and autumn, the role of clouds became less important, and in summer they blocked less than 10% of GHI and DNI. The effects of clouds were generally stronger over higher altitudes. The cycle of attenuation by aerosols was not as pronounced as the corresponding cycle of attenuation by clouds. Monthly average attenuation of GHI by aerosols ranged from 5% to 10%, with maximum values in spring and summer. Attenuation of GHI by aerosols was generally stronger, on average by 1-2%, at the south side of Troodos Mountains relative to the north. Depending on location and time of the year, the monthly average attenuation of DNI by aerosols ranged from 15% to 35% and was generally stronger in spring. At the south side of Troodos mountains aerosols blocked 10% more DNI relative to the north side. In summer, when SSR was maximum, aerosols constituted the main attenuator of GHI and DNI.
Dust aerosols constituted a large fraction of the aerosol mixture in Cyprus. In spring their average contribution to the overall aerosol load could be of the order of 60%. Furthermore, the short-and the long-term variability of aerosols depended strongly on the corresponding variability of dust. Thus, they had a strong impact on the production of solar energy over the island. On an annual basis, dust reduced GHI by 3% and DNI by 15%, and had a significant impact on energy production and the profit from PV and CSP installations. Indicatively, for a 500 kW installation of horizontal PV panels the average FL, solely due to dust, was 2800 EUR while for 2-axis solar tracking PV it was 17,000 EUR. For the same installations, the corresponding FL was 15,000 and 17,000 EUR due to clouds, and 6000 EUR and 32,000 EUR due to total aerosols. In reality, dust impacts energy production even stronger, since mechanisms such as the deposition of dust on the panels [106,107] and increase of atmospheric temperature due to the emission of infrared radiation by dust particles (and subsequently heating of PV) [108] were not taken into account in the above analysis.
Despite the relatively small area of the island of Cyprus, the complex terrain affected meteoclimatic conditions resulting to non-negligible spatial variability of GHI and DNI. Maximum GHI and DNI were estimated for the north of Troodos Mountains. Spatial variability in GHI was estimated to 10%. Spatial variability of DNI was larger, of 15%. The estimated average levels of GHI are in reasonable agreement with measured GHI reported in other studies for specific locations [16,17]. Comparison of the estimated GHI and DNI with corresponding estimates from CM SAF, resulted to a good agreement for lower altitude stations and larger differences (up to 20%) for higher altitudes. These differences can be attributed to differences in the AOD used for the simulations and to the way that CM SAF-SARAH2.1 takes into account changes in elevation [57]. Further validation is however out of the scope of the present study.
It is obvious that future changes in the climatological levels and physical and optical properties of aerosols and clouds over the area would significantly affect the production of solar energy. In 2004-2017 we found that the GHI increased by 1.8%, mainly because of decreasing attenuation by clouds. Projections of regional and climate models are however still uncertain regarding future changes in aerosols and clouds, and thus it is very difficult to estimate the magnitude and the direction of these changes [109,110]. Continuous monitoring of SSR and aerosol levels and properties in Cyprus are thus necessary for the accurate detection of long-term changes in these parameters, and subsequently future PV and CSP installation planning in an optimal way. Synergies of satellite and ground-based sensors are necessary in order to achieve high quality measurements with wide spatial coverage (e.g., [111,112]). Satellite retrievals are currently available on a global scale, with high spatial resolution, and as shown in Section 2 with a relatively good accuracy. Ground-based measurements are however necessary for the validation and the further improvement of satellite products.
The results of Section 3.4 also reveal the usefulness of SSR forecasts for Cyprus, since the day-to-day variability of solar irradiation can be large, even for cloudless days, due to the corresponding large variability of aerosol load. Accurate forecasts can provide very useful information for the optimization of solar energy production. For example, installations' maintenance can be planned for low-production days. Furthermore, energy providers can optimize energy management policies if they have accurate estimates of the total amount of energy entering the network in the coming hours or days. Again, complementary ground-based measurements are necessary for validation of the forecasts (e.g., [69]). Data Availability Statement: PollyXT lidar observations (level 0 data, measured signals) are in the PollyNET data base (https://www.tropos.de/en/research/projects-infrastructures-technology/ coordinated-observations-and-networks/pollynet (accessed on 11/06/2021)). AERONET observations were downloaded from the AERONET data base (https://aeronet.gsfc.nasa.gov/ (accessed on 11/06/2021)). The MIDAS data set is available at https://doi.org/10.5281/zenodo.4244106 (accessed on 11/06/2021) [49]. Climatologies of GHI and DNI are freely available at https://doi.org/10.5281/ zenodo.4737072 (accessed on 11/06/2021) [113].

Acknowledgments:
The ERATOSTHENES Centre of Excellence (www.eratosthenes.org.cy (accessed on 10 April 2021)), which was established after receiving funding by the Government of the Republic of Cyprus through the Directorate General for the European Programmes and the EXCELSIOR EU H2020 Widespread Teaming program with Grant Agreement No 857510 (www.excelsior2020.eu (accessed on 10 April 2021)). The MIDAS data set has been developed in the framework of the DUST-GLASS project (grant no. 749461; European Union's Horizon 2020 Research and Innovation programme under the Marie Skłodowska-Curie Actions). We thank AERONET for their continuous efforts in providing high-quality measurements and products and AERONET-Europe for providing an excellent calibration service. AERONET-Europe is part of the ACTRIS project under grant agreement no. 654109 and 739530 from the European Union's Horizon 2020 research and innovation program. We also thank the LACROS and PollyNET team for their well-organized website and data visualization. The authors gratefully acknowledge the NOAA Air Resources Laboratory (ARL) for the provision of the HYSPLIT transport and dispersion model and for the provision of Global Data Assimilation System (GDAS) data used in this publication. S.K. and K.P. would like to acknowledge the European Commission project EuroGEO e-shape (grant agreement No 820852). Finally, we thank the two anonymous reviewers and the editor for their constructive comments.

Conflicts of Interest:
The authors declare no conflict of interest.