Improved estimation of global solar radiation over rugged terrains by the disaggregation of Satellite Applications Facility on Land Surface Analysis data (LSA SAF)

This paper presents a new method to predict global solar radiation over irregular terrain, named Estimation of global solar RADiation (ERAD). The method is based on the disaggregation of Satellite Applications Facility on Land Surface Analysis (LSA SAF) data using a digital elevation model and is applied in Italy with a time step of 1 min and a spatial resolution of 200 m. A quantitative assessment of ERAD is performed in comparison with three other standard methods (Mountain Microclimate Simulation Model [MTCLIM], LSA SAF and Copernicus Atmosphere Monitoring Service [CAMS]) using measurements taken in 43 stations located in Italy or in the surrounding countries, in the years 2005–2016. Such assessment concerns the irradiance incoming on a horizontal surface, which is measured by ground radiation sensors and is summarized by means of four accuracy statistics (i.e. mean absolute error [MAE], root mean square error [RMSE], coefficient of determination [R2] and mean bias error [MBE]). Overall, the average daily global solar radiation estimates obtained by ERAD have RMSE and R2 about 25 W·m−2 and 0.943, respectively. These statistics are similar to those of LSA SAF and better than those of CAMS and, above all, MTCLIM. The bias analysis by elevation ranges shows a slight ERAD overestimation over plains and hills and a slight underestimation over mountains. An additional qualitative assessment shows how the ERAD radiation estimates are more spatially detailed than those of the other methods and are redistributed on inclined surfaces consistently with expectations.


| INTRODUCTION
Downwelling Surface Shortwave radiation Flux (DSSF) is the solar irradiance reaching the Earth's surface per unit of time and area (Geiger et al., 2008a). The solar flux, consisting predominantly of ultraviolet-visible and nearinfrared (300-4,000 nm) radiation, affects the weather, the climate and most biological processes in terrestrial and marine ecosystems. Therefore, solar radiation is a crucial factor in global warming mechanisms, glacier retreat, water resource and carbon budgeting (Wild, 2016). In addition, the prediction of solar radiation is decisive to estimate evapotranspiration (Aguilar et al., 2010), to predict crop and forest yield  and for numerous solar energy applications such as photovoltaic systems (Demain et al., 2013).
Outside the atmosphere, the solar flux is almost constant and is referred to as solar constant (I 0 ). The solar flux, before reaching the Earth's surface, interacts with atmospheric constituents and is partially reduced by the phenomena of scattering, absorption and reflection. In clear sky conditions, the atmospheric transmittance mainly depends on water vapour and aerosol characteristics (Lefèvre et al., 2013). Instead, in cloudy or partially cloudy sky conditions, the prevailing factors are the optical and physical properties of clouds (Oumbe et al., 2009). At ground level, the incoming solar radiation depends on various astronomical and geographical factors including the terrain characteristics (elevation, slope and aspect), which are usually derived from a digital elevation model (DEM). In morphologically complex areas, the grid resolution of the DEM is a factor to consider with due attention, since it affects the estimation accuracy of solar radiation (Huang and Zhao, 2017).
The global solar radiation (I g ) that reaches the surface can be split into direct or beam (I b ), diffuse (I d ) and reflected (I r ) components (Iqbal, 1983) and can be measured by specific sensors (e.g. pyranometers) at weather stations. However, the stations that measure solar radiation are few because of the high cost of sensor installation and maintenance. One of the most famous solar radiation station networks is the World Radiation Data Centre (WRDC: http://wrdc.mgo.rssi.ru) -World Meteorological Organization (WMO) consisting of about 1,600 stations distributed all over the globe. These stations often collect solar radiation in the form of daily or monthly averages but usually show low accuracy and reliability. An exception is represented by the World Radiation Monitoring Center (WRMC) -Baseline Surface Radiation Network (BSRN) (Driemel et al., 2018), although there are only 64 and 11 stations over the globe and in the European area, respectively. To overcome the problem of the absence of stations, many indirect methods have been developed based on the correlation between solar radiation and other more easily measurable meteorological or geographical parameters, such as sunshine duration (Suehrcke et al., 2013). A typical indirect model is the Mountain Microclimate Simulation Model (MTCLIM) proposed by Thornton and Running (1999) which essentially uses daily temperature range and rainfall observations to estimate incoming solar radiation.
Among the most common methods for estimating solar radiation are those based on meteorological satellites, such as the Meteosat Second Generation (MSG). Geostationary satellites allow the diurnal cycle of solar radiation, or of cloudiness, to be captured due to the high temporal scanning frequency of the acquired images (15 min for MSG) (Journée and Bertrand, 2010). Conversely, the drawbacks of these satellite data are due to low spatial resolution (5 km in the Italian area) and to the consequent poor accuracy in rugged areas.
Classically, the methods to retrieve DSSF from satellite data can follow three main schemes: empirical, physical and semi-empirical (Sengupta et al., 2015). Empirical models are based on statistical relationships between the information derived simultaneously from satellite data and ground observations (Tarpley, 1979). Physical approaches use satellite data to solve radiative transfer models (RTMs). These last methods are computationally more expensive and require a detailed and always updated knowledge of the optical properties of the atmosphere (Urraca et al., 2017). A series of methods, adopted by many researchers and referred to as Heliosat methods, have evolved towards a semi-empirical, or hybrid, approach. From the first version of the Heliosat method (Cano et al., 1986) various improvements were proposed in successive versions, characterized by increasing numbers (Heliosat-1, Beyer et al., 1996;Heliosat-2, Rigollier et al., 2004;Heliosat-3, Schroedter-Homscheidt et al., 2006). All these methods use physical models to describe, in cloudless conditions, the interaction of solar radiation with gases and aerosols in the atmosphere, while satellitederived data are used to retrieve the cloudiness conditions empirically. Another semi-empirical approach is implemented by the European Organization for the Exploitation of Meteorological Satellites (EUMETSAT) -Satellite Applications Facility on Land Surface Analysis (LSA SAF; Trigo et al., 2011). The LSA SAF service provides, amongst others, DSSF products derived from MSG data. In this case, the empirical observations are related both to the presence and optical characteristics of the clouds and to the land surface albedo (Geiger et al., 2008a;2008b).
In the last few years, many semi-empirical methods tend to include more physical descriptions. The convergence to physical methods is feasible due to both the improvements obtained in the RTM calculation times by using a series of discrete pre-calculated solutions saved in look-up tables and the higher availability of ancillary information related to the physical properties of the atmosphere (Urraca et al., 2017). An example of these methods is the new version of Heliosat, referred to as Heliosat-4, based on a library of RTMs and distributed by the Copernicus Atmosphere Monitoring Service (CAMS) (Qu et al., 2017).
As the spatial and temporal distribution of surface solar radiation depends on topographic effects and shadows due to the horizon, spatial disaggregation methods have been demonstrated to be efficient to improve the spatial detail of existing, low resolution satellite datasets (Ruiz-Arias et al., 2010). In particular, such methods are aimed at adding information related to the effects of topography, thus enhancing the accuracy of solar radiation estimates over irregular terrains. Following this research line, a procedure for predicting solar radiation over rugged terrains, named Estimation of global solar RADiation (ERAD), is currently put forward. This semi-empirical disaggregation method, which is similar to others proposed in the literature, performs a downscaling based on a DEM (Ruiz-Arias et al., 2010;Haurant et al., 2012;Moreno et al., 2013). In particular, ERAD can improve the DSSF LSA SAF product, which has spatial and temporal resolutions of 5 km and 30 min, respectively, producing solar radiation estimates with a spatial resolution of 200 m and a time step of 1 min. The new disaggregation method is applied over a wide study area which includes Italy, and the estimates obtained are validated against a consistent dataset of daily ground observations taken from 2005 to 2016. An intercomparison is also performed versus radiation estimates obtained from three other standard methods, MTCLIM, LSA SAF and CAMS. Finally, the improvements produced by the new method over the other algorithms are highlighted by examining some examples of how the respective solar radiation estimates vary in space and time.

| RADIATION ESTIMATION METHODS
The next sections briefly introduce the three standard methods used for intercomparison, followed by a detailed description of the new ERAD algorithm.

| MTCLIM
MTCLIM is an indirect method for the simulation of daily solar radiation based on meteorological data usually collected in weather stations. The method uses the minimum and maximum daily temperature together with the daily precipitation and other geographical information such as latitude, elevation, slope and aspect (Thornton et al., 2000). The application of MTCLIM over large regions with complex morphologies requires the preliminary extrapolation-interpolation of meteorological variables, such as those obtained by the DAYMET algorithm (https://daymet.ornl.gov) with 1 km resolution (Thornton et al., 1997) or the European Observation database (https://www.ecad.eu), with 0.25 resolution (Haylock et al., 2008). MTCLIM has a universal valence, being usable over a wide range of climates without reparameterization on a site-by-site basis.
Indeed, when applied to BSRN station data, the method produced reasonable results over most sites and under most geographical and climate conditions (Bohn et al., 2013). The version of MTCLIM used in the present study is 4.3, although the empirical coefficients related to daily temperature range are those defined in the original version of the method (Thornton and Running, 1999).

| LSA SAF
The current MSG satellites are supported by EUMETSAT in collaboration with the European Space Agency. Additionally, EUMETSAT supervises and coordinates the activities of LSA SAF for routinely providing some solar radiation products, such as MSG DSSF (LSA-201), available at the website http://landsaf.ipma.pt. The DSSF products, obtained by the approach described by Geiger et al. (2008a), have the same projection and spatial resolution (3 km at nadir) as the original MSG imagery. Overall, DSSF mainly depends on the solar constant, the cosine of the solar zenith angle, the transmittance of the atmosphere (which is determined by cloud cover), the atmospheric absorption and scattering and the surface albedo. For the SAF products, the distribution of some atmospheric gases (water vapour and ozone) is derived from numerical weather forecasts of the European Centre for Medium-Range Weather Forecasts and the Total Ozone Mapping Spectrometer, respectively. The land surface albedo is derived from MSG, too (Geiger et al., 2008a). In summary, the DSSF products are obtained by calculating the effective transmittance of the atmosphere, both in clear and cloudy sky conditions, and multiplying it for the solar constant. In the clear sky case, MSG data are not used and the transmittance is determined from a physical model of the atmosphere with visibility fixed to 20 km and a continental-type aerosol. Instead, in the cloudy sky case, a simple RTM is employed and the MSG observations are a fundamental input. In this case, in fact, more radiation is reflected by the clouds (brighter pixels in MSG images) and less radiation reaches the ground. However, some approximations are made, such as the one for which the whole image pixels are covered by a homogeneous cloud layer (Geiger et al., 2008a). Therefore, the cloud transmittance of each pixel is related to cloud albedo by an empirical cloud absorption factor. The validation of the DSSF LSA SAF product yielded good results in different study regions (Journée and Bertrand, 2010;Cristóbal and Anderson, 2013

| CAMS
CAMS uses the Heliosat-4 method based on MSG imagery to produce solar radiation estimates, mostly at the site level, covering European and African regions since 2004. The Heliosat-4 method is a physical approach based on two libraries of RTMs that use look-up tables: the McClear and McCloud models (Qu et al., 2017). One of the advantages of this method is that it allows the radiation received at ground level in all sky conditions to be estimated. The McClear model provides the clear sky irradiance, which is multiplied by a modification factor, while the McCloud model takes into account the effects of radiation extinction due to clouds. The main input variables (aerosol optical properties, water vapour and ozone) are updated every 3hr by CAMS on a global scale. Ground albedo is determined using Moderate Resolution Imaging Spectroradiometer (MODIS) imagery. The McCloud model uses the physical optical properties of the clouds (e.g. optical depth, type and coverage of cloud) derived from MSG observations, every 15 min and at 5 km of resolution, to calculate the cloudy sky irradiance. The Heliosat-4 model generally exhibits satisfactory performance (Marchand et al., 2019) although some limitations are known (Qu et al., 2017). The CAMS products are related to the solar radiation incoming on a horizontal surface of a certain site, for the actual weather conditions, at different time steps (from 1 min to 1 month). If not included in the user request, the CAMS service provides the altitude of the site on the basis of some DEMs available, including the Shuttle Radar Topography Mission data with 100 m of spatial resolution (CAMS, 2019). The current solar radiation product is downloaded through the website http://www.soda-pro.com.

| ERAD
ERAD is a semi-empirical method based upon a disaggregation procedure applied to the DSSF LSA SAF product and a DEM of the study area. This spatial disaggregation method consists in the redistribution of the estimated solar radiation evaluating carefully all topographic effects. The resolution of the used DEM must obviously be much higher than that of the original satellite data (5 km), so that the estimated solar radiation offers highly improved spatial detail over rugged terrains (Moreno et al., 2013). ERAD provides the solar radiation for each point of the DEM and 1 min time steps, similarly to the procedure of Haurant et al. (2012). As LSA SAF provides half-hourly solar radiation at a resolution of about 5 km, the disaggregation method concerns both spatial and temporal scales. The calculation of the daily solar radiation average integrates all estimates of the daytime, which requires some approximations. The most significant of these is the assumption of a homogeneous spatial distribution and a constant temporal variation of all variables considered in the radiation estimation process for each DEM point. A classical linear interpolation method is applied to the variables sampled at different time intervals compared to the basic 1 min step. The application of this method also requires the calculation of numerous astronomical features, such as the true and apparent solar position in the sky, or other geographical features such as the sky view factor (SVF), for each point of the DEM. All these features are defined on the basis of standard reference guidelines (e.g. Iqbal, 1983;Flint and Childs, 1987;Michalsky, 1988;Sengupta et al., 2015). ERAD determines the solar radiation in two phases as described in the following sections.

| Half-hourly data processing
ERAD has as main input the half-hourly DSSF images provided by LSA SAF (LSA-201) product. An inversion of the physical model described by Geiger et al. (2008a) is applied to estimate the actual transmittance of the atmosphere by the clearness index (k t ) defined as the ratio (Skartveit et al., 1998): where I gh and I 0h are the global and extraterrestrial global radiation, estimated on a horizontal surface. In ERAD the extraterrestrial global radiation (I 0h ) is calculated as (Iqbal, 1983): where I 0 is the solar constant over the whole solar spectrum (1,367 WÁm −2 ) and the expression in square brackets corresponds to the Earth's orbit eccentricity factor, in which DOY is the day of the year. Instead, LSA SAF estimates k t by integrating the radiation in a limited range of the solar spectrum, 300-4,000 nm, corresponding to a solar constant of 1,358 WÁm −2 (Geiger et al., 2008a).
The cosine of apparent solar zenith angle (θ SZA ) can be expressed as a function of some astronomical parameters such as the DOY declination (δ), the latitude (Φ) and the hourly angle (ω) calculated at the centre point of each DEM pixel, in accordance with the acquisition times of the MSG overpasses: where Φ is known and δ and ω are calculated following Michalsky (1988). The clearness index (k t ) is a dimensionless number between 0 and 1, and has a high value under clear sky and a low value under cloudy conditions. In some circumstances, such as in cloud enhancement situations or high albedo due to snowy conditions, k t could exceed unity. A maximum threshold of k t , equal to 0.9, is therefore imposed to avoid these situations. Moreover, k t is calculated only when the solar zenith angle is less than 85 to avoid miscalculation problems when the sun is near the horizon.
At regional or local scale, the distribution of solar radiation is affected by elevation (z) which represents a topographic effect. In fact, due also to the reduction of the optical path length, an increase in altitude generally produces a rise of the global radiation incoming on the surface. The elevation correction is made for each point of the image considered. As the solar global radiation value I g (z 0 ) derived from the DSSF product is referred to sea level (z 0 ) (Moreno et al., 2013), a z-correction is applied according to Wahab et al. (2010). Hence, given the global radiation at z 0 , the global radiation at elevation z is obtained from the following empirical function: where τ(z 0 ) = -ln[I g (z 0 )/I 0 ] and β = 1.2. Consequently, for each elevation z of DEM point, a corrected clearness index k t is defined as: where k t (z 0 ) = I g (z 0 )/I 0 and z 0 = 0. This obviously implies that every DEM point is associated with the nearest point on the k t DSSF image. As the elevation effect on the global radiation is mainly due to the beam component of solar radiation, the elevation correction is applied only in prevailing clear sky conditions. These conditions are met when the correct clearness index k t exceeds a threshold value of 0.65, in accordance with Ruiz-Arias et al. (2010).
After the elevation correction, ERAD uses a decomposition model to subdivide the global solar radiation received on a horizontal plane into beam and diffuse radiation components. Among the numerous decomposition models proposed in the literature (Bertrand et al., 2015), the model of Skartveit et al. (1998) is used. This model is among the most widespread and has recently been found to be one of the most efficient in a comparison with 13 other models for a highly rugged area near Bolzano (northern Italy) (Laiti et al., 2018). The model is based on three main input predictors: the solar zenith angle, the clearness index and a temporal variability index (hourly in the original model and semi-hourly in ERAD).
In this decomposition model, the clearness index depends also on the local land surface albedo, which is derived from the daily albedo product provided by LSA SAF dataset. Therefore, in the ERAD procedure, the Skartveit model allows the diffuse component of the solar radiation to be derived by the cloudiness index (k d ), defined as the ratio between diffuse (I dh ) and global (I gh ) solar radiation that reaches a horizontal surface: The cloudiness index (k d ) is obtained with the same temporal resolution of DSSF product (30 min) and therefore k d is assumed to remain constant in this time interval.

| One minute data processing
Since the position of the sun in the sky changes continuously throughout the day, the simulation of the other topographic effects on radiation distribution is performed using a 1 min time step. For this step, the global (I gh ), diffuse (I dh ) and beam (I bh ) solar radiation incoming on a horizontal surface are calculated for each DEM point on the basis of the previously determined parameters (k t (z), k d ) by means of the following equations: Transposition models are then used to convert the solar radiation on the horizontal surface to that on the tilted surface (Padovan and Del Col, 2010) as determined by the DEM. These models define how the various components of solar radiation arrive on the terrain from the sky and the surrounding land surface. SVF is a key factor to reach this objective since it takes into account the light conditions in relation to the DEM.
The first condition calculated for each time step (1 min) of the day is if the site (or point) considered is directly illuminated by the Sun or if it is in shade. Given the complexity of the actual lighting conditions, it is convenient to treat each of the solar radiation components separately and, lastly, sum the beam, diffuse and reflected components to obtain the daily global solar radiation incoming on the tilted surface.
The beam radiation component, on a tilted (T) surface (I bT ), follows a purely geometric relation with respect to the horizontal component (I bh ) (Iqbal, 1983). In fact, it depends on the following ratio: where θ is the incidence angle of the Sun's light rays compared to the normal at the tilted surface. cos(θ) of Equation (10) can be obtained using the following expression: where S and A are the slope and aspect, respectively, and ψ is the solar azimuth which can be derived from the formula: and α is the apparent solar elevation. Hence, the beam radiation component (I bT ) can finally be computed as: The s b factor in Equation (13), which will be described later, takes into account the Sun's position to determine the partial, or total, shading conditions for each DEM point.
For the sky diffuse radiation component the type of solar radiation distribution needs to be defined. For example, if an isotropic distribution model is assumed, the intensity of radiation is supposed to be uniform over the sky dome. The diffuse radiation incident on a tilted plane, with slope S, is usually approximated by multiplying the diffuse radiation on a horizontal surface (I dh ) by transposition factor or view factor to the sky of the diffuse component, which is given by [1 + cos(S)]/2 (Demain et al., 2013). ERAD instead uses the more accurate Hay anisotropic model, for which the diffuse radiation (I dh ) is divided into two components (Iqbal, 1983): an isotropic component (I di ), which implies a uniform radiation distribution from the sky, and a circumsolar anisotropic component (I dc ). The latter is related to the incoming radiation from the area surrounding the solar disc, which can be analysed as beam irradiation. As suggested from the Hay model, the two diffuse components are weighted according to an isotropy index, defined from the ratio I bh /I 0h . Consequently, the isotropic radiation can be evaluated by I di = I dh (1 -I bh /I 0h ) while the circumsolar radiation is evaluated by I dc = I dh I bh /I 0h . Moreover, the diffuse radiation received by a tilted surface also depends on some topographic factors. In fact, taking into account the SVF for each point of the DEM, the isotropic diffuse radiation will be limited by the horizon, and therefore from the SVF value (0 < SVF < 1), except in the case of a wide flat area for which SVF = 1. Hence, similarly to Aguilar et al. (2010), given the diffuse solar radiation on a horizontal surface (I dh ), the diffuse and circumsolar radiation components on a tilted surface, I diT and I dcT respectively, are obtained from the following equations: where SVF is a function of slope (S), aspect (A) and local horizon angle (α hor,i ), in each of N d directions i, equally distributed on the plane selected starting from the centre of each DEM point. More specifically, SVF is obtained by integrating in all directions i, with azimuth i(360/N d ), the following expression, equivalent to that suggested by Dozier and Frew (1990): where N d is the number of directions i (72 in the present study), the summation is evaluated for each direction i, and γ is the difference between the azimuth of each direction i and the aspect (A) of the surface considered. The s b and s dc factors in Equations (13) and (15) define the shadowing conditions for the beam and diffuse-circumsolar components, respectively. In fact, the shading from the Sun depends on the Sun's position referred to the horizon. This phenomenon is usually taken into account by a discrete Heaviside step function, whose values are 0, 0.5 or 1 when the Sun is below, halfway or above the horizon, respectively (Ruiz-Arias et al., 2010). A similar function is currently used for the s b and s dc factors, which differ from the classical Heaviside function when the Sun is partially obscured from the horizon. In this case the value of s b , or s dc , is not 0.5 but is variable between 0 and 1, equal to the ratio between the area of the solar segment which is visible and the total area of the solar disc. The circumsolar radiation instead corresponds to a region around the solar disc (solar aureole). In accordance with the measurements of Flint and Childs (1987), the circumsolar region, viewed from each DEM point, is defined by a cone of sunlight with an opening angle of 4 ; the same cone, relative to the solar disc alone, is on average 0.52 .
The last component of global solar radiation incoming on a tilted surface is the ground-reflected radiation. As ground albedo and topography can vary also over a short distance, their interaction can lead to a high variability of the radiation field, in terms of both intensity and spatial-temporal distribution. Since the surface reflection characteristics are not easy to find, the reflected radiation component is usually assumed to be isotropic. In this case, each surface point reflects a constant radiance, as a perfectly diffuse reflector, which depends on the view factor to ground: (1 − cos(S))/2 (Demain et al., 2013). In the present study, knowing the global solar radiation on a horizontal surface (I gh ), the reflected radiation on a tilted surface (I rT ) coming from the surrounding terrain is estimated using the following equation (Aguilar et al., 2010): where ρ is the regional albedo derived from the LSA SAF dataset. The term in square brackets in Equation (17) represents the terrain transposition factor, or view factor to the ground of the reflected component, named ground view factor (GVF). In the case of a horizontal surface (S = 0), the GVF is complementary to SVF, so that GVF = 1 − SVF. The estimation of the ground albedo is fundamental for the assessment of reflected radiation, but the assumptions of isotropy for this radiation and of daily constant distribution for ρ are both unrealistic. In particular, ρ shows a high rate of daily and seasonal variability due to various factors, such as changes in soil water content, vegetation status, snow cover (Demain et al., 2013). The global solar radiation on a tilted surface corrected for all topographic effects (I gT ) is therefore obtained for each minute as the sum of all solar components derived from Equations (13)- (17): The global solar radiation on a horizontal surface is obviously obtained using S = 0 in each of the above defined components. In the case of daily outputs, the

| Study area and data
The four radiation estimation methods were assessed in an area which corresponds to the land surface of Italy and its surroundings (Figure 1a,b). This area is characterized by a wide variety of morphological, vegetation cover and climatic features. Due to the presence of the Alps chain in the north and the Apennines chain along the ridge of the peninsula, the Italian landscape is prevailingly mountainous and hilly. Most of the country benefits from a typical Mediterranean climate characterized by warm dry summers and mild wet winters. Figure 1a shows the DEM of the study area with the position of the 43 weather stations currently used. The DEM lies between 35.05 and 47.18 N and 6.50 and 19.71 E and is derived from Shuttle Radar Topography Mission data (SRTM-3, http://srtm.csi.cgiar.org) after a resampling to 7.5 arc second (about 200 m in the study area) from the original resolution of 3 arc second. The DEM area extends outside the Italian State and includes various weather stations in the surrounding countries. In addition to elevation, other morphological parameters such as slope (S) and aspect (A) are directly derived from the DEM by standard algorithms. The DEM is also used for calculation of the SVF, which is defined as the fraction of sky visible above a certain ground observation point (Dozier and Frew, 1990) (Figure 1b). Table 1 indicates the elevations, the geographical coordinates and the periods of data collection of the weather stations used in the study. The stations are ordered with increasing elevation and are divided into three categories (plain, hill and mountain) if their altitude is less than 200, between 200 and 800 and higher than 800 m above sea level, respectively. Apart from some stations of the LaMMA Consortium, Sesto Fiorentino, Livorno and Grosseto, most of these stations belong to the WRDC-WMO network; two exceptions are Payerne and Sonnblick, which belong to the BSRN.

| Performance evaluation
The performances of the estimation methods considered (MTCLIM, LSA SAF, CAMS and ERAD) were evaluated both quantitatively and qualitatively. The first assessment was carried out versus daily global solar radiation observations taken at the 43 reference ground stations over the period 2005-2016. Since all these observations are referred to horizontal surfaces, their use for the current purpose had to rely on some assumptions which are discussed in Section 5.
The four methods used different drivers (input data). Both MTCLIM and ERAD used the available DEM. Additionally, MTCLIM required daily temperature and rainfall measurements, while the other methods used half-hourly MSG imagery. Each estimated data series therefore showed a different number of missing values, which were due to various causes (data processing system failures, the application of filtering criteria, etc.). For example, the average daily global solar radiation could not be computed when some diurnal cycles of the MSG images were missing and therefore not representative of all daytime. The data series estimated by MTCLIM, LSA SAF, CAMS and ERAD for all ground stations consisted of 94,635, 59,240, 94,393 and 95,069 daily samples, respectively. The MTCLIM, CAMS and ERAD datasets had about 48% missing data in the period examined, while LSA SAF dataset had about 37% additional missing data. Finally, the common dataset of measured and estimated daily solar radiation data consisted of 57,534 daily samples, which were used for the accuracy assessment.
The accuracy of each method was summarized by means of common statistics, that is, the mean bias error (MBE), the mean absolute error (MAE) and the root mean square error (RMSE), computed using the following equations: where e i and o i are the estimated and observed daily values, respectively, for every day sampled (i) on all data series (N) considered. To obtain dimensionless errors, these statistics were also transformed into percentages (%), dividing each value by the average of the ground observations and multiplying by 100. The accuracy assessment was completed by calculating the coefficient of determination (R 2 ) between measurements and estimates. All these statistics were computed both considering all ground stations and these aggregated by elevation (plain, hill and mountain). The performances of the radiation methods considered were also qualitatively evaluated by visually analysing the spatial and temporal variations of the respective solar radiation estimates. In the case of MTCLIM, high spatial resolution daily radiation maps could be obtained only for a region of Italy, Tuscany, where the DAYMET interpolation algorithm had been applied to a great number of ground meteorological observations and the same DEM used in the present study. For LSA SAF the visual analysis was obviously limited only to the comparison between the estimates of solar radiation on a horizontal surface; the average annual map was produced by integrating all available daily data for which no SAF product (semi-hourly) during the diurnal cycle was missing. CAMS limits the number of daily requests and prevents the acquisition of solar radiation estimates at area level. Actually, CAMS provides also a gridded product, whose spatial resolution (0.2 ), however, is incompatible with that of the DEM currently used. Thus, only ERAD could be applied to the entire Italian territory to produce high spatial resolution (200 m) solar radiation estimates over rugged terrains taking into account all topographic effects; this exercise was carried out for an exemplary year (2015), yielding 30 min images which were finally aggregated to an annual map.

| Quantitative evaluation of estimated radiation
All 43 ground stations currently considered are placed in areas where solar radiation is not affected by the surroundings reliefs; for all stations, in fact, the SVF is greater than 0.94. Tables 2 and 3 summarize the results of the statistical evaluation of the solar radiation estimates obtained by the four models at these stations. More particularly, Table 2 shows the numbers of observations and the daily global solar radiation averages of the measured data series for all stations and for these aggregated 16.8 (9.9) 19.3 (11.4) 16.7 (9.9) Plain 32.9 (18.6) 11.6 (6.6) 13.4 (7.6) 12.9 (7.3) Hill 30.4 (19.7) 11.5 (7.5) 14.6 (9.5) 12.  by elevation (plain, hill and mountain), while Table 3 shows the accuracy statistics obtained by each method referred to both the original and percentage values.
These results are confirmed by Figure 2, which displays the scatterplots of daily global solar radiation measured and estimated by the four models for all stations. The estimates of MTCLIM show the greatest dispersion around the 1:1 line, which testifies to the worst performance of this model. The results are better for LSA SAF and CAMS, but the former model shows a clear underestimation for the highest radiation values. Overall, ERAD has the smallest dispersion and no clear over/underestimation pattern.
The results are slightly different when considering the observations divided into the three elevation ranges (Table 3); the interpretation of these results, however, is complicated by the incomplete concordance of the four accuracy statistics, particularly the MBE.

| Qualitative evaluation of estimated radiation
As previously mentioned, the qualitative assessment of the solar radiation products concerned the Italian territory for the year 2015. Figure 3a surface. The estimates of the two methods have a similar spatial distribution, but the former model predicts lower radiation values particularly over mountain areas. Figure 4a,b shows the maps of mean annual daily global solar radiation estimated by MTCLIM (a) and ERAD (b), respectively, in a region of complex topography (Tuscany). MTCLIM produces generally reasonable radiation estimates, but with some anomalies which are probably associated with areas having a low density of the meteorological station network in 2015. This defect is not visible in the map obtained by ERAD, which simulates solar radiation as expected considering all described topographic effects.
Both Figures 3a,b and 4a,b confirm the expected trends of solar radiation related to latitude and elevation. There is an increase of mean solar radiation for lower latitudes due to different solar elevation angles and a decrease of radiation for mountain areas due to increased cloudiness and, in narrow and deep valleys, to the shading by mountain slopes.   The performances of MTCLIM and ERAD are further analysed by examining the annual radiation evolutions estimated for three pairs of sites situated over northsouth slopes in a Tuscany mountain area (Figure 4b,c). Table 4 shows the main morphological features of these sites and the daily global solar radiation averages estimated by the four models for 2015.
LSA SAF and CAMS are both incapable of differentiating the solar radiation incident on different slopes. In contrast, both MTCLIM and ERAD can reasonably reproduce the higher solar radiation expected for southern slopes. The latter model, however, predicts higher radiation differences between contrasting slopes, that is, a higher effect of topography. Figure 5a,b more clearly illustrates the annual trends predicted by the last two models for the second site pair. Both models simulate a prevalence of direct radiation incoming on the southern slope, with consequently higher day-to-day variability with respect to the northern slope, where diffuse radiation prevails. In accordance with what is observed above, these features are more pronounced for ERAD, that is, this model shows a higher sensitivity to the contrasting topographic characteristics of the two sites ( Figure 5b).

| DISCUSSION
The current study has put forward a new method to predict global solar radiation over rugged terrains, ERAD. The method performs a disaggregation of LSA SAF data using a DEM with a spatial resolution of 200 m. Being based on LSA SAF data, ERAD presents most of the approximations brought by this product but provides more spatially and temporally detailed radiation estimates related to the resolution of the DEM used and the time step applied for data processing.
The results of ERAD have been inter-compared with those of three other standard models (MTCLIM, LSA SAF and CAMS) using daily solar global radiation observations taken in 43 stations in or around the Italian peninsula during the years 2005-2016. As is the case for almost all other solar radiation measurements taken worldwide, these observations are not optimal for the current purpose, being referred to horizontal surfaces. The model validation had therefore to rely on the assumption that an accurate estimation of the main radiation components (direct and diffuse) in this situation is a fundamental pre-requisite for a similar prediction over tilted surfaces (Ruiz-Arias et al., 2011). According to the theory presented in Section 2.4.1, the estimation of these components is intrinsically associated with that of global solar radiation, which has actually been assessed versus ground observations. The validity of this indirect evaluation strategy, which has been used by other authors (see for example Moreno et al., 2013), is obviously dependent on the soundness of the physical theory applied for redistributing the estimated radiation components in all sky-terrain conditions. Particular attention was therefore devoted to using a straightforward theoretical background capable of taking into account the impact of numerous interacting astronomical, geographical, topographic and atmospheric factors on radiation simulation.
Based on these premises, the performances of ERAD and the other models have been quantitatively analysed by the use of standard accuracy statistics. Overall, the improvement produced by ERAD is particularly evident with respect to MTCLIM and marginal with respect to LSA SAF and CAMS. Considering all 43 ground stations, ERAD yields the best accuracy, with a MAE of 16.7 WÁm −2 (9.9%), RMSE of 25.3 WÁm −2 (14.9%), R 2 of 0.943 and MBE of 2.0 WÁm −2 (1.2%). This analysis was completed by an assessment in three altitude belts, which indicates some peculiarities in the performance of the four methods. MTCLIM shows a clear radiation underestimation in plain areas, which probably derives from the presence of some coastal sites where the MBEs are strongly negative. For example, considering only the coastal stations of Pantelleria, Livorno and Nice, the MBEs are −35.2, −34.3 and −29.7 WÁm −2 respectively. As noted by Bohn et al. (2013), these patterns probably depend on the low daily temperature ranges usually found in coastal areas, which affect MTCLIM simulation of diffuse radiation.
LSA SAF shows a general tendency to underestimate solar radiation, which is most evident over mountains, while ERAD shows a marginal overestimation over plains and hills, which is reversed at higher altitudes. These differences can be attributed to the use of slightly different radiation constants by LSA SAF and ERAD (1,358 and 1,367 WÁm −2 , respectively) and to ERAD accounting for the dependence of atmospheric transmittance on altitude, which is neglected by LSA SAF. The latter factor is obviously most influential over mountains, and this explains the almost complete correction of the radiation underestimation caused by LSA SAF in these cases. These results are similar to those obtained by Moreno et al. (2013) in a study conducted over Peninsular Spain; also these authors, in fact, found a clear radiation underestimation brought by LSA SAF in mountain areas, which was fixed by the application of an elevation correction.
Next, a qualitative evaluation was carried out using exemplary annual series of solar radiation maps produced by three of the models (MTCLIM, LSA SAF and ERAD). This analysis confirms the expected latitude and altitude trends of the solar radiation estimates, which originate from the astronomical and topographic factors mentioned previously. The same analysis confirms that ERAD, as well as MTCLIM, can take into account the effect of surface tilting, which is not considered by LSA SAF and CAMS. The performances of MTCLIM and ERAD over rugged terrains have therefore been analysed by examining the radiation evolution on some site pairs having contrasting north-south aspects. Such analysis shows that ERAD is more sensitive than MTCLIM to topographic effects, that is, more clearly differentiates the solar radiation incoming on slopes with different aspects. In particular, the greater differentiation of north and south slopes predicted by ERAD is due to higher fractions of direct radiation estimated on clear days, which also leads to greater radiation drops on cloudy days and, consequently, to higher day-to-day radiation variability.
While a direct assessment of these different radiation trends is prevented by the mentioned lack of measurements on tilted surfaces, it can be noted that the limits of MTCLIM are widely recognized in the literature and have been confirmed by the current quantitative analyses. These limits are mostly related to the operational nature of this method, which uses standard daily meteorological observations (i.e. temperature range and rainfall) to model atmospheric conditions and consequently estimate direct and diffuse radiation fractions. These considerations lead reasonably to presume a higher accuracy of the radiation evolution predicted by ERAD over rugged terrains.
Such expected better performance is also related to the use of a 1 min time step to calculate all model variables, which corresponds to a temporal frequency much higher than that of the original LSA SAF product (half an hour). This choice should particularly affect the characterization of some topographic factors, such as the SVF, theoretically leading to a higher accuracy of the solar radiation estimates. The actual impact of this choice on the radiation estimation accuracy, however, has not been currently ascertained and could be the subject of future investigations.
In general, our analyses have confirmed the critical importance of the main DEM characteristics for disaggregating solar radiation over rugged terrains. It must in fact be recalled that the proper definition of surface slopes and aspects is dependent on the spatial resolution and accuracy of the DEM used (Kang-tsung and Bor-wen, 1991). The current choice of a medium resolution DEM (200 m) is dictated by both theoretical constraints and operational reasons. The low spatial resolution of the original satellite product (5 km), in fact, poses limits to the prediction of atmospheric variability in much smaller areas (CAMS, 2019). At the same time, the resolution chosen allows the straightforward production of daily radiation image series over the study area for a relatively long time period.
The improvements brought by ERAD do not concern other limits of the LSA SAF product, such as the assumption of isotropic and constant albedo within each 5 km pixel. The approximations brought by such an assumption should be of minor importance in most situations, and the same is true for the improvements brought by the possible use of higher spatial resolution albedo estimates, such as those produced by the MODIS system (Carrer et al., 2010).
Finally, in the near future, ERAD will be able to use, in place of the current DSSF LSA-201 product, the new DSSF LSA-207 product, which became operational from April 3, 2020. The latter product, available every 15 min, shows some important improvements particularly concerning the definition of the aerosol scheme used in a cloudless atmosphere (Ceamanos et al., 2014a;2014b;Carrer et al., 2019a;2019b). Moreover, this new product directly provides the diffuse component of solar radiation derived from a standard decomposition model (Reindl et al., 1990). These characteristics are expected to improve the accuracy of LSA SAF data and, consequently, that of ERAD estimates.

| CONCLUSIONS
The quantification of solar radiation is crucial to characterize the global phenomena related to climate and climate change (Wild, 2016). At a more local level, the prediction of the solar radiation incoming on complex morphologies is decisive to improve evapotranspiration estimates (Aguilar et al., 2010), to predict crop and forest yield  and to plan the generation of electric power by photovoltaic systems (Demain et al., 2013).
The new method currently proposed to predict solar radiation over rugged terrains, ERAD, disaggregates the LSA SAF product by modelling the distribution of solar radiation over a digital elevation model having a spatial resolution of 200 m. The theoretical bases of ERAD have been thoroughly presented, followed by a testing of the model in comparison with other three standard methods (MTCLIM, LSA SAF and CAMS).
The results of this intercomparison indicate that, over horizontal surfaces, the solar radiation estimates produced by ERAD are notably more accurate that those of the other algorithm specifically conceived for operational application on rugged terrains, MTCLIM. The radiation estimates produced by ERAD have an accuracy similar to those of LSA SAF and CAMS, but higher spatial and temporal resolutions (200 m and 1 min). Consequently it can be presumed that also over tilted surfaces the new model is capable of providing solar radiation estimates having enhanced details in space and time. This quantitative assessment has been completed by a qualitative evaluation of the variability of the radiation estimates produced by the four methods, which supports the same conclusion.
ERAD is therefore reasonably efficient in reproducing the spatio-temporal variations of solar radiation over rugged areas. In particular, this method is suitable for achieving the main practical objective of the current investigation, which is the completion of an existing database containing medium-term series of temperature and rainfall observations for the entire Italian national territory (Maselli et al., 2012;Fibbi et al., 2016).