Comparison of Surface Solar Irradiance from Ground Observations and Satellite Data (1990–2016) over a Complex Orography Region (Piedmont—Northwest Italy)

: Climate Monitoring Satellite Application Facility (CM SAF) surface solar irradiance (SSI) products were compared with ground-based observations over the Piedmont region (north-western Italy) for the period 1990–2016. These products were SARAH-2.1 (Surface Solar Radiation DataSet—Heliosat version 2.1) and CLARA-A2 (Cloud, Albedo and Surface Radiation dataset version A2). The aim was to contribute to the discussion on the representativeness of satellite SSI data including a focus on high-elevation areas. The comparison between SSI averages shows that for low OCI (orographic complexity index) stations, satellite series have higher values than corresponding ground-based observations, whereas for high OCI stations, SSI values for satellite records are mainly lower than for ground stations. The comparison between SSI anomalies highlights that satellite records have an excellent performance in capturing SSI day-to-day variability of ground-based low OCI stations. In contrast, for high OCI stations, the agreement is much lower, due to the higher uncertainty in both satellite and ground-based records. Finally, if the temporal trends are considered, average low-elevation ground-based SSI observations show a positive trend, whereas satellite records do not highlight significant trends. Focusing on high-elevation stations, the observed trends for ground-based and satellite records are more similar with the only exception of summer. This divergence seems to be due to the relevant role of atmospheric aerosols on SSI trends. Margherita (4560 m a.s.l.) and Gran Vaudala (3272 m a.s.l.), which presented for ground-based observations values of 171 and 184 Wm -2 , respectively, compared to 128 (139) Wm − 2 and 135 (138) Wm − 2 of SARAH-2.1 (CLARA-A2), for the two sites respectively. The average annual SSI over the ﬁve highest station sites exhibited the same behaviour, with values of 145, 143 and 174 Wm − 2 for SARAH-2.1, CLARA-A2 and the ground stations, respectively.


Introduction
Surface solar irradiance (SSI) is the fraction of solar radiation that reaches the Earth's surface [1]. This variable has a central role in the surface energy balance, driving a large number of physical and chemical processes in the atmosphere [2,3]. It plays an essential role in a wide range of sectors, including agriculture, energy production, tourism, and even health care [4][5][6][7]. SSI, and in particular the proportion between direct and diffuse solar irradiance, also influences global plant productivity and the land carbon sink [8]. Moreover, because solar radiation is a key component of the surface-energy balance, it provides about 75% of the melt energy on most glaciers [9], affecting river runoff and thus freshwater resources fundamental for the human communities. It is therefore very important to improve the datasets and the methods that allow investigating the spatial distribution and the temporal evolution of this variable.
Available studies document a decrease in SSI, referred to as "global dimming", of about 3-9 Wm −2 from the 1950s to the 1980s and a subsequent increase, referred to as "brightening period", of about 1-4 Wm −2 from the beginning of the 1980s [10,11]. This pattern is highlighted both by studies focusing on specific regions [12][13][14][15][16][17] as well as by studies focusing on worldwide datasets [18]. However, a continuous decreasing trend in SSI even after the 2000s has been documented by some studies focusing on countries with developing economies [19].
SSI variations are mostly driven by changes in atmospheric transparency due to variations in cloudiness and/or changes in aerosol concentrations. The relative contribution of clouds and aerosols to SSI variations is however not yet clear, especially when comparing studies regarding different areas and covering different periods. A better understanding of the factors influencing SSI variations is often hampered by the low availability of long-term high-quality station records for both SSI [20,21] and aerosols/clouds especially over oceans, remote or sparsely populated areas, and mountainous regions, where installing and maintaining meteorological stations is more difficult. This problem also reflects a high uncertainty when SSI and its driving factors are estimated in climate models [22].
SSI data from meteorological stations can be integrated with data from satellite remote sensing: on a monthly time-scale, the uncertainty of these data may be in the range of the uncertainty of ground-based observations [23]. Although satellite data have only been available since the 1980s, they have a higher coverage compared to ground-based measurements, especially over remote areas where very often no station data are available. Widely used satellite data for Europe are the gridded SSI fields SARAH-2.1 [24]-an extension in time of SARAH-2 [25]-and CLARA-A2 [26], which are available among climate products from the EUMETSAT Satellite Application Facility on Climate Monitoring (CM SAF-https://www.cmsaf.eu). The validation of these products is extensively discussed both in CM SAF validation reports [24,26] and in scientific papers [27,28]. It highlights an excellent quality of the satellite data showing that they effectively capture SSI spatial and year-to-year variability of corresponding ground stations, although CLARA-A2 records suffer from missing data before 1992 that cause increasing data uncertainty. Moreover, CM SAF SSI products give evidence of a good decadal stability.
An extensive validation of SARAH-JRC (SARAH version used by the European Commission Joint Research Centre for photovoltaic energy calculations) and CLARA-A2 datasets over Europe has also been performed by Urraca et al. [29] for the period from 2005 to 2015.
At EUMETSAT CM SAF, the production of long-term datasets is pursued. The ultimate aim is to make the resulting datasets suitable for the analysis of climate variability, and potentially for the detection of climate trends. It is therefore very important to investigate whether SARAH-2.1 and CLARA-A2 records correctly capture decadal time scale variations such as the slow increase that many ground stations showed in the brightening period. This issue has therefore been investigated extensively by Pfeifroth et al. [28] for Europe: they found a general good agreement even though some areas turn out to be more problematic. Specifically, for northern and southern Europe, trends shown by station data are not completely reproduced by satellite data, especially in summer in southern Europe. For this area, the authors suggested that the underestimation of the trend shown by satellite data could be due to aerosols being considered as constant (monthly climatology) while in the period under analysis they showed a significant reduction. The disagreement between SSI trends captured by SARAH-2 and CLARA-A2 and corresponding trends from ground-based observations in southern Europe has been further investigated by Montero-Martin et al. [30], focusing on 12 stations over the Iberian Peninsula and considering the period 1985-2015. They reported that the two satellite datasets substantially underestimated the trends found for the considered 12 in situ stations in the brightening period. Specifically, although the ground-based SSI data highlighted significant trends at most stations Remote Sens. 2020, 12, 3882 3 of 26 (slopes between −0.5 and 6.5 Wm −2 decade −1 ), the corresponding SARAH-2 and CLARA-A2 data did not show any significant trend: the slopes of the satellite records turned out to be in the intervals between 0.2 and 2.8 Wm 2 decade −1 and −0.4 and 3.8 Wm 2 decade −1 for SARAH-2 and CLARA-A2, respectively. The relevant role of aerosols in the comparison between satellite and ground-based SSI records has also been highlighted by Wang et al. [31] for China. They reported a rather large bias for both SARAH-E (data from Meteosat East satellites) and CLARA-A2 records that largely decreased if the comparison between satellite and ground-based records was performed for rural stations only.
Within this context, the present work aims to provide a new contribution to the comparison between SSI data from ground-based observations and CM SAF SARAH-2.1 and CLARA-A2 products. It focuses both on low-and high-elevation sites where satellite datasets can be fundamental because of the low availability of station data. Specifically, the work focuses on an area in north-western Italy which is covered by a very dense network of SSI stations covering the 1990-2016 period that has been subjected to a detailed homogenisation by Manara et al. [32]. Besides its excellent data availability, this area is also interesting because it includes stations with considerably different levels of aerosol concentrations, ranging from the highly polluted Po Plain to the unpolluted Alpine high-elevation sites. Moreover, in this area the emission of pollutants has strongly decreased over the last decades, causing a significant trend in all variables connected to atmospheric transparency with the only exception of high-level stations that fall outside the atmospheric layer directly influenced by ground-level pollutant emissions [14,15,[32][33][34]).
The paper is organised in five sections. After the Introduction, Section 2 presents the databases used and the methods adopted for the analyses. Then, Section 3 presents the main results, Section 4 discusses the results presented in Section 3, and finally some conclusive remarks are given in Section 5.

Ground-Based Observations
In this work we used a subset of 57 records extracted from the dataset of 1990-2016 ground-based SSI observations presented by Manara et al. [32]. The considered stations were in the area 6.6 • -9.0 • E, 44.1 • -46.2 • N and they had elevations between 77 m and 4560 m. All the stations were included in the meteorological network of the Piedmont environmental agency (Arpa Piemonte-http://www.arpa. piemonte.it/) ( Figure 1 and Table 1).
Remote Sens. 2020, 12, 3882 4 of 27 2016, whereas some gaps remained for the period 1990-2000 (about 0.4% of the data considering only the series selected for this paper) because for some mountain stations adequate reference series were not always available within the gap filling procedure. Station metadata give evidence that some station locations do not fulfil World Meteorological Organization requirements because they are affected by shading in the first hours after sunrise or before sunset. This problem induces significant bias in some station records; it mainly affects some urban stations and some stations in mountain areas. Another issue that has to be considered when station records are compared to satellite data is the homogenisation procedure: in fact, as explained above, this procedure is generally applied to station data with the only goal of removing breaks and it does not take into account that if a biased period is considered as the reference period, then it will increase the bias of a record. In case of a break, homogenisation methods can in fact either adjust the data before the break or those after the break (independently from which of them are correct) and the resulting records, although different in terms of absolute values, are equivalent for the assessment of temporal trends. Therefore, when focusing on the comparison between stations and satellite records, station metadata (including information on the applied adjustments) should be considered in order to correctly interpret the obtained results. CLARA-A2 grid points (blue circles). The filled circles identify the grid points considered for the analyses. The orography of the region is represented too. Table 1. Coordinates and elevation of the 57 station sites ordered according to the maximum between SARAH-2.1 and CLARA-A2 orographic complexity indexes (OCI). The table also shows the associated orographic complexity index and the yearly average surface solar irradiance (SSI) in the period 2001-2016 from the ground-based, the SARAH-2.1, and CLARA-A2 records. In the last column, it is reported if a station's record is considered in Figures 2 and 3.  The data were recorded with first-class radiometers usually checked for calibration one or two times per year by means of a reference instrument (see Manara et al. [32] for full details). After preliminary quality checks, the station records were subjected by Manara et al. [32] to a gap-filling procedure in Remote Sens. 2020, 12, 3882 5 of 26 order to estimate daily missing values. This procedure aimed to avoid bias in the monthly averages due to gaps in the daily series. Specifically, each daily missing value was estimated by means of the most correlated reference series using a relationship based on the assumption of the constancy of the ratio between the test series (the series to be filled) and the reference series. Moreover, the reference series should belong to the same climatic area of the test series and the data that are missing in the test series should be available.
The monthly station records were then subjected by Manara et al. [32] to a homogenisation procedure in order to identify and eliminate non-climatic signals. In particular, the Craddock test was applied, comparing the test series with a reference series to identify possible breaks. In order to be sure that the break was ascribable to the test series and not to the reference series, the Craddock test [35,36] was applied to more than one reference series. Then, the adjustment factors to apply to the test series were calculated using only those reference series that proved to be homogeneous in a sufficiently long period centred on the break. It should be noted that the adjustment factors are usually applied to the part of the series before the break (the old one), considering the part after the break as homogeneous (the more recent one) in order to simplify the updating of the series when new data are available. This approach assures the homogeneity in terms of temporal evolution but not in terms of normal values. The monthly correcting factors were also transferred to the daily series in order to make them homogeneous as well. Finally, we subjected the monthly records to a procedure aiming at completing the records over the entire investigated period [32] using a methodology similar to that adopted for the daily series. The data filling procedure enabled the completion of all the records over the period 2001-2016, whereas some gaps remained for the period 1990-2000 (about 0.4% of the data considering only the series selected for this paper) because for some mountain stations adequate reference series were not always available within the gap filling procedure.
Station metadata give evidence that some station locations do not fulfil World Meteorological Organization requirements because they are affected by shading in the first hours after sunrise or before sunset. This problem induces significant bias in some station records; it mainly affects some urban stations and some stations in mountain areas. Another issue that has to be considered when station records are compared to satellite data is the homogenisation procedure: in fact, as explained above, this procedure is generally applied to station data with the only goal of removing breaks and it does not take into account that if a biased period is considered as the reference period, then it will increase the bias of a record. In case of a break, homogenisation methods can in fact either adjust the data before the break or those after the break (independently from which of them are correct) and the resulting records, although different in terms of absolute values, are equivalent for the assessment of temporal trends. Therefore, when focusing on the comparison between stations and satellite records, station metadata (including information on the applied adjustments) should be considered in order to correctly interpret the obtained results.

SARAH Dataset
SARAH, the Surface Solar Radiation DataSet-Heliosat, is a satellite-based database of several variables linked to surface solar radiation. SARAH first edition (SARAH-1.0) was issued by EUMETSAT CM SAF in 2015 and it covered the period 1983-2013, whereas SARAH-2.0 (1983-2015) was issued in 2017. SARAH-2.1 was then issued in 2019, updating SARAH-2.0 to 2017. The surface solar radiation records are obtained after the combination of the following two steps [37]: (i) the observations obtained from the visible channels of the Meteosat Visible Infra-Red Imager (MVIRI) and the Spinning Enhanced Visible and Infrared Imager (SEVIRI) instruments on board EUMETSAT's geostationary Meteosat satellites are processed using a climate version of the Heliosat algorithm to obtain the effective cloud albedo [38]. Specifically, the Heliosat algorithm uses the reflectance or radiance measurements to determine the effective cloud albedo comparing the satellite measurement and the corresponding observation under clear-sky conditions. A brighter pixel suggests that a higher fraction or thicker clouds are present. This algorithm can present some problems, for example over bright surfaces such as snow; (ii) clear sky radiation is then calculated by means of the clear sky part of the SPECMAGIC method (SPECtral Mesoscale Atmospheric Global Irradiance Code-http://gnu-magic.sourceforge.net/) [39]. This method is based on an LUT (lookup table) approach where discrete pre-computed results of a radiative transfer model are stored in order to reduce runtime computation. Then, the value for a given atmospheric state can be extracted for each satellite pixel and time and combined with the extra-terrestrial incoming solar flux density. The main processes considered are Rayleigh scattering, scattering and absorption by aerosols (aerosol monthly climatology at 1 • × 1 • resolution), the absorption by water vapor (reanalysis monthly values at 0.75 • × 0.75 • resolution, topographically corrected to 0.15 • resolution), ozone, and the reflection by the Earth's surface.
SARAH-2.1 data cover the region between ±65 • longitude and ±65 • latitude. SSI records can be downloaded from CM SAF web site as monthly and daily means, and as 30 min data on a regular latitude/longitude grid with a spatial resolution of 0.05 • × 0.05 • degrees. As stated in the validation report [40], the comparison with the Baseline Surface Radiation Network (BSRN) data of the monthly means shows a mean bias of 2 Wm −2 , a mean absolute error of 5.1 Wm −2 , and a correlation of the monthly anomalies of 0.93. For the daily means, the mean bias is comparable to the value derived for the monthly means while the mean absolute error is 11.8 Wm −2 .
In this study, we used only the SARAH-2.1 grid points closest to the 57 station sites ( Figure 1). They are 54 grid points, because 3 grid points were associated to more than one station site. For each SARAH-2.1 grid point, we also considered all the grid points of the 30 arc-second resolution GTOPO30 digital elevation model [41] falling within the SARAH-2.1 grid box and we calculated the average and the standard deviation of their elevations. The grid point standard deviation was then used as an index representing the orographic complexity (OCI) associated with the considered SARAH-2.1 grid point, and finally the SARAH-2.1 grid point OCIs were used, together with corresponding CLARA-A2 grid point indexes (see next section), to sort the stations in Table 1.
SARAH-2.1 records have a very low fraction of missing data in the investigated period (about 0.2% for daily data). However, we also subjected these records to the same data-filling procedure applied to ground-based observations aiming at obtaining monthly records which are not biased by missing values. This procedure enabled the completion of all the records over the considered period.

CLARA Datset
CLARA, the Cloud, Albedo and Surface radiation dataset from AVHRR data, is another satellitebased database of several variables linked to surface solar radiation. CLARA first edition (CLARA-A1) was issued by EUMETSAT CM SAF in 2012 and it covered the period 1982-2009, whereas CLARA-A2  was issued in 2017 [26]. The CLARA-A2 dataset, an improved and extended version of CLARA-A1, is based on data from the advanced very high resolution radiometer (AVHRR) on board the polar orbiting NOAA satellites, as well as on board the MetOp polar orbiters operated by EUMETSAT since 2006. To estimate the SSI from the AVHRR measurements, the cloud mask based on the NWC SAF polar platform system (PPS) cloud processing software [42] is used to distinguish clear and cloudy pixels. For clear-sky pixels, the SSI is derived using auxiliary parameters (e.g., water vapor, surface albedo) with the MAGIC (mesoscale atmospheric global irradiance code) clear-sky radiative transfer model [43]. Under cloudy situations, broadband solar fluxes are estimated using the visible reflectance as measured by the AVHRR instrument. A look-up table (also depending on auxiliary data) is used to relate these reflected solar fluxes to the atmospheric transmissivity, which is transferred to SSI using the extra-terrestrial solar flux. The daily mean SSI is derived from the instantaneous satellite information using the method from Diekmann et al. [44], which considers the diurnal solar cycle for the averaging calculations. For the spatial averaging from the resolution of the satellite pixel (approx. 0.05 • ) to the 0.25 • grid, a minimum of 20 valid daily means are required. The monthly means are only calculated if more than 20 valid daily means are available. More detailed information can be found in Karlsson et al. [45] and references therein.
Remote Sens. 2020, 12, 3882 7 of 26 CLARA-A2 data are available globally with a spatial resolution of 0.25 • × 0.25 • and at daily and monthly temporal resolutions. They cover the period 1982-2018, even if the initial period has a rather large number of missing data because only one satellite was available until the early 1990s. As stated in the validation report, CLARA-A2 shows a bias of −1.6 Wm −2 at monthly resolution when compared against the BSRN data, a mean absolute error of 8.8 Wm −2 , and a correlation of the anomalies of 0.87. Similarly, at daily resolution CLARA-A2 shows a bias of −1.7 Wm −2 , a mean absolute error of 18.6 Wm −2 and a correlation of the anomalies of 0.9. As for SARAH-2.1, for CLARA-A2 we used only the grid points that were closest to the 57 station sites (Table 1). They were 34 grid points, because 15 grid points were associated to more than one station site. Moreover, for CLARA-A2 too, we used the GTOPO30 digital elevation model to obtain, for each considered grid point, the average elevation and the OCI using the same procedure adopted for SARAH-2.1 (see Table 1). When sorting the station sites for orographic complexity, we used the highest value between the CLARA-A2 and SARAH-2.1 indexes. It was usually the OCI from CLARA-A2, because of the lower resolution of the CLARA-A2 grid.
CLARA-A2 records had a greater fraction of missing data than SARAH-2.1 (about 5.4% for daily data) and this fraction was particularly relevant before 2001. We subjected the records to the data-filling procedure, already used in Sections 2.1 and 2.2, aiming at obtaining monthly records which were not biased by missing values. This procedure enabled the completion of all the records over the period 2001-2016, whereas a relevant fraction of gaps remained for the period 1990-2000 (12.4% of the data) because adequate reference series were not available within the gap filling procedure.

Other Datasets
As well as the ground-based stations and satellite SSI datasets, we also used Atmospheric Optical depth data from the Modern-Era Retrospective Analysis for Research and Applications, Version 2 (MERRA-2-https://gmao.gsfc.nasa.gov/reanalysis/MERRA-2/data_access/ [46]) to investigate the evolution of atmospheric transmissivity over the considered period. This atmospheric reanalysis, produced by NASA's Global Modelling and Assimilation Office (GMAO), has a horizontal resolution of 0.5 • × 0.625 • . Moreover, we used 500 m-resolution climatologies of length of ground snow coverage per year calculated from MODIS satellite data using the methodology described by Gafurov and Bardossy [47]. These values have been used to calculate corresponding climatologies for SARAH-2.1 and CLARA-A2 grid cells. Specifically, for each SARAH-2.1 and CLARA-A2 grid cell, the corresponding length of snow coverage is calculated as the average of all MODIS grid cells falling within the selected SARAH-2.1 and CLARA-A2 grid.

Error Indexes and Trend Analysis
The comparison of SSI values from the different data records is performed by first analysing the average differences over the period 2001-2016 between ground-based and satellite observations, and of the two satellite datasets (Section 3.1). Specifically, the comparison was performed focusing on the entire dataset when the two satellite datasets were compared and on a subset of 32 station sites (last column of Table 1) when ground-based and satellite records were considered. In both cases, the records were clustered in two groups according to their OCI (stations with OCI lower than or equal to 153 m and stations with OCI higher than 153 m).
We selected a subset of series to compare ground-based and satellite observations (both SARAH-2.1 and CLARA-A2) because a significant fraction of the ground-based records, even though homogenous in terms of temporal trend, may have had biased SSI values in terms of normal values (see Section 2.1). Indeed, we eliminated the records falling in at least one of the following three categories: • records with metadata giving evidence of shading; • records for which the 2001-2016 SSI normal values changed by more than 5% because of the homogenisation procedure; Remote Sens. 2020, 12, 3882 8 of 26 • records with seasonal average SSI below the following thresholds: 60 Wm −2 for winter, 180 Wm −2 for spring, 240 Wm −2 for summer, 100 Wm −2 for autumn, and 150 Wm −2 for the year. These thresholds were considered because metadata were not complete for all stations and it was therefore necessary to add a further condition to identify the records which may have been affected by shading problems not documented by the metadata.
Twenty-five records were excluded by these criteria, while 32 records were evaluated as unbiased in terms of normal values. For these 32 records, we also checked the effect of shading from the orography falling within a distance of 10 km by means of a 30 m-resolution Digital Elevation Model (DEM), resampled to 90 m to increase computation speed [48]. Specifically, for each station site we considered the days falling in the middle of each month and we first divided the period between sunrise and sunset into 48 parts. Then, for each of them we verified whether direct radiation could reach the station and finally we used this result within a clear-sky model, based both on Lambert-Beer's law and on a simple estimation of diffuse radiation (see Manara et al. [33] for full details on this model), to estimate by how much shading can reduce SSI. The results highlighted that, even with a clear sky (i.e., the condition more critical for shading issues), in no one of the station sites did topographic shading reduce SSI by more than 3%. However, any shading effect can affect only climatological values and not the year-to-year variability nor long term trend.
All SSI records were then transformed into daily, monthly, seasonal, and annual anomaly records (additive anomalies) with respect to the corresponding normals (daily, monthly, seasonal, annual) of the period 2001-2016 and these anomaly records were used to investigate the agreement of the variability between SSI from ground-based and satellite observations (Section 3.2). The agreement was investigated focusing on mean absolute differences as well as on common variances, considering different temporal resolutions (from the day to the year) and clustering the entire dataset (57 stations) and the selected subset of records (32 stations) according to the station OCI. The mean absolute differences, usually indicated as station-site satellite mean absolute errors, give evidence of the accuracy of the satellite records, whereas the common variances show how much of the variability of the ground-based records is captured by the satellites.
The aim of Sections 3.1 and 3.2 is to analyse: (i) the bias between the different data sources considering that both type of measures can be affected by errors; and (ii) to analyse the agreement of the anomaly records at different temporal resolutions, investigating in the first case whether a systematic difference exists between the different datasets, and in the latter how much of the variability explained by the ground-based observations can be reproduced by satellite observations. Finally, the SSI data records from the different data sources are compared for trends (Section 3.3). Specifically, monthly, seasonal and annual anomaly records representative of low-elevation station sites (elevation lower than or equal to 600 m-39 station sites) and high-elevation station sites (station elevation higher than 600 m-18 station sites) are calculated by simply averaging the corresponding monthly, seasonal and annual anomaly records. The trends are then evaluated with the Theil-Sen [49,50] method and the corresponding significances are evaluated by means of the Mann-Kendall non-parametric test [51]. Specifically, the Theil-Sen estimator considers the median of the slopes of all lines passing through each pair of points. With respect to the simple linear regression, it is less sensitive to outliers. The Mann-Kendall non-parametric test is used to determine whether a time series has an upward or downward trend, and being a non-parametric test it can be used for all distributions. For trends, we focused on station site elevations rather than on OCI because we aimed at highlighting whether the agreement of the temporal behaviour shown by the two types of observations depends on the considered altitude characterised by different cloudiness conditions and aerosol concentrations.

Comparison of Average SSI Values from Different Data Records
SARAH-2.1 data gave evidence of higher yearly average SSI values over the 57 considered sites than CLARA-A2 ones, with an average difference of 6.7 Wm −2 over the period 2001-2016 ( Table 2). The seasons with the highest differences between SARAH-2.1 and CLARA-A2 average SSI values were spring and summer (8.6 and 8.3 Wm −2 , respectively), whereas in winter the difference reduced to 3.9 Wm −2 ( Table 2). Focusing only on the 30 station sites with low orographic complexity (maximum OCI = 153 m), the average difference between SARAH-2.1 and CLARA-A2 SSI values was very stable along the year and all seasonal averages fell within the interval 5.0-5.4 Wm −2 ( Table 2). These station sites also provided evidence of a rather low spatial variability of the differences between SARAH-2.1 and CLARA-A2 SSI values. Considering the yearly averages of these differences, we obtained a standard deviation of 1.4 Wm −2 , whereas considering the seasonal averages the standard deviations ranged from 1.0 Wm −2 in winter to 2.9 Wm −2 in summer ( Table 2). In contrast, focusing on the remaining 27 station sites (171 m ≤ OCI ≤ 702 m), the differences between SARAH-2.1 and CLARA-A2 SSI values exhibited a much less coherent pattern both in terms of seasonality and in terms of spatial distribution. In this case, the seasonal average station site differences ranged from 2.6 Wm −2 in winter to about 12 Wm −2 in spring and in summer, whereas the standard deviations of the station site average differences ranged from 6.6 Wm −2 in autumn to 14.0 Wm −2 in summer, with a rather high value on the yearly scale as well (10.5 Wm −2 - Table 2). The results of the comparison between satellite records and ground-based records using the subset of 32 records evaluated as unbiased in terms of normal values (see last column of Table 1 and Section 2.5) are shown in Figures 2 and 3. They show that for station sites with low OCI (OCI ≤ 153 m), the differences between SSI values recorded by ground stations and corresponding satellite values were rather stable: SARAH-2.1 SSI annual averages were 6.9 ± 2.4 Wm −2 higher than corresponding values from ground stations, whereas for CLARA-A2 this value decreased to 1.7 ± 2.3 Wm −2 . The average differences between SSI values recorded by satellite values and corresponding ground stations were lowest in winter (5.1 ± 2.6 and 0.2 ± 2.5 Wm −2 for SARAH-2.1 and CLARA-A2, respectively) and highest in autumn (9.1 ± 1.9 and 3.7 ± 1.4 Wm −2 for SARAH-2.1 and CLARA-A2, respectively). For station sites with high OCI, the situation was completely different and Figures 2 and 3 give evidence of a complex pattern; most station sites had higher SSI values for ground-based observations and a few others had higher values for satellite data.      Figure 4 provides a quick overview of the differences between the spatial patterns of yearly average SSI values from the three discussed data sources. As for Figures 2 and 3, also in this case, for ground-based observations we focused only on the selected subset of 32 stations.  Table 1) are shown too.

Low OCI Station Sites
Considering the station sites characterised by low OCI, both for daily and monthly resolution, the mean absolute differences between SARAH-2.1 and ground-based records, with values of 11.8 Wm −2 and 5.0 Wm −2 (considering all data), respectively, were very close to SARAH-2.1 mean absolute errors obtained by the comparison with the BSRN dataset [40]. At monthly resolution, CLARA-A2 records showed mean absolute differences from ground-based records (5.2 Wm −2 considering all months), comparable to those obtained between SARAH-2.1 and ground-based records, whereas at daily resolution the mean absolute differences were larger (15.2 Wm −2 considering all days), even though the daily mean difference obtained considering all days turned out to be lower than CLARA-A2 mean absolute errors obtained by the comparison with the BSRN dataset (18.6 Wm −2 ) [52]. Similarly, the common variances between monthly satellite and ground-based records (84% and 83% for SARAH-2.1 and CLARA-A2, respectively, considering all months-95% confidence interval (CI95) ranging from 80% to 87% for SARAH-2.1 and from 80% to 86% for CLARA-A2) were rather similar to those presented in SARAH-2.1 and CLARA-A2 validation reports. The difference between SARAH-2.1 and CLARA-A2 performances, when compared with ground-based records, was, however, lower than that reported by the comparison with the BSRN dataset. Moving from monthly to daily resolution, the common variances increased (92.1% (CI95: 91.7-92.4%) and 87.5% (CI95: 86.9-88.0%) for SARAH-2.1 and CLARA-A2 considering all days), giving evidence of an excellent performance of the daily satellite records, especially of SARAH-2.1, in capturing SSI variability measured by ground stations. On the contrary, moving from monthly to seasonal resolution, common variances decreased with seasonal records showing common variances equal to 79% (CI95: 71-85%) and 79% (CI95: 70-85%) when all seasons were considered for SARAH-2.1 and CLARA-A2, respectively. The common variances further reduced to 73% (CI95: 50-86%) and 69% (CI95: 39-86%) for SARAH-2.1 and CLARA-A2, respectively, when annual resolution was considered.
If the analysis was limited to the 18 low OCI station sites selected from the subset of 32 stations we identified in Section 2.4 (data not shown in Table 3), the agreement between satellite and groundbased records slightly increased and monthly (daily) mean absolute differences reduced to 4.9 Wm −2 (11.2 Wm −2 ) for SARAH-2.1 and to 5.1 Wm −2 (14.6 Wm −2 ) for CLARA-A2. In this case, common variances also slightly increased. The improvement in the agreement between satellite and groundbased records was, however, very low, highlighting that the full set of station sites seems to be adequate to study the coherence of SSI anomaly records, even if it is not completely appropriate to investigate the agreement between their average SSI values as discussed in Section 2.5.  Table 1) are shown too. Table 3 shows the common variances between ground-based anomaly records (considering all investigated sites) and satellite records, averaged over low and high OCI station sites and considering different temporal resolutions, from the day to the year. Moreover, it shows the corresponding mean absolute differences. Table 3. The table shows the common variances and the mean absolute differences (MAD) between satellite and ground station datasets using the 57 investigated sites. Both common variances and mean absolute differences are calculated separately for low and high OCI station sites for different periods and considering daily, monthly, seasonal, and annual anomaly records. The anomalies are calculated with respect to the period 2001-2016.

Low OCI Station Sites High OCI Station Sites
Common Variance MAD (Wm −2 ) Common Variance MAD (Wm −2 )

Low OCI Station Sites
Considering the station sites characterised by low OCI, both for daily and monthly resolution, the mean absolute differences between SARAH-2.1 and ground-based records, with values of 11.8 Wm −2 and 5.0 Wm −2 (considering all data), respectively, were very close to SARAH-2.1 mean absolute errors obtained by the comparison with the BSRN dataset [40]. At monthly resolution, CLARA-A2 records showed mean absolute differences from ground-based records (5.2 Wm −2 considering all months), comparable to those obtained between SARAH-2.1 and ground-based records, whereas at daily resolution the mean absolute differences were larger (15.2 Wm −2 considering all days), even though the daily mean difference obtained considering all days turned out to be lower than CLARA-A2 mean absolute errors obtained by the comparison with the BSRN dataset (18.6 Wm −2 ) [52]. Similarly, the common variances between monthly satellite and ground-based records (84% and 83% for SARAH-2.1 and CLARA-A2, respectively, considering all months-95% confidence interval (CI95) ranging from 80% to 87% for SARAH-2.1 and from 80% to 86% for CLARA-A2) were rather similar to those presented in SARAH-2.1 and CLARA-A2 validation reports. The difference between SARAH-2.1 and CLARA-A2 performances, when compared with ground-based records, was, however, lower than that reported by the comparison with the BSRN dataset. Moving from monthly to daily resolution, the common variances increased (92.1% (CI95: 91.7-92.4%) and 87.5% (CI95: 86.9-88.0%) for SARAH-2.1 and CLARA-A2 considering all days), giving evidence of an excellent performance of the daily satellite records, especially of SARAH-2.1, in capturing SSI variability measured by ground stations. On the contrary, moving from monthly to seasonal resolution, common variances decreased with seasonal records showing common variances equal to 79% (CI95: 71-85%) and 79% (CI95: 70-85%) when all seasons were considered for SARAH-2.1 and CLARA-A2, respectively. The common variances further reduced to 73% (CI95: 50-86%) and 69% (CI95: 39-86%) for SARAH-2.1 and CLARA-A2, respectively, when annual resolution was considered.
If the analysis was limited to the 18 low OCI station sites selected from the subset of 32 stations we identified in Section 2.4 (data not shown in Table 3), the agreement between satellite and ground-based records slightly increased and monthly (daily) mean absolute differences reduced to 4.9 Wm −2 (11.2 Wm −2 ) for SARAH-2.1 and to 5.1 Wm −2 (14.6 Wm −2 ) for CLARA-A2. In this case, common variances also slightly increased. The improvement in the agreement between satellite and ground-based records was, however, very low, highlighting that the full set of station sites seems to be adequate to study the coherence of SSI anomaly records, even if it is not completely appropriate to investigate the agreement between their average SSI values as discussed in Section 2.5.
Considering the annual cycle of the common variances between ground-based and satellite records and focusing both on daily and monthly resolution, a clear minimum was evident in winter, the season in which air pollution has a major role. Focusing on seasonal records, in winter there was also a rather low common variance between ground-based and satellite SSI anomalies, with values of 64% (CI95: 38-82%) and 60% (CI95: 28-82%), for SARAH-2.1 and CLARA-A2, respectively. Winter therefore turned out to be the most problematic season in terms of the agreement between ground-based and satellite SSI records.

High OCI Station Sites
Considering the station sites characterised by high OCI, the common variance between satellite SSI anomalies and ground-based ones was much lower than that obtained for stations with low OCI, as expected due to the higher uncertainty both in satellite and ground-based records. In this case, the common variance between satellite and ground-based SSI anomalies was generally higher for CLARA-A2 than for SARAH-2.1. This behaviour was particularly evident in winter, where SARAH-2.1 and ground-based records exhibited common variances lower than 50% at all temporal resolutions. The better agreement of CLARA-A2 SSI anomaly records with corresponding ground-based records for high OCI station sites was also evident from lower mean absolute differences than those observed for SARAH-2.1.
For high OCI station sites, the improvement in the agreement between satellite and ground-based records obtained when the analysis was limited to the sites we identified in Section 2.5 was higher than for low OCI station sites (data not shown in Table 3), especially in winter.

Trends of SSI Records from Different Data Sources
At low-elevation station sites, the average SSI record from ground-based observations showed a clear positive trend in the period 1990-2016 ( Figure 5 and Table 4). This trend was evident (p-value ≤ 0.05) in spring (5.7 Wm 2 decade −1 ), summer (6.4 Wm 2 decade −1 ) and autumn (4.1 Wm 2 decade −1 ), as well as when considering monthly and yearly data (4.0 Wm 2 decade −1 and 3.8 Wm 2 decade −1 , respectively), whereas in winter, the sign of the trend was negative, but the trend was not significant. Conversely, the average SSI records from satellite observations (both SARAH-2.1 and CLARA-A2- Figure 5 and Table 4) did not highlight significant positive trends (p-value ≤ 0.05). Considering a lower threshold for trend significance (p-value ≤ 0.1), there was evidence of a positive trend only for the average autumn low-elevation SARAH-2.1 record (3.4 Wm −2 decade −1 ), whereas for the seasonal CLARA-A2 records there was no evidence of significant positive trends, with a significant negative trend in winter (−5.5 Wm −2 decade −1 ). These results generally did not change when the average SSI records were calculated focusing only on the 29 station sites that, besides being at low-elevation, also had a low OCI ( Table 4). The only relevant difference concerned the positive trend of the summer CLARA-A2 record that became significant (p-value ≤ 0.05). Moreover, in this case the trend of the SARAH-2.1 monthly record became positive when a p-value ≤ 0.1 was considered as the significance threshold.
The average SSI record from ground-based observations generally had a positive trend also focusing on high-elevation station sites ( Figure 6 and Table 4), even though the trend was significant only in summer (6.3 Wm 2 decade −1 ) and autumn (5.0 Wm 2 decade −1 ) and the increase observed for the annual record (3.1 Wm 2 decade −1 ) was lower than at low-elevation station sites. For high-elevation station sites, in autumn there was also a positive trend for the average satellite records (both SARAH-2.1 and CLARA-A2) and the SARAH-2.1 record also had a positive trend in spring (p-value ≤ 0.1). As for low-elevation sites, the CLARA-A2 record gave evidence of a significant negative trend in winter (−4.5 Wm 2 decade −1 ).    The different behaviour of ground-based and satellite low-elevation SSI anomaly records is highlighted even more clearly by focusing on the differences between them (Figures 7-10 and Table 5). Average annual and seasonal series (thin lines) of the differences between ground-based and SARAH-2.1 low-elevation SSI anomaly records. The averages were calculated for both lowelevation station sites and low-elevation station sites with a low OCI. In the case of records without gaps, the series were plotted together with an 11-year window, 3-year standard deviation Gaussian low-pass filter (thick lines).

Figure 7.
Average annual and seasonal series (thin lines) of the differences between ground-based and SARAH-2.1 low-elevation SSI anomaly records. The averages were calculated for both low-elevation station sites and low-elevation station sites with a low OCI. In the case of records without gaps, the series were plotted together with an 11-year window, 3-year standard deviation Gaussian low-pass filter (thick lines).

Figure 8.
Average annual and seasonal series (thin lines) of the differences between ground-based and CLARA-A2 low-elevation SSI anomaly records. The averages were calculated for both lowelevation station sites and low-elevation station sites with a low OCI. In the case of records without gaps, the series were plotted together with an 11-year window, 3-year standard deviation Gaussian low-pass filter (thick lines).

Figure 8.
Average annual and seasonal series (thin lines) of the differences between ground-based and CLARA-A2 low-elevation SSI anomaly records. The averages were calculated for both low-elevation station sites and low-elevation station sites with a low OCI. In the case of records without gaps, the series were plotted together with an 11-year window, 3-year standard deviation Gaussian low-pass filter (thick lines).
for SARAH-2.1 data these ratios show a rather different yearly cycle with respect to ground-based observations with up to 25% lower winter SSI values for high-elevation sites with respect to lowelevation ones. For CLARA-A2 data, some differences are evident too with a less pronounced yearly cycle than for ground-based observations and with a minimum falling in spring.
It will be interesting in the future, when the data become available, to compare the Piedmont ground-based series with SARAH-3 data records, the new version of SARAH-dataset that will also include a scheme to better identify and consider snow coverage. In addition, the spatial representativity of the station measurements also needs to be considered when comparing local measurements with gridded datasets, which always represent a larger area (25 km 2 and 625 km 2 in the case of SARAH-2.1 and CLARA-A2, respectively) especially when areas with complex orography are considered [54]. However, the spatial resolution of the considered grid should be considered together with the temporal resolution of the data. In fact, for temporal resolution higher than the day, if the spatial resolution is of a few kilometres (e.g., 5 × 5 km) the representativeness error can be considered smaller. Differently, the representativeness error cannot this period in Italy, affected atmospheric transparency more strongly at low-than at high-elevations. This effect is extensively discussed for the entire Italian territory by Manara et al. [34] in a study focusing on atmospheric visibility records. It is therefore not surprising that satellite records do not give evidence of the same brightening signal of ground-based observations. Interestingly, the differences in the trends of SSI records from ground-based and satellite observations are more evident at low-elevation, whereas high elevation sites show a similar behaviour to low-elevation sites only in summer: in this season the boundary layer is highest and air pollution concerns a much deeper layer than in the other seasons [56][57][58].
The results obtained in this study are in agreement with those reported in literature even if they  High-elevation station sites 2.4 ± 0.5 + + 4.9 ± 1.9 1.5 ± 0.8 + For low-elevation station sites (Figures 7 and 8), the records of these differences gave evidence of a clear increase in all seasons as well as for the monthly and yearly data. Specifically, the trend of the differences between ground-based and CLARA-A2 anomalies exhibited a p-value ≤ 0.01 in all cases, whereas when the differences between ground-based and SARAH-2.1 anomalies were considered, trend significance was lower in autumn and the trend was not significant in winter. Interestingly, in this case, focusing only on the 29 station sites that, besides being at low-elevation, also have a low OCI, gave slightly more significant trends for the SARAH-2.1 average records, with an autumn p-value ≤ 0.01 and winter p-value very close to the 0.05 threshold. Considering the differences between the ground-based and the satellite yearly anomaly records, the trend was slightly lower than 3 Wm 2 decade −1 for SARAH-2.1 and slightly higher than 3 Wm 2 decade −1 for CLARA-A2. These slopes caused a variation of the differences between ground-based and satellite SSI of about 7 Wm −2 for SARAH-2.1 and about 8 Wm −2 for CLARA-A2, that correspond to about 5% of the yearly average SSI over the investigated area.
Focusing on high-elevation station-sites (Figures not shown), the differences between the trends of ground-based and satellite average SSI anomaly records were generally lower than for low-elevation station-sites (Table 5). In this case, the season with the strongest differences in the trends was summer. In this season, the trends of the differences between ground-based and satellite anomaly records exhibited p-values ≤ 0.01.

Discussion
The analysis of the differences between the different datasets showed that on one side, the comparison between the two satellite datasets (Section 3.1) over the period 2001-2016 gives evidence of higher yearly SSI averages for SARAH-2.1 than for CLARA-A2. This problem is less evident for station sites with low OCI (orographic complexity index-Sections 2.2, 2.3 and 2.5) while it is more evident for those with high OCI. The latter station sites also give evidence of a much less coherent pattern both in terms of seasonality and in terms of spatial distribution. On the other hand, starting from a subset of ground-based observations considered un-biased in terms of normal values (Section 2.5), the comparison between satellite records and corresponding ground-based records showed that for the low OCI station sites, the satellite series have higher annual SSI averages than ground-based series. In contrast, for the high OCI station sites the same comparison showed a complex pattern, with most station sites having lower SSI values in the satellite records and a few others having lower values for ground-based observations. The most important factor decreasing SSI for satellite records seems to be snow coverage. Underestimation of satellite SSI is, in fact, particularly pronounced for the sites with highest elevation (e.g., for Capanna Margherita and Gran Vaudala stations the average lengths of ground snow coverage per year calculated using the nearest grid-point from MODIS satellite data are 365 and 311 days, respectively) and in winter and spring (in this season, for example, SSI averages over the five highest station sites are 159, 160 and 215 Wm −2 for SARAH-2.1, CLARA-A2 and ground stations, respectively) when surface albedo is maximum. On the contrary, the underestimation is less evident in autumn and summer when albedo decreases because snow coverage is lower. The effect of ground snow coverage on SARAH-2.1 and CLARA-A2 SSI is also evident (Figure 9) when satellite yearly average SSI (averages over the 2001-2016 period- Figure 4) is plotted versus the corresponding average length of snow coverage (see Section 2.4). The problem of snow coverage is that it is not always correctly distinguished from clouds, resulting in an enhanced cloud optical depth and subsequently in an underestimation of satellite SSI records, in particular under clear-sky conditions [28,53]. At the same time, high surface albedo from surrounding slopes can probably bias ground observations too because of reflected radiation increasing the measured values. The combination of these two effects therefore causes very high differences between SSI ground-based observations and corresponding satellite values for some station sites. Ground-based observations can however also have significantly lower SSI than satellite ones for other reasons, e.g., because of shading issues.
Despite potentially relevant problems, areas with complex orography ground-based observations show a more consistent behaviour than satellite ones. This is evident from the station-site average yearly ( Figure 4) and seasonal (Figures not shown) values over the considered period. More evidence of this can be obtained by plotting the yearly cycles of the ratios between SSI averages over high-elevation station sites and low-elevation station sites with low OCI: we calculated these ratios for the three datasets over the 2001-2016 period focusing on the sites corresponding to the subset of 32 stations considered un-biased in terms of normal values ( Figure 10). For ground-based observations, as expected, these ratios showed higher values in winter when SSI at low-elevation is often reduced by low atmospheric transparency in the boundary layer, whereas in summer low-and high-elevation values are much more similar as expected, because high-elevation sites, beside higher clear-sky radiation, also have higher cloudiness due to frequent convective clouds. On the contrary, for SARAH-2.1 data these ratios show a rather different yearly cycle with respect to ground-based observations with up to 25% lower winter SSI values for high-elevation sites with respect to low-elevation ones. For CLARA-A2 data, some differences are evident too with a less pronounced yearly cycle than for ground-based observations and with a minimum falling in spring.
It will be interesting in the future, when the data become available, to compare the Piedmont ground-based series with SARAH-3 data records, the new version of SARAH-dataset that will also include a scheme to better identify and consider snow coverage.
In addition, the spatial representativity of the station measurements also needs to be considered when comparing local measurements with gridded datasets, which always represent a larger area (25 km 2 and 625 km 2 in the case of SARAH-2.1 and CLARA-A2, respectively) especially when areas with complex orography are considered [54]. However, the spatial resolution of the considered grid should be considered together with the temporal resolution of the data. In fact, for temporal resolution higher than the day, if the spatial resolution is of a few kilometres (e.g., 5 × 5 km) the representativeness error can be considered smaller. Differently, the representativeness error cannot be neglected when lower temporal or spatial scales are considered, even if the site is located at the centre of the grid [55].
Beside snow coverage, shading might also be a relevant factor explaining the differences in SSI values recorded by ground stations and corresponding satellite values, even though the effect seems to be rather limited for the 32 station sites (Section 2.5) for which we considered the SSI records as unbiased. This is evident for station site Borgofranco d'Ivrea, which has almost the same SSI when winter ground-based and SARAH-2.1 data are compared, whereas in summer SSI from SARAH-2.1 data is more than 20 Wm −2 higher. This station is in fact located in a valley surrounded by mountains both on the eastern and western sides, reducing the hours of sunshine during the day. Shading is indeed very difficult to eliminate in areas with complex orography and station sites which have good exposure may suffer from it at least in parts of the year. We have indeed checked the effect of shading from surrounding orography (see Section 2.5), however the effect of surrounding mountains can go beyond the shading they cause directly and may include shadings due to clouds that in Piedmont are more frequent over mountains than over surrounding plain areas.
The comparison of the ground-based and satellite anomaly records (Section 3.2) showed that for the low OCI stations the common variances between monthly satellite and ground-based records are about 84% and 83% for SARAH-2.1 and CLARA-A2, respectively. Moving from monthly to daily resolution, the common variances increase to 92.1% and 87.5% for SARAH-2.1 and CLARA-A2, respectively, giving evidence of an excellent performance of the daily satellite records in capturing SSI variability measured by ground-based stations. The better performance of the daily records with respect to the monthly ones could be due to the fact that the variability of daily anomalies is mainly due to meteorological factors, where the role of cloudiness is predominant and thus well captured by satellites, while the variability of monthly (and seasonal and yearly) records depends on other factors (e.g., aerosol concentrations) that are less accurately considered within the procedure used to retrieve SSI satellite data. However, the better performance of the daily records in terms of common variances depends also on the higher variability of the daily records with respect to the monthly, seasonal, and annual ones.
Considering the station sites characterised by high OCI, the common variance between satellite SSI anomalies and ground-based ones turned out to be much lower, as expected, due to the higher uncertainty both in satellite and ground-based records. Interestingly, the common variance between satellite SSI anomalies and ground-based ones turned out to be generally higher for CLARA-A2 than for SARAH-2.1 especially during winter.
Finally, trend analysis (Section 3.3) showed that the average SSI records from ground-based observations, at low-elevation station sites, gives evidence of a clear positive trend in the 1990-2016 period, whereas the average SSI records from satellite observations (both SARAH-2.1 and CLARA-A2) do not highlight significant positive trends. The average SSI record from ground-based observations generally has a positive trend also focusing on high-elevation station sites, even though the trend is significant only in summer and autumn and the increase observed for the annual record is lower than at low-elevation station sites. For high-elevation station sites, the differences in the trends of ground-based and satellite observations are less evident than for low-elevation station sites.
The different behaviour of ground-based and satellite observations supports the hypothesis of Manara et al. [32] on the relevant role of atmospheric aerosols on SSI trends over the investigated area. Satellite records are in fact derived by an approach that considers constant aerosol levels (see Sections 2.2 and 2.3) and therefore they do not take into account the significant reduction in aerosols in the period 1990-2016 [32]. A detailed quantification of this reduction by means of ground-based measures is not easy, because aerosol concentration data are available only for the last two decades. Moreover, the spatial coverage of aerosol measurements is highly inhomogeneous and very limited in high-elevation areas. Thus, we consider atmospheric optical depth (AOD) data in this discussion. Specifically, we consider all MERRA-2 grid cells (see Section 2.4) falling over the area considered in this study and we focus on the trend of the record obtained calculating the average yearly AODs for low-and high-elevation grid cells and considering the difference between these AODs. This record gives evidence of a significant negative trend (−0.006 ± 0.002 per decade-p-value ≤ 0.05) in the period 1990-2016 highlighting that the effect of the reduction in air pollutant emissions, which occurred over this period in Italy, affected atmospheric transparency more strongly at low-than at high-elevations. This effect is extensively discussed for the entire Italian territory by Manara et al. [34] in a study focusing on atmospheric visibility records.
It is therefore not surprising that satellite records do not give evidence of the same brightening signal of ground-based observations. Interestingly, the differences in the trends of SSI records from ground-based and satellite observations are more evident at low-elevation, whereas high elevation sites show a similar behaviour to low-elevation sites only in summer: in this season the boundary layer is highest and air pollution concerns a much deeper layer than in the other seasons [56][57][58].
The results obtained in this study are in agreement with those reported in literature even if they are not always easy to compare considering the different observed area, the different period covered, and sometimes the different version of satellite datasets. Specifically, Pfeifroth et al. [28] reported higher average values (bias equal to 2.7 Wm −2 ) for SARAH-2 than for ground-based stations as in this study, although they report lower average values for CLARA-A2 (bias equal to −0.4 Wm −2 ). However, it should be considered that they were calculated over the 1983-2015 period and for all of Europe. The reported trends are in agreement with stronger positive values for station data than for satellite data (both for SARAH-2.1 and CLARA-A2) in the southern part of Europe, especially during summer. Focusing only on Spain, Montero-Martin et al. [30] reported lower common variances (69% and 64% for SARAH-2.1 and CLARA-A2, respectively) than those reported in this study, even though it was always higher than 50% for both satellites, while for SSI trends they again reported lower positive trends for satellite data than for ground-based data especially in summer. Similar results are also reported by Alexandri et al. [59] for Greece, where they found an overestimation of about 7.1 Wm -2 for SARAH with respect to ground-based series.

Conclusions
The analyses presented in this paper aimed at comparing SSI CM SAF SARAH-2.1 and CLARA-A2 products with a dataset published detailing surface solar radiation from ground-based stations over the Piedmont region (north-western part of Italy) for the period 1990-2016 [32]. The goal was to provide a new contribution to the discussion on the representativeness of these satellite data, also including a focus on high elevation areas where satellite data are a key resource because the availability of station data is generally rather low.
The analyses focused onto: (i) the comparison of the average surface solar irradiance of the records from the different data sources; (ii) the comparison of the corresponding anomaly records that was evaluated both in terms of common variances and mean absolute differences; and (iii) the comparison of the trends of the records from the different data sources.
In spite of some open issues, CM-SAF SARAH-2.1 and CLARA-A2 surface solar irradiance products turned out to give evidence of an excellent ability in capturing SSI signals shown by ground-based stations. The three main results highlighted by the comparison between ground-based and satellite records turned out to be: (i) a tendency of satellite data-especially of SARAH-2.1-to be higher than corresponding ground-based observations for low-elevation areas; (ii) a tendency of satellite data to be lower than corresponding ground-based observations for high-elevation areas mainly because of snow coverage not always correctly distinguished from clouds; and iii) the inability of satellite observations to correctly capture the brightening signal highlighted by ground-based observations. This last result, already reported by several papers focusing on other southern areas of Europe, is probably caused by relevant changes in atmospheric aerosol levels that SARAH-2.1 and CLARA-A2 surface solar irradiance records are not able to capture because they are based on constant aerosol climatologies. Other factors, also related to station data issues, cannot, however, completely be ruled out in explaining the different trends of ground-based and satellite observations. Problems in station data, especially in areas with complex orography, can naturally also not be ruled out in explaining the other differences between ground-based and satellite records highlighted by our analyses. The analyses presented in this study have in fact highlighted that the comparison of surface solar irradiance records from different data sources is a complex issue that has to be performed, critically evaluating the strengths and the weaknesses of each of them in relation to the objectives of the analysis. For example, a study has shown that short-and long-wave satellite data can be successfully used on a glacier surface to estimate summer glacier melting if they are not available from on-glacier surface measurements [60], even if satellite data at high-elevations present some limitations.