Effect of Aerosols, Tropospheric NO2 and Clouds on Surface Solar Radiation over the Eastern Mediterranean (Greece)

In this work, the effect that two basic air quality indexes, aerosols and tropospheric NO2, exert on surface solar radiation (SSR) is studied, along with the effect of liquid and ice clouds over 16 locations in Greece, in the heart of the Eastern Mediterranean. State-of-the-art satellite-based observations and climatological data for the 15-year period 2005–2019, and a radiative transfer system based on a modified version of the Santa Barbara DISORT Atmospheric Radiative Transfer (SBDART) model are used. Our SSR simulations are in good agreement with ground observations and two satellite products. It is shown that liquid clouds dominate, with an annual radiative effect (RE) of −36 W/m2, with ice clouds (−19 W/m2) and aerosols (−13 W/m2) following. The radiative effect of tropospheric NO2 is smaller by two orders of magnitude (−0.074 W/m2). Under clear skies, REaer is about 3–4 times larger than for liquid and ice cloud-covered skies, while RENO2 doubles. The radiative effect of all the parameters exhibits a distinct seasonal cycle. An increase in SSR is observed for the period 2005–2019 (positive trends ranging from 0.01 to 0.52 W/m2/year), which is mostly related to a decrease in the aerosol optical depth and the liquid cloud fraction.


Introduction
Surface solar radiation (SSR) affects many climatic elements on a planetary scale, such as evapotranspiration [1], the hydrological cycle [2,3], photosynthesis and plant productivity [4,5], and the global (land and ocean) energy budget [6,7], largely affecting temperature and precipitation [8]. Apart from climate, SSR also affects solar energy production [9], which is vital to the economy.
The importance of solar radiation has led to the monitoring of SSR levels since the first decades of the 20th century, the first station being established in 1923 in Sweden. In the decades following, many more stations were established in Europe and around the globe, shaping global networks such as the Global Energy Balance Archive (GEBA) [10], the Baseline Surface Radiation Network (BSRN) [11], the Word Radiation Data Center (WRDC), and regional/national networks. Obviously, these networks cannot offer coverage over every single spot on Earth, a gap filled by satellite measurements during the last four decades (since the early 1980s) [12].
Various satellite missions played a vital role in improving our knowledge on the Earth's energy balance. Down-and up-welling solar radiation at the top of the atmosphere (TOA) is measured directly by satellite sensors, while surface solar radiation fluxes are estimated indirectly, involving the synergistic use of radiative transfer models with satelliteand ground-based measurements of various physical parameters [13][14][15][16][17][18][19].
Among other regions, various aspects of the SSR have been studied so far, by different groups over the Mediterranean Basin, which is a climate change hot spot [20]. In the last decade, the role of aerosols, clouds, and water vapor on SSR levels has been studied, mostly at a coarse spatial resolution, for all-and clear-sky conditions, based on satellite A core product of this research is the Moderate Resolution Imaging Spectroradiometer onboard Aqua satellite (MODIS/Aqua) collection 6.1 (C6.1) level-2 joint atmosphere product (MYDATML2). More specifically, single-pixel AOD 550 , cloud optical, microphysical and macrophysical properties, and water vapor column data retrieved from the MODIS/Aqua satellite sensor [33] were acquired for the period 2005-2019. Aqua, part of the A-train satellite constellation, has been flying at a sun-synchronous polar orbit at 705 km above the Earth since May 2002. MODIS/Aura is characterized by a wide spectral range (0.41 µm to 14.5 µm in 36 channels or bands), a broad swath (2330 km) and a spatial resolution of 250, 500 and 1000 m (depending on band), providing almost daily global coverage. MODIS/Aqua crosses Greece, which is our focus area here, at around 11:30 UTC.
For clouds, the following cloud optical, microphysical, and macrophysical properties from MODIS/Aqua C6.1 [45] were acquired: cloud fraction (Cloud_Fraction; CF), cloud phase (Cloud_Phase_Optical_Properties; CP), cloud optical thickness (Cloud_Optical_thickness; COT), cloud effective radius (Cloud_Effective_Radius; CER) and cloud top pressure (Cloud_Top_Pressure; CTP). CF is determined from the 1 km MODIS pixels, which are identified as cloudy and probably cloudy from the cloud mask product [46][47][48][49], and is given at a spatial resolution of 5 km × 5 km along with the other cloud properties. The C6.1 CP used here is retrieved by the daytime-only algorithm that uses a combination of visible (VIS), shortwave infrared (SWIR), and infrared (IR) channels (Marchant et al. [50] and references therein for details). The cloud optical and microphysical product algorithm makes use of six VIS, near-IR (NIR), SWIR, and mid-wave IR (MWIR) MODIS channels, as well as several thermal IR channels [45]. COT and CER are retrieved simultaneously based on the bispectral solar reflectance method [51]. CTP is retrieved by the cloud top properties algorithm, which relies on CO 2 -slicing channels (13-14 µm spectral region) and two IR window channels [52,53].
In addition to aerosol-and cloud-related parameters, total column precipitable water (Precipitable_Water_Near_Infrared_ClearSky; PW) data from the MODIS NIR algorithm [47] were acquired. The algorithm makes use of the ratios of water vapor absorbing bands (within the 0.94 µm water vapor band) with atmospheric window bands at 0.86 and 1.24 µm. The data are also available at a 5 km × 5 km spatial resolution.

CALIOP/CALIPSO Aerosol Profiles and MACv2 Aerosol Climatology
To get the full image of the aerosol loading and optical properties in the area, aerosol extinction coefficient profiles at 532 nm (Extinction_Coefficient_532_Mean) from the Cloud-Aerosol Lidar with Orthogonal Polarization onboard the Cloud-Aerosol Lidar and Infrared Pathfinder Satellite Observations satellite (CALIOP/CALIPSO) version 4.20 (v4.20) monthly level-3 product were acquired along with AOD, single scattering albedo (SSA) and asymmetry parameter (ASY) data at 8 wavelengths (350, 450, 550, 650, 1000, 1600, 2200 and 3000 nm) from the Max Planck Aerosol Climatology version 2 (MACv2) [54]. Similarly to Aqua, CALIPSO, also part of the A-train satellite formation (until September 2018), has been flying at a sun-synchronous polar orbit since April 2006. It crosses the equator at around 13:30 LT shortly after Aqua and has a 16-day repeat cycle [55]. CALIOP onboard CALIPSO is an Nd:YAG elastic backscatter and polarization lidar transmitting a linearly polarized laser light at 532 nm and 1064 nm [56]. CALIOP provides profiles of aerosol-and cloud-related properties from the ground up to 30 km [55]. The daytime (temporally coincident with MODIS/Aqua) CALIOP/CALIPSO data used here cover the period 6/2006-12/2019 being available at a 2 • (latitude) × 5 • (longitude) horizontal resolution and a 60 m vertical resolution (up to 12 km). The level-3 data are compiled from level-2 CALIOP/CALIPSO extinction coefficient data after applying various quality filters [57,58]. The MACv2 climatology data are available at a 1 o × 1 o spatial resolution and refer to the year 2005. MACv2 is compiled by fusing observational data from the Aerosol Robotic Network (AERONET) and the Marine Aerosol Network (MAN) with simulations from 14 different AeroCom (phase1) models [54].

OMI/Aura Tropospheric and Stratospheric NO 2
Another core product of this research is the Ozone Monitoring Instrument onboard Aura satellite (OMI/Aura) QA4ECV tropospheric (tropospheric_no2_vertical_column) and stratospheric (stratospheric_no2_vertical_column) NO 2 vertical column density (VCD) version 1.1 (v1.1) dataset [59]. OMI is a nadir-looking, pushbroom ultraviolet/visible (UV/VIS) solar backscatter grating spectrometer that has been flying onboard Aura satellite at a sun-synchronous polar orbit since July 2004, measuring the Earth radiance spectrum from 270 to 500 nm at a spectral resolution of around 0.5 nm [60]. Aura is part of the A-train satellite constellation with an equator overpass time 8 min (15 min before May 2008) after that of Aqua (13:38 LT). With its broad swath (2600 km) OMI provides almost daily global coverage at a spatial resolution of 13 km × 24 km (at nadir).
The QA4ECV algorithm, which is described in detail in Boersma et al. [59], comprises the following three basic steps: 1) the total NO 2 slant column density (SCD) is calculated using the differential optical absorption spectroscopy (DOAS) method, 2) the stratospheric contribution to the SCD is calculated from the TM5-MP chemistry transport model (CTM) with data assimilation and it is subtracted from the total SCD to obtain the tropospheric SCD, and 3) the tropospheric/stratospheric SCDs are converted to tropospheric/stratospheric VCDs by dividing by the corresponding air mass factors (AMFs).
The data used here are gridded at a 0.125 • × 0.125 • spatial resolution on a monthly basis covering the period 2005-2019. While averaging, each single-pixel observation is weighted by its fractional area within the grid cell. Pixels with cloud radiance fraction greater than 50% (cloud fraction greater than about 20%) and solar zenith angles higher than 80 • are filtered out. The monthly gridded tropospheric NO 2 data are available on the Tropospheric Emission Monitoring Internet Service (TEMIS; http://www.temis.nl, accessed on 30 June 2021).

CM SAF, CERES and Ground-Based SSR Data
For evaluation purposes, SSR satellite-based observations were acquired from the Satellite Application Facility on Climate Monitoring Surface Solar Radiation Data Set-Heliosat-edition 2.1 (CM SAF SARAH-2.1) product [65]. SARAH-2.1 is an extension of SARAH-2 [66], being the latest release of a series of CM SAF SSR products, which are based on observations from the Meteosat First and Second Generation (MFG/MSG) geostationary satellite series covering Africa, Europe, and the Atlantic Ocean (Mueller et al., [18] and references therein). The effective cloud albedo calculated with the Heliosat algorithm is used in combination with the clear-sky surface radiation model SPECMAGIC to derive SSR (in the 200-4000 nm wavelength region) from the Meteosat Visible-InfraRed Imager (MVIRI) and Spinning Enhanced Visible and Infrared Imager (SEVIRI) sensors onboard MFG and MSG satellites based on a modified Monitoring Atmospheric Composition and Climate (MACC) aerosol climatology [18]. The CM SAF SARAH-2.1 SSR data are available on a 30 min, daily and monthly basis at a spatial resolution of 0.05 • × 0.05 • spanning the period 1983-2017. In this work, monthly data are used for the period 2005-2017. Evaluation of the SARAH-2.1 monthly product over Europe [66] has shown good agreement with ground-based pyranometric observations (mean bias of +2.7 W/m 2 , mean absolute bias of 6.9 W/m 2 and a Pearson's correlation coefficient R of 0.89).
Another satellite-based SSR dataset used here for evaluation comes from the Clouds and the Earth's Radiant Energy System Energy Balanced and Filled edition 4.1 (CERES EBAF edition 4.1) level-3b product [67,68]. The data are available on a monthly basis at a 1 • × 1 • spatial resolution. The product is mostly based on measurements from the CERES sensors onboard Terra and Aqua satellites, MODIS measurements from the same satellites and geostationary imagers. CERES measures filtered radiances in the shortwave (SW; 0.3-5 µm), total (TOT; 0.3-200 µm) and window (wavelengths between 8 and 12 µm) spectral regions with a spatial resolution of~20 km at nadir [69]. Unfiltered SW, longwave (LW) and window radiances are calculated [70] and finally are converted to fluxes using empirical angular distribution models (ADMs) [71]. For the calculation of surface fluxes cloud data from MODIS and geostationary satellites are utilized along with aerosol properties from MATCH (Model of Atmospheric Transport and Chemistry) that assimilates MODIS AODs, atmospheric profile data from the Goddard Earth Observing System version 5.4.1 (GEOS-5.4.1) reanalysis, and ancillary atmospheric profile data from other satellite sensors (i.e., CALIPSO, Cloudsat and AIRS/Aura) [68]. The radiative transfer calculations are constrained by the CERES TOA fluxes. The uncertainty of the 1-degree CERES EBAF edition 4.1 SSR product used here is less than 5 W/m 2 , which is within the uncertainty of ground measurements [68].
Finally, for evaluation purposes, hourly pyranometric data from three ground stations with quality-assured long time series [27] are utilized here. The first station (urban) is located in the center of Thessaloniki within the university campus of the Aristotle University of Thessaloniki (40.63 • N, 22.96 • E). The measurements come from a CM21 pyranometer (Kipp and Zonen, Delft, the Netherlands) [72]. The instrument measures solar radiation in the spectral range 305-2800 nm with an uncertainty of ±2%. The second station (urban) is located in the center of Athens, Greece at the facilities of the National Observatory of Athens (37.97 • N, 23.72 • E). The measurements were taken from an Eppley precision spectral pyranometer that measures solar radiation in the wavelength range from 285 to 2800 nm with an uncertainty of ±1-2% [73]. The third station (rural) is located in a coastal area (35.33 • N, 25.28 • E) 15 km east from the city of Herakleion at the premises of the Hellenic Center for Marine Research (HCMR) on the island of Crete. The measurements come from a CM11 pyranometer that measures solar radiation in the spectral range 305-2800 nm with an uncertainty of ±3% (Alexandri et al., [27] and references therein). The ground-based observations from Thessaloniki and Athens used here cover the period 2005-2019 while the data from HCMR cover the period 2005-2013.

Ancillary Data
Along with the basic data products, in this work a number of ancillary datasets is utilized. OMI/Aura O 3 total column data from the 1 • × 1 • daily gridded dataset OMTO3d version 3 (v3) [74] were acquired for the period 2005-2019. The OMI science team produces this level-3 product by gridding and averaging only good-quality level-2 total column ozone orbital swath data retrieved by the TOMS version 8 algorithm. Also, satellite-based monthly gridded (0.25 • × 0.25 • resolution) surface albedo (SAL) data from the CM SAF Cloud, Albedo, Radiation AVHRR-based edition 2.1 (CLARA-A2.1) product [75] were acquired for the period 1/2005 to 6/2019. The surface albedo data are retrieved from the Advanced Very High Resolution Radiometer (AVHRR) instruments onboard the National Oceanic and Atmospheric Administration (NOAA) and Meteorological Operational Satellite Program of Europe (MetOp) satellites and cover the spectral range 0.25-2.5 µm. Another CM SAF product used for the scopes of this research is the CLoud property dAtAset using SEVIRIedition 2.1 (CLAAS-2.1)-cloud fractional cover (CFC) [76]. More specifically, daytime CFC data were acquired on a monthly basis at a spatial resolution of 0.05 • × 0.05 • and monthly diurnal data (for 11:00 and 12:00 UTC) which are available at a spatial resolution of 0.25 • × 0.25 • . The CFC measurements are retrieved by the SEVIRI instruments onboard the MSG geostationary satellite series. Finally, the surface topography data used here are from the Global Multi-resolution Terrain Elevation Data 2010 (GMTED2010) product [77] and the population data for the year 2020 from the Gridded Population of the World version 4 (GPWv4.11) dataset [78].

Data Aggregation
Within the framework of the project "Long-term variability of solar radiation in Greece: effect of air quality and clouds" all the datasets presented above have been spatially and temporally homogenized. The MODIS/Aqua aerosol, cloud, and water vapor data have been gridded on a monthly basis at a standard 0.125 • × 0.125 • grid covering the Eastern Mediterranean for the period 2005-2019. The spatial resolution of the grid was selected as the best compromise between the pixel size (at nadir) of the different core satellite observations used here (from MODIS/Aqua clouds 5 km × 5 km to OMI/Aura NO 2 13 km × 24 km). For clouds, taking into account the MODIS/Aqua cloud phase parameter their optical, microphysical and macrophysical properties are calculated separately for liquid (CP: 2 or 4) and ice clouds (CP: 3). While gridding, each aerosol, cloud and water vapor single-pixel observation is weighted by its fractional area within the grid cell following the gridding methodology of the Royal Netherlands Meteorological Institute (KNMI) for the OMI/Aura tropospheric and stratospheric NO 2 data [32]. The weighted pixel values are averaged on a daily basis for each grid cell and finally the daily gridded data are averaged on a monthly basis. All the other datasets were gridded (at the same standard 0.125 • × 0.125 • grid) on a monthly climatological basis (12 values per grid cell) using bilinear interpolation.

A System for Radiative Transfer Calculations-SBDART-NO 2
A system that implements radiative transfer calculations at any spatial and temporal resolution using gridded atmospheric properties data was developed (see Figure 1). It comprises an interactive data language (IDL) interface and an updated version of the wellestablished Santa Barbara DISORT Atmospheric Radiative Transfer (SBDART) model [79].
The interface gathers all the information from the gridded data and builds all the necessary input parameters for the radiative transfer model. Then, it executes SBDART for the following three sky conditions: clear, liquid and ice cloud-covered skies.
The system uses AOD 550 from MODIS/Aqua and climatological MACv2 AOD values at 8 different wavelengths (Section 2.1.2) to account for the spectral dependence of AOD. SSA and ASY come from MACv2. The MACv2 parameters characterize the aerosol mixture over each location. Climatological CALIOP/CALIPSO aerosol extinction profiles account for the shape of the vertical profile of aerosols. The CALIOP/CALIPSO profile is scaled accordingly in every simulation to make sure that the total AOD is equal to that of MODIS/Aqua. For clouds, MODIS/Aqua CF, COT, CER and CTP data for liquid and ice clouds separately are used. Cloud bottom pressure (CBP) is calculated from CTP for low (CTP > 680 mbar), medium (440 mbar < CTP < 680 mbar), and high clouds (CTP < 440 mbar) using the corresponding cloud depth (CD) values (100, 50 and 50 mbar) from Hatzianastassiou and Vardavas [80] (CBP = CTP + CD). CTP and CBP are converted to cloud top and bottom height and temperature using the McClatchey et al. [81] standard atmospheres. From the cloud top and bottom temperatures the mean temperature at the cloud layer (T c ) is derived. The relative humidity at the cloud layer, which is used by SBDART to define the water vapor profile, is calculated from T c following Korolev and Isaac [82]. MODIS/Aqua PW is used for water vapor column while broadband surface albedo is from CM SAF CLARA-A2.1.
SBDART is mostly oriented on aerosol-and cloud-related studies and the spectral absorption of NO 2 in the UV/Vis is missing from the original version. To account for the radiative effect of NO 2 , we updated SBDART by incorporating into the model the NO 2 cross sections from Burrows et al. [83] in the spectral range from 231 to 794 nm with a spectral resolution of 20 cm −1 . On top of that, the updated version of SBDART (hereafter also denoted as SBDART-NO 2 ) was further modified so as to account for tropospheric and stratospheric NO 2 , separately. The predefined SBDART tropospheric and stratospheric profiles are shifted accordingly (multiplied by a factor) so that the tropospheric and stratospheric NO 2 columns are equal to the satellite-based ones from OMI/Aura. The ability of SBDART-NO 2 to determine the radiative heating produced by NO 2 within the atmosphere has been validated. More specifically, the calculated radiative heating values produced by 0.5 DU of NO 2 distributed within the first 1 km of the atmosphere under different solar zenith angle and surface albedo scenarios were identical to the ones appearing in Vasilkov et al. [30] for the same specifications. The radiative transfer calculations are implemented with SBDART-NO 2 for 185 wavelengths in the range from 250 to 4000 nm with a decreasing resolution of 4-10 nm in the UV/Vis and 10-60 nm thereafter. The daily average solar zenith angle (SZA d ) is derived following Vardavas and Koutoulaki [84].
As the radiative transfer calculations here are implemented on a monthly basis, the monthly mean SZA d is used as input by SBDART. The results for clear, liquid and ice cloudcovered skies are then multiplied by the monthly mean daylength [84,85] and weighted by the corresponding clear-sky/cloud fractions to obtain the monthly mean all-sky daily SSR (e.g., Hatzianastassiou et al. [15]). Prior to weighting, a correction is applied to the MODIS/Aqua cloud fraction data as they characterize only a very specific time range in the afternoon and not the daytime cloud cover. The correction is obtained by multiplying by climatological monthly ratios of daytime cloud cover and cloud cover at 11-12:00 UTC (overpass time of Aqua) from CM SAF CLAAS-2.1. By comparing SSR simulations with and without this correction with ground-based observations we found that the correction reduces the bias by 3-4%. The radiative transfer system is capable of producing gridded solar radiation fields at the top of the atmosphere, within the atmosphere and at the surface similarly to other global and regional systems [15,17,19,21]. Following the study of Benas et al. [86], here we show a demonstration for SSR for the period 2005-2019 for 16 selected stations around Greece, located in the heart of the climatically sensitive region of the Eastern Mediterranean. The locations (see Table 1 for details) are divided into urban and rural taking into account population (GPWv4.11) and air quality criteria (MODIS/Aqua AOD 550 and OMI/Aura tropospheric NO 2 ; Figure 2), and cover the whole country. Also, the 16 locations are embracing stations where there are or there may be solar radiation measurements in the future, e.g., stations of the Hellenic Network of Solar Energy (HNSE) [87], the PANhellenic GEophysical observatory of Antikythera (PANGEA) [88] station and the Olympus Research Center (EKO) station.

Evaluation of SSR Data, Radiative Effects, and Trends
The SSR data produced by the radiative transfer system described above are evaluated against ground-based observations from three stations in Greece, and CM SAF SARAH-2.1 and CERES EBAF edition 4.1 satellite-based observations for the 16 selected locations in Greece. Then the radiative effect of aerosols, tropospheric NO 2 , and liquid and ice clouds is calculated by running the system five times, with the following conditions: (1) with all the parameters, (2) without aerosols, (3) without tropospheric NO 2 , (4) without liquid clouds and (5) without ice clouds. The radiative effect of each parameter i (RE i ) is defined as the difference between the SSR from the simulation with all the parameters and the SSR from the simulation without i (SSR witout-i ) as shown in Equation (1).
The trends of the radiative effects of the aerosols, tropospheric NO 2 and clouds, and the related parameters are examined for the period 2005-2019 using the well-established method of Weatherhead et al. [89] as described in detail in Alexandri et al. [27]. The monthly timeseries are fitted by using a model with a linear trend and a Fourier-based seasonal component for the annual cycle (Equation (2)).
Y t is the monthly mean of month t, X t is the period of time (in fractional years) after January 2005 (t/12), A is the monthly mean for January 2005 and B is the annual trend. The seasonal component contains the amplitudes a n and b n , T is the period (one year) and N t is the remainder. The remainder represents the difference between the modeled and the measured value and is given by the following formula: N t = ϕN t−1 + ε t . ϕ is the autocorrelation in the remainder and ε t the white noise. The precision of the trend σ B is given as a function of ϕ, the length m of the dataset in years, and the variance σ N of the remainder for small autocorrelations (Equation (3)).
The trend ω is considered statistically significant at the 95% confidence level if |ω/σ B | is greater than 2. Only trends for timeseries with at least 8 months per year are presented here following previous studies [90,91]. Specifically for the radiative effect, which takes negative values, the trends refer to its magnitude (absolute values of the radiative effect) rather than its algebraic values, for easier interpretation of the results (e.g., decreasing AODs over the years would lead to a negative trend of the radiative effect magnitude).

SSR Data Evaluation
To gain confidence in the monthly solar radiation patterns and the corresponding radiative effects calculated by our radiative transfer system, our SSR simulations (SBDART) are validated against monthly SSR observations from three ground stations with qualityassured long time series [27], and satellite-based observations from CM SAF SARAH-2.1 and CERES EBAF edition 4.1 for the 16 locations investigated here. As shown in Figure 3a, SSR is underestimated in our simulations by −7.5% (−14.5 W/m 2 ) compared to the groundbased observations, with a very good correlation (regression line slope and Pearson's correlation coefficient R of~1). The high R values mentioned here show, to a great extent, that our SSR simulations capture the seasonal cycle of SSR well [92]. To gain a better insight into the quality of our simulations, relative to the ground-based observations, we also compare the deseasonalized values (Figure 3b). The correlation is again very good, with R being 0.78, similarly to other studies comparing deseasonalized SSR data [92,93]. Apart from the normalized mean bias (NMB) and the mean bias (MB) referenced here, other metrics, i.e., regression line, mean absolute bias (MAB), and root mean square error (RMSE), are given in Figure 3. A negative bias appears in the comparison against the CM SAF and CERES data (−9.1% and −7.2%) as well, again with a very good correlation (Figure 3c-f) for both the original and the deseasonalized timeseries. The ability of our system to simulate the monthly SSR values is similar to that of other systems (e.g., Papadimas et al. [21]), which are based on a monthly mean daily SZA and daylength. Papadimas et al. [21] simulations underestimated the SSR by −5.3% compared to the monthly observations from the Global Energy Balance Archive (GEBA) for the Mediterranean Basin, with a regression line slope of 0.95 and R of 0.95. Specifically for the station of Thessaloniki, our SSR simulations agree with the ground-based observations better than theirs (NMB of −6.5% vs. −10.2%) ( Figure S1 in the electronic supplement of this paper). An NMB of −10.2% compared to the daily instantaneous ground-based SSR observations from HCMR was shown in another study that focused on the radiative effect of aerosols [86]. These results were obtained using a different version of the radiative transfer model used by Papadimas et al. [21], with aerosol input data from MODIS/Aqua. Our SSR simulations over HCMR exhibit better agreement with the ground-based observations (NMB of −6% with a very good correlation) (see Figure S1). However, it has to be noted that our simulations are monthly and do not refer to instantaneous SSR values.
The evaluation of our SSR simulations against ground-based observations, and against the CM SAF SARAH-2.1 and CERES EBAF edition 4.1 satellite-based products for each of the 16 locations of interest, is shown in Figures S2, S4 and S5, respectively. The evaluation of the CM SAF and CERES data against the ground-based observations is shown in Figure  S3. From all those figures, we conclude that our monthly SSR simulations are in good agreement with the ground-based observations and other satellite-based products, which gives confidence in the radiative effect calculations presented in the next paragraphs.

Annual
The annual all-sky radiative effect of aerosols, tropospheric NO 2 , and liquid and ice clouds over the 16 locations examined here is shown in Figure 4. It is obvious that liquid clouds dominate, with an average radiative effect (RE lc ) of −36 W/m 2 . Ice clouds exhibit a smaller radiative effect (RE ic ) of −19 W/m 2 , with the radiative effect of aerosols (RE aer ) following (−13 W/m 2 ). It has to be highlighted that the radiative effect of tropospheric NO 2 (RE NO2 ) is smaller by two orders of magnitude (−0.074 W/m 2 ) compared to that of aerosols and clouds and hence is negligible.
Aerosols exhibit a stronger radiative effect, by~2 W/m 2 , over the urban locations (14 W/m 2 ) compared to the rural locations (12 W/m 2 ). The annual RE aer ranges from −18 W/m 2 (Patras) to −8 W/m 2 (Olympus). It has to be highlighted here that even over the Olympus location, a high-altitude background area entailing the top of Mount Olympus (2917 m), RE aer is still −8 W/m 2 , which is very close to other rural (Xanthi and Ioannina; −10 W/m 2 ) and even urban locations (Kozani; −10 W/m 2 ). The annual RE aer values for all the 16 locations of interest can be found in Figure 4a and in Table 2. The values presented in this work are in accordance to previous studies in the area. Papadimas et al. [21] reported an average RE aer value of −17.4 W/m 2 over the Eastern Mediterranean, in the electronic supplement of their paper. For the same area, and based on climatological calculations, Alexandri et al. [27] suggested that aerosols reduce incoming shortwave solar radiation bỹ 18 W/m 2 . However, both these studies refer to later periods. Hence, taking into account the decreasing AOD trends in the area (see Section 3.3.2), the smaller RE aer values we obtain here make sense.
The radiative effect of tropospheric NO 2 over urban locations (−0.113 W/m 2 ) is double that over rural locations. The highest (still negligible) RE NO2 values appear over the following three pollution hot-spot locations in the area (see Figure 2d (Figure 4b). Athens is the capital city, and industrial and commercial center of Greece, with a population of approximately four million people in its metropolitan area. Thessaloniki is the second largest city, with a population of approximately one million people in its metropolitan area, while the greater area of Kozani is a power-producing center entailing lignite-burning power plants. However, the tropospheric NO 2 levels in those areas are still quite low (~4 × 10 15 molecules/cm 2 ) to exert a significant radiative effect on an annual basis compared to the heavily polluted areas of Eastern China where tropospheric NO 2 levels can be up to~10 times higher [32].  Table 2). The annual RE lc and RE ic for all the 16 locations of interest can be seen in Figure 4c,d, respectively. In general, RE lc and RE ic exhibit a latitudinal south-to-north increasing gradient in accordance to previous results in the area [22,94]. Pyrina et al. [94], using a radiative transfer model and satellite-based aerosol and cloud observations at a coarse spatial resolution, calculated an average net surface shortwave cloud radiative effect of −44.5 W/m 2 over the Eastern Mediterranean. This is close to the combined (liquid and ice) cloud radiative effect (−55 W/m 2 ) we calculated here, taking into account that in their analysis they included oceanic areas where cloud coverage is quite limited compared to the land areas, and hence the cloud radiative effect is smaller ( Figure 5). It is obvious that clouds dominate, and hence RE aer and RE NO2 are expected to be different under clear and cloudy skies. In Figure 6, the radiative effect of aerosols and tropospheric NO 2 over the 16 locations studied here is shown for clear, liquid cloudcovered, and ice cloud-covered skies. We see that the average RE aer for all the locations under clear skies (−18 W/m 2 ) is about 3-4 times larger than for liquid (−5 W/m 2 ) and ice (−6 W/m 2 ) cloud-covered skies; for individual locations, it can be up to 4-5 times larger (e.g., Olympus). For tropospheric NO 2 , the average RE NO2 for all the locations under clear skies (−0.095 W/m 2 ) is double that for liquid (−0.045 W/m 2 ) and ice (−0.051 W/m 2 ) cloud-covered skies. RE aer and RE NO2 are smaller under cloudy conditions, as part of the solar radiation is already reflected by clouds before reaching the aerosol and tropospheric NO 2 layers, which, on a climatological basis, are primarily residing below the clouds.

Seasonal
The magnitude of the radiative effect of aerosols, tropospheric NO 2 , and liquid and ice clouds exhibits a seasonal variability that is similar to that of AOD 550 , tropospheric NO 2 VCD, liquid cloud fraction and ice cloud fraction, as shown in Figures S5-S8, for each one of the 16 locations. This is clearer in the case of aerosols and clouds, and less distinct for tropospheric NO 2 where RE NO2 is nearly negligible. Specifically for Athens and Thessaloniki, the two urban centers of Greece, the seasonal variabilities of all the aforementioned parameters are shown in a single comparative figure (Figure 7).
The seasonal variability of RE aer and AOD 550 for each one of the stations is shown in Figure S5. It is well known that AOD 550 in the area peaks in the spring and summer, as a result of the transport of dust from the Sahara Desert and the Middle East, pollution aerosols from Central Europe, and biomass burning aerosols from Southwestern Russia and Eastern Europe [26]. Other factors, such as the low precipitation amounts and the photochemical production of secondary organic aerosols in the summer, also play a substantial role. RE aer peaks in the spring and summer, with maximum values up to −30 W/m 2 (e.g., Rhodes, Preveza, Thessaloniki) appearing in August over all the 16 locations. During the winter, especially in December and January, RE aer generally drops below −5 W/m 2 and in some cases reaches values around −1 W/m 2 (e.g., Olympus, Kozani, Ioannina). In Figure 7a,b, we see that over Athens, AOD 550 exhibits higher values in the spring compared to Thessaloniki and other locations in the northern part of Greece, where AOD 550 peaks in the summer. This is related to the fact that Northern Greece is closer to the biomass burning hot spots of the area (biomass burning season in July-August), while Southern Greece is closer to the Sahara Desert (dust transport peaks in the spring). In Figure 7c,d, distinct winter/summer tropospheric NO 2 maxima/minima are observed over Athens and Thessaloniki, with VCDs ranging 2-3 to 5-6 × 10 15 molecules/cm 2 . This is most probably related to the increased consumption of fossil fuel and biomass burning for heating purposes in the winter, and the decreased traffic in the summer. On the contrary, over island locations where tourism peaks in the summer, we observe either an opposite seasonal variability with a clear summer maximum (e.g., Rhodes) or a secondary summer peak (e.g., Mytilene) ( Figure S6). The radiative effect of tropospheric NO 2 is overwhelmed by the effect of other parameters, especially aerosols, water vapor and clouds, peaking in July over all the locations examined here (Figure 7c,d and Figure S6). The highest RE NO2 value appears over Athens in July (−0.215 W/m 2 ), which is still very low to play a substantial role in the local radiative budget.  The liquid cloud fraction in the area exhibits a winter peak (values up to 0.5) with a summer minimum, especially in July-August ( Figure S7), being almost zero over specific locations, mostly in Southern Greece (e.g., Rhodes). RE lc exhibits a quite similar seasonal variability with liquid cloud fraction, with clear summer minima over all the locations. The lowest RE lc levels appear in July over the island locations of Rhodes (−1 W/m 2 ) and Antikythera (−3 W/m 2 ). The highest RE lc appears in May over the mountainous location of Olympus (−123 W/m 2 ). Over Olympus, RE lc takes values higher than −40 W/m 2 throughout the year. As shown in Figure 7e,f, over Athens and Thessaloniki, RE lc peaks in the winter-spring (values up to −40 to −50 W/m 2 ), with a clear summer minimum.
The seasonal variability of the ice cloud fraction is similar to the liquid cloud fraction, with a winter peak (values up to 0.3) and a summer minimum when the presence of ice clouds is very limited, over all the 16 locations ( Figure S8). In general, RE ic exhibits quite a similar seasonal variability with ice cloud fraction, with clear summer minima over all the locations and peaks in the spring. In July, RE ic is zero over the island/coastal locations of Rhodes/HCMR, while for others it is lower than −2 W/m 2 in July and/or August  Figure S8).

SSR Trends
To gain insight into the SSR changes during the fifteen-year period (2005-2019), trends are calculated for each one of the 16 locations, as described in detail in Section 2.2.3. The trends of the SSR values produced by our system (SBDART), from the CM SAF SARAH-2.1 (2005-2017) and the CERES EBAF edition 4.1 satellite-based products, and from the two ground stations (Athens and Thessaloniki) covering the whole time period, appear in Figure 8. The SSR timeseries, along with the trends and relative trends (trends relative to the first year, as in Alexandri et al. [27]), and the corresponding uncertainties (±1σ) for each station, are shown in Figures S9-S11, correspondingly.
It is shown, in Figure 8a and Table 3, that positive SBDART SSR trends ranging from 0.01 (Rhodes) to 0.52 W/m 2 /year (Olympus) appear over all the locations except for two, which are Mytilene (−0.27 W/m 2 /year or −0.14%/year) and HCMR (−0.12 W/m 2 /year or −0.06%/year). Only over the two locations with the strongest trends, Volos (0.78 W/m 2 /year or 0.45%/year) and Thessaloniki (0.67 W/m 2 /year or 0.40%/year), are the trends statistically significant at the 95% confidence level. The SBDART SSR trends, in general, agree with those of CM SAF ( Figure 8b) and CERES (Figure 8c). Differences in the magnitude of the trends are expected, due to the special characteristics of the satellite-based products (sensors with different spatial and temporal resolution), and the different methods and assumptions involved in the calculations.
Our SBDART SSR trends are in better agreement with the CERES trends than the CM SAF. Apart from the slightly different period (2005-2017) covered by the CM SAF data, it is also well known that this product is based on climatological aerosol data [18], and hence misses the trends induced by aerosols compared to our SBDART SSR and CERES (both use MODIS aerosol data). On the other hand, due to its high spatial (0.05 • × 0.05 • ) and temporal (30 min) resolution, this product probably obtains a better representation of cloud radiative effects [27]. The CERES EBAF edition 4.1 product also incorporates cloud information from geostationary satellites, and the radiative calculations are resolved on an hourly basis. However, its low spatial resolution (1 • × 1 • ) will definitely impact the trend calculations.
Compared to the ground-based observations from Athens and Thessaloniki, it is shown, in Figure 8d and Figure S12, that the strong positive SSR trend over Athens (0.66 W/m 2 /year or 0.34%/year) is captured better by our SBDART SSR than CM SAF and CERES, but, contrary to the ground-based trends, statistical significance is not indicated. Over Thessaloniki, ground-based observations exhibit a statistically insignificant positive trend (0.26 W/m 2 /year or 0.14%/year), which is better captured by CERES compared to SBDART and CM SAF. The evaluation of the SSR trends against the data from two ground stations only does not allow for safe conclusions; however, it is pretty enhancing that all the three independent datasets capture the sign of the trends.

Trends of the Aerosol, Tropospheric NO 2 and Clouds Radiative Effect
The contribution of the aerosol, tropospheric NO 2 , and cloud radiative effect on the SSR trends is reflected in the RE aer , RE NO2 , RE lc and RE ic magnitude trends, which are shown in Figure 9 and in Table 3. The corresponding timeseries, along with the trends, the relative trends, and their uncertainties for each station, are shown in Figures S13-S16, correspondingly. As the radiative effect takes negative values, it is underlined once more here (see also Section 2.2.3) that the trends refer to the magnitude of the radiative effect rather than its algebraic values. Table 3. Trend (in W/m 2 /year) and relative to the first-year trend (in parentheses; in %/year) of the SSR (data produced in this work based on SBDART) and the radiative effect of aerosols, tropospheric NO 2 , liquid clouds and ice clouds for the period 2005-2019 with the corresponding ±1σ. Red/blue color denotes positive/negative trends, respectively. Bold underlined characters are used for statistically significant trends at the 95% confidence level. Negative RE aer magnitude trends are found over all the locations, except for HCMR. The strongest statistically significant RE aer magnitude trends to appear over the urban locations of Patras (−0.35 W/m 2 /year or −1.74 %/year) and Thessaloniki (−0.31 W/m 2 /year or −2.05 %/year), with Xanthi, Argos, Volos, Orestiada and Antikythera following (see Table 3). The negative RE aer magnitude trends are driven by the negative AOD 550 trends shown in Figure 10a and Figure S17, and in Table 4. The strongest statistically significant AOD 550 trends are seen over Patras, Thessaloniki, and Xanthi (−0.005/year or −2 to −2.5%/year). The AOD 550 trend values for all the 16 locations of interest shown in Table 4 (ranging from −0.005/year to zero) are in accordance to the average trend value of −0.0022/year for the Eastern Mediterranean, given in a previous study using C6 MODIS/Aqua data for the period 2002-2015 [90].

Location
The RE NO2 magnitude, despite being negligible as shown in Section 3.2, exhibits a negative trend over all the locations examined here, except for one (Antikythera) (see Figure 9b and Figure S14, and Table 3). The strongest statistically significant trends appear over the urban locations of Kozani (−0.008 W/m 2 /year or −3.91%/year), Athens (−0.008 W/m 2 /year or −3.44%/year) and Thessaloniki (−0.004 W/m 2 /year or −2.47%/year). The negative RE NO2 magnitude trends reflect the negative trends of the tropospheric NO 2 VCD appearing over all the locations, except for Antikythera (Figure 10b and Figure S18, and Table 4). The strongest statistically significant trends are seen over Kozani (−0.21 × 10 15 molecules/cm 2 /year or −3.95%/year), Athens (−0.19 × 10 15 molecules/cm 2 /year or −3.57%/year) and Thessaloniki (−0.11 × 10 15 molecules/cm 2 /year or −2.89%/year). However, it is obvious that tropospheric NO 2 cannot affect the SSR trends over the area examined here. The RE lc magnitude exhibits negative trends over all the 16 locations. However, the trends are statistically significant over the urban location of Volos only where the strongest trend is found (−0.80 W/m 2 /year or −1.78%/year), with Olympus and Thessaloniki following (Figure 9c and Figure S15, and Table 3). The negative RE lc magnitude trends are not statistically significant except for one location, despite the fact that statistically significant liquid cloud fraction trends are observed over nine locations ( Figure S19 and Table 4). The strongest liquid cloud fraction trends are seen over Volos (−0.006/year or −1.88%/year), with Xanthi and Olympus following (−0.004/year or −0.8 to 1.2%/year). The effect of the negative liquid cloud fraction trends on RE lc is partly counterbalanced by the positive liquid cloud optical thickness trends over some locations (Ioannina, Kozani, Mytilene, Olympus, Orestiada, Pylos, Thessaloniki, and Volos) (see Figure 10c,e and Figure S21, and Table 4). The liquid cloud optical thickness trends are statistically insignificant, ranging from -0.236/year or −1.24%/year (Thessaloniki) to 0.188/year or 1.03%/year (Xanthi) ( Table 4).  Table 4. Trends of AOD 550 (in 1/year), tropospheric NO 2 VCD (in 10 15 molecules/cm 2 /year), liquid cloud fraction (CFL; in 1/year), liquid cloud optical thickness (COTL; in 1/year), ice cloud fraction (CFI; in 1/year) and ice cloud optical thickness (COTI; in 1/year) for the period 2005-2019 with the corresponding relative to the first-year trends (in parentheses; in %/year) and ±1σ. Red/blue color denotes positive/negative trends, respectively. Bold underlined characters are used for statistically significant trends at the 95% confidence level. Aerosol-and cloud-related parameters are from MODIS/Aqua and tropospheric NO 2 from OMI/Aura. Contrary to liquid clouds, the RE ic magnitude exhibits positive trends over all the 16 locations except for one (Rhodes). The trends are statistically significant over the rural location of Argos only. Argos exhibits the strongest trend (0.46 W/m 2 /year or 3.57%/year), with Pylos and Olympus following (Figure 9d and Figure S16, and Table 3). The positive RE ic magnitude trends follow the positive ice cloud fraction trends. In line with RE ic , statistical significance is indicated only over Argos, which is the location with the strongest trend (0.003/year or 2.13%/year) along with Volos (0.003/year or 1.70%/year) ( Figure S20 and Table 4). Similarly to the case of liquid clouds, the effect of the positive ice cloud fraction trends on RE ic is partly counterbalanced by the negative ice cloud optical thickness trends over some locations (Kozani, Mytilene, Olympus, Orestiada, Patras, Preveza) (see Figure 10d,f and Figure S22, and Table 4). The ice cloud optical thickness trends are statistically insignificant except for Pylos, ranging from −0.218/year or −1.01%/year (Olympus) to 0.615/year or 5.53%/year (Pylos) ( Table 4).

Location
Concluding, among the parameters that are studied here, AOD 550 and liquid cloud fraction drive the positive SSR trends observed in the area. These parameters generally decreased over the area within the period 2005-2019, while, on the contrary, the ice cloud fraction increased, exhibiting an opposite radiative effect change. Also, for some stations, the liquid and ice cloud optical thickness has increased, and for some others it has decreased during the same period, suppressing or boosting the changes in the radiative effect of clouds induced by changes in the liquid and ice cloud fractions. Our results support previous regional and global studies that underlined the importance of cloud optical thickness in the determination of SSR trends and the diming/brightening effect (e.g., Alexandri et al. [95]; Kambezidis et al. [96]; Hatzianastassiou et al. [97]).

Discussion
Here, various aspects related to the quality of the simulations and the interpretation of our results are discussed. This discussion is relevant to numerous other studies using similar radiative transfer systems. In Section 3.1, it is shown that our monthly SSR simulations are in good agreement with ground-and satellite-based observations exhibiting a negative bias. Cronin [98] stressed that the use of a daytime-average zenith angle may lead to high biases in the amount of radiation reflected by the clouds. Indeed, our clear-sky SSR values are very close to the corresponding clear-sky values from CM SAF and CERES. Other tests (e.g., the use of AERONET ground-based AOD observations instead of MODIS, and the use of artificial biases in other aerosol and cloud optical properties and atmospheric parameters) did not remove the negative bias. On the other hand, by running our system for the three evaluation stations on an hourly basis (24 simulations per month), we obtained a decreased bias (NMB of -3.5% relative to ground observations, the correlation still being very good). All these give us an indication that the systematic negative bias appearing in our SSR data, and in the data from similar radiative transfer systems [21], is largely due to the fact that our simulations are done on a monthly basis, assuming a monthly mean daily SZA and daylength. Another reason for the negative bias might be the representation of clouds within a day by using afternoon measurements (despite correcting for the diurnal variability). This systematic negative bias is not expected to affect the radiative effect calculations substantially, while, on the other hand, the implementation of monthly simulations is much less time consuming compared to daily or hourly simulations, allowing for more tests and easier fine-tunning of the radiative transfer system. While our study focuses on the average radiative effects on monthly timescales, the instantaneous radiative effects of the parameters examined here could be much higher or lower. For example, Kosmopoulos et al. [99], studying an extreme dust event, showed that the radiative effect induced by dust on SSR was up to −150 W/m 2 over Southern Greece, which is more than 10 times the average all-aerosol radiative effect (−13 W/m 2 ) shown in this work. As far as clouds are concerned, Tzoumanikas et al. [100] showed that the instantaneous cloud radiative effect over Thessaloniki can reach values up to −900 W/m 2 . This is more than 15 times the average radiative effect of clouds that we calculate here for Thessaloniki (−54 W/m 2 ). It is obvious that by using monthly data as the input, one will rule out such extreme values. On top of that, the common practice to use of monthly means as input in the radiative transfer simulations, at some extent, will impact the relative importance of the various parameters, i.e., aerosols, tropospheric NO 2 , and clouds, for the calculated solar radiation fields. For example, under cloudy conditions, satellite-based AOD retrievals are unavailable and hence the AOD monthly mean values come from cloud-free scenes, which could bias the RE aer values. This is an unavoidable limitation of satellite observations and will always exist, despite the fact that we follow the optimal way of aggregating the data here (temporal averaging of spatially aggregated observations), as discussed in Schutgens et al. [101].
Another point that has to be stressed here is that the radiative effect of each parameter (see Equation (1)), except for the parameter per se, depends on the amount of solar radiation reaching the Earth's atmosphere, which in turn varies according to the season of the year and the geographical latitude. This means that even if a parameter was constant throughout the year, the radiative effect would be higher in the summer compared to the winter and over locations at lower latitudes, following the solar radiation variability. This is one of the basic reasons for the deviation observed between the seasonal variability of the parameters and their radiative effects (see Figure 7 and Figures S5-S8). One should have this in mind when comparing radiative effect values from different locations and different seasons. A way to account for this could be to present the "relative" radiative effect values (100 times the radiative effect divided by the corresponding SSR value; in %). For example, for the location of Thessaloniki, RE lc exhibits low values in the winter despite the high liquid cloud fraction, as the incoming solar radiation levels are low compared to the other seasons (Figure 7f). On the other hand, the magnitude of the relative RE lc exhibits a seasonal variability very similar to that of the liquid cloud fraction ( Figure S23). However, this is a climate-oriented study, and hence we select to present the absolute radiative effect values following previous studies in the area [21,22,27,86].

Conclusions
In this work, the direct radiative effect of aerosols and tropospheric NO 2 on SSR is estimated under different sky conditions, along with the radiative effect of different cloud types (liquid and ice), focusing on 16 locations in Greece, in the Eastern Mediterranean. For the scopes of this study, the latest release of various satellite and climatological products is used for the period 2005-2019. The core data are aerosol, cloud, water vapor, and tropospheric/stratospheric NO 2 satellite observations from the MODIS/Aqua C6.1 MYDATML2 and the OMI/Aura QA4ECV v1.1 product. These, and the rest of the datasets, were gridded on a monthly basis at a standard 0.125 • × 0.125 • grid covering the Eastern Mediterranean. The gridded data over the 16 locations are used as the input in a system that was developed based on the SBDART radiative transfer model. The original SBDART code is updated here (SBDART-NO 2 ) so as to account for the absorption of NO 2 in the UV/Vis, and for the tropospheric and stratospheric NO 2 profiles. The main findings of this study are summarized in the following:

•
Our monthly SSR simulations are in good agreement with the ground-based observations from three stations with quality-assured long time series (underestimation of −14.5 W/m 2 or −7.5% with a very good correlation) and two satellite products, CM SAF SARAH-2. Contrary to liquid clouds, the RE ic magnitude exhibits positive trends over all the locations except for one (Rhodes). Argos exhibits the strongest trend (0.46 W/m 2 /year or 3.57%/year), being also the only location with a statistically significant trend. The positive RE ic magnitude trends follow the positive ice cloud fraction trends (ranging from zero to 0.003/year depending on the location); • The liquid and ice cloud optical thickness has increased for some stations, while for others it has decreased during 2005-2019, suppressing or boosting the changes in the radiative effect of clouds, which are induced by changes in the liquid and ice cloud fractions.
This study highlights the driving role of clouds and aerosols in the solar radiation levels reaching the Earth's surface. Tropospheric NO 2 may exert a significant radiative effect on SSR on an annual basis, only over the most heavily polluted areas of the world (e.g., Eastern China). However, under episodic events (e.g., extended forest or agricultural fires), the instantaneous radiative effect of tropospheric NO 2 might be considerable on a local or regional scale. This is something that merits further investigation.