Two decades of distributed global radiation time series across a mountainous semiarid area (Sierra Nevada, Spain)

The main drawback of the reconstruction of high-resolution distributed global radiation (R_g) time series in mountainous semiarid environments is the common lack of station-based solar radiation registers. This work presents 19 years (2000–2018) of high-spatial-resolution (30 m) daily, monthly, and annual global radiation maps derived using the GIS-based model proposed by Aguilar et al. (2010) in a mountainous area in southern Europe: Sierra Nevada (SN) mountain range (Spain). The model was driven by in situ daily global radiation measurements, from 16 weather stations with historical records in the area; a 30 m digital elevation model; and 240 cloud-free Landsat images. The applicability of the modeling scheme was validated against daily global radiation records at the weather stations. Mean RMSE values of 2.63 MJ m⁻² d⁻¹ and best estimations on clear-sky days were obtained. Daily R_g at weather stations revealed greater variations in the maximum values but no clear trends with altitude in any of the statistics. However, at the monthly and annual scales, there is an increase in the high extreme statistics with the altitude of the weather station, especially above 1500 m a.s.l. Monthly R_g maps showed significant spatial differences of up to 200 MJ m⁻² per month that clearly followed the terrain configuration. July and December were clearly the months with the highest and lowest values of R_g received, and the highest scatter in the monthly R_g values was found in the spring and fall months. The monthly R_g distribution was highly variable along the study period (2000–2018). Such variability, especially in the wet season (October–May), determined the interannual differences of up to 800 MJ m⁻² yr⁻¹ in the incoming global radiation in SN. The time series of the surface global radiation datasets here provided can be used to analyze interannual and seasonal variation characteristics of the global radiation received in SN with high spatial detail (30 m). They can also be used as cross-validation reference data for other global radiation distributed datasets generated in SN with different spatiotemporal interpolation techniques. Daily, monthly, and annual datasets in this study are available at https://doi.org/10.1594/PANGAEA.921012 (Aguilar et al., 2021).

1 Introduction

High-mountain areas in semiarid environments present singular characteristics due to the continuous interaction of alpine conditions in the summits with the surrounding semiarid climate. They play a key role as water providers during the warm and dry season when they often constitute the only water source for many rivers. Here, water fluxes from the snowpacks show a shift from the predominant partition between snowmelt and sublimation usually found in colder and wetter climates on an annual and seasonal basis (Herrero and Polo, 2016). This shift is caused by the radiation balance that enhances sublimation during cold and dry periods and intense snowmelt rates during late winter and spring in these areas (MacDonell et al., 2013; Liu et al., 2019). However, weather stations are not always equipped to monitor the global radiation or their components, and, moreover, they are seldom found in high altitudes, especially over 1500 m a.s.l., which makes it difficult to accurately assess not only the solar radiation temporal regime but also the spatial patterns of solar radiation fields in high-mountain areas. This impacts the availability of data for studies in mountains dealing with climate and hydrology, global warming, ecosystem services provided by the snow areas, and environmental and social and economic impacts on-site and downstream (Yang et al., 2010; M. Liu et al., 2012; Tang et al., 2019). It is not surprising that many mountain regions are identified as biodiversity hotspots around the world, with Mediterranean and other semiarid to arid regions being highly represented (Myers et al., 2000; O'Farrell et al., 2010; Hewitt, 2011; Pauli et al., 2012).

There are several research papers on solar radiation estimations from routine ground-based observations in high-altitude regions (Dubayah and van Katwijk, 1992; Dubayah, 1994; Tovar et al., 1995; Oliphant et al., 2003; Tovar-Pescador et al., 2006; Yang et al., 2006, 2010; Batlleìs et al., 2008; Bosch et al., 2008; Sheng et al., 2009; Aguilar et al., 2010; Mamassis et al., 2012; Chen et al., 2013; Zhang et al., 2020). All of them insist on the need to consider topographic effects and advise of the errors that simple interpolation/extrapolation techniques can create. Radiation data obtained from a dense and properly maintained weather station network in mountainous areas are rarely available, and therefore, modeling techniques need to be applied. M. Liu et al. (2012) state that the most difficult issue in solar radiation modeling in data-sparse regions is cloud accounting, due to the rapid spatially and temporally changing weather conditions and the three-dimensional structure of clouds. This complexity adds to the heterogeneity resulting from shadowing and reflection due to steep topography (Dubayah, 1992; Batlleìs et al., 2008; Mamassis et al., 2012; Chen et al., 2013; Zhang et al., 2019, 2020).

According to Dubayah and Rich (1995), as solar radiation models become more complex, they can be more difficult to use, mainly because of the requirement for additional input data. In fact, the complexity of physically based solar radiation formulations for topography and the lack of the data needed to drive such formulations led in the past to the lack of suitable modeling tools (Dubayah, 1994). Thus, it is important that the models allow for some flexibility regarding the component of radiation calculated and the input data needed.

Excluding traditional interpolation methods, there are two major methods for solar radiation modeling, namely satellite-derived solar radiation estimates and geographic information system (GIS)-based solar radiation models. Satellite-derived solar radiation models provide a wide spatial and temporal coverage but low spatial resolution when dealing with pixels with a strong topographic gradient. By contrast, GIS-based models calculate the incoming solar radiation for each cell of a digital elevation model (DEM) and allow for higher spatial resolutions including topographic effects. In the past decades, several models based on GIS have been proposed (e.g., Dubayah and Rich, 1995; Fu and Rich, 2000a, 2002; Wilson and Gallant, 2000; Goldberg and Häntzschel, 2002; Šùri and Hofierka, 2004; M. Liu et al., 2012; Zhang et al., 2019, 2020). Required input data include digital elevation values and atmospheric attenuation parameters that are commonly estimated from ground-based measurements and/or satellite data (Dubayah, 1994).

The aim of this study was to generate the spatiotemporal distribution of global solar radiation in a high-mountain semiarid area in southern Spain with a modeling scheme that reconstructs time map series from the usually available weather datasets. For this purpose, a GIS-based topographic solar radiation model (Aguilar et al., 2010) was applied in Sierra Nevada (SN) (Spain), a high mountain range running west–east parallel to the Mediterranean coastline with influence from both the sea and the African continent to the south and the continental conditions to the north. The accuracy of solar radiation estimates by the model was evaluated in terms of the error in the approximation to observed data. This study site is a high-value environmental area declared a “biosphere reserve” by UNESCO in 1986 due to the exceptional presence of endemisms (Heywood, 1995; Blanca et al., 1998; Anderson et al., 2011; Cañadas et al., 2014). In addition, SN is also included in the global change observatories network given its singular location between two seas and two continents, as well as its extreme topographic gradients (Bonet-García et al., 2015).

This paper presents 19 years of daily, monthly, and annual solar radiation maps with high resolution (30 m) over SN. The huge number of members involved in the management of this area make this information valuable in different fields, such as hydrology – the crucial role of the energy budget in the hydrological cycle over this area; and ecology – ecological communities' behavior and development clearly link with the amount of energy available; production systems downstream – as hydropower facilities and traditional to tropical crop systems from the top to downhills. In addition, these datasets directly contribute, or are relevant for many studies that could do so, to 2 of the 23 unsolved problems in hydrology (UPH) recently posed by Blöschl et al. (2019) in a participatory analytical discussion among the scientific community: UPH 16 “How can we use innovative technologies to measure surface and subsurface properties, states and fluxes at a range of spatial and temporal scales?” and UPH 5 “What causes spatial heterogeneity and homogeneity in runoff, evaporation, subsurface water and material fluxes (carbon and other nutrients, sediments), and in their sensitivity to their controls (e.g., snowfall regime, aridity, reaction coefficients)?”.

Figure 1 Location of the study site in southern Spain (left). Digital elevation model (DEM) and weather stations in Sierra Nevada (SN) (right). The numbers correspond to the station codes.

2 Study site

The Sierra Nevada mountain range (SN) is located 35 km north of the Mediterranean Sea (Fig. 1) and constitutes a mountainous area of the Natura 2000 network. Elevations rise up from 262 to 3479 m a.s.l. in a 4583.72 km² area that runs parallel to the sea. High altitudinal gradients are representative of the area, with variation in elevation of about 3400 m in less than 40 km of horizontal distance and a mountain climate in the summits surrounded by Mediterranean climate in the lower areas. Thus, the interaction of such conditions creates a strong heterogeneity in terms of soil types, landforms, and vegetation species that determine a complex hydrological response in the area and many endemic species (Heywood, 1995; Blanca et al., 1998; Anderson et al., 2011). The rainfall regime is highly variable, even in consecutive years, with annual cumulative values in the period (1960–2000) that range between 200 mm in dry years to 1000 mm in wet years, with an average value of 510 mm (Pérez-Palazón et al., 2015). The temperature regime is also heterogeneous, with values of 26, 12.5, and 0.4 ^∘C for maximum, mean, and minimum daily temperature in the same period.

The snow presence becomes relevant from November above 2000 m a.s.l. and extends up to spring with conditions that make the activity of a major ski resort in the area possible. However, in some winters, mild episodes can be found in January and February that melt most of the snow much earlier than the mean end of the snow season in the area (Herrero et al., 2009; Herrero and Polo, 2012; Pimentel et al., 2013). Because of its singular characteristics and fragile environment, Sierra Nevada receives international recognition as a biosphere reserve (1986), a national park (1999), an important bird area (2003), a special area of conservation (2012) and one of the international global change observatories in mountain areas. These environmental protection figures together with the different and numerous members involved in the management of such a unique area have determined the strong effort in data collection in the last years to advance the knowledge of the different aspects that determine the dynamics of this natural system. Moreover, global warming impacts threaten the environmental values of this system but also the associated ecosystem services and social and economic activities due to the estimated shift of the snowfall regime (Pérez-Palazón et al., 2018).

3 Data

3.1 Input data

A digital elevation model (DEM) with 30 m spatial resolution and 1 m vertical precision was used in this study (Fig. 1). The DEM was provided by the Andalusian regional administration, and it was generated by digital stereo correlation of aerial photographs of the Spanish National Plan for Aerial Orthophotography. The DEM is used to calculate the slope, aspect, sky view factor, and terrain configuration maps that are used in the modeling process (Dozier and Frew, 1990).

Table 1 Information of the weather stations included in this study: elevation, z (m a.s.l.); code; data length, as initial and final dates of the time series; number of initially available daily records, N_o (days); number of available daily records after the quality check, N (days); rate of days for cloudy, N_CI<0.3 (%), partially cloudy, $N_{\mathrm{0.3}<\mathrm{CI}<\mathrm{0.6}}$ (%), and clear-sky conditions, N_CI>0.6 (%); and maximum, R_{go_max} (MJ m⁻² d⁻¹), mean, R_{go_mean} (MJ m⁻² d⁻¹), and minimum, R_{go_min} (MJ m⁻² d⁻¹), daily global radiation observed values. The selected descriptors for sky conditions and global radiation correspond to registered data after quality check.

Meteorological input data are the longest available in situ daily global radiation (R_go) of 16 weather stations over the area (Fig. 1 and Table 1). The extent of the records in all weather stations (N_o in Table 1) was considered long enough to carry out the evaluation process dating from February 2000 for the oldest station (608 in Table 1). Twelve out of the 16 weather stations are located above 1500 m a.s.l. and 7 of them above 2000 m a.s.l. (Fig. 1). The stations belong to four different organizations: the Department of Agriculture, Fisheries and Environment of the Andalusian Government (601–608 in Table 1), the Water and Environment Agency (1001 and 1002 in Table 1), the National Parks Organization (853–860 in Table 1), and the Guadalfeo Monitoring Network (802–804 in Table 1) described in Polo et al. (2019). Pyranometers used to collect the data were of different natures but all of them with a characteristic range of around 0.35–1.1 µm: Skye SP1110 (stations 601, 602, 604, and 608), Kipp & Zonen SP-Lite pyranometer (station 802), Hukseflux LP02 (station 803), Hukseflux NR01 (stations 1001, 1002, and 804), and Kipp & Zonen CNR1 net radiometer (stations 853, 854, 855, 857, 858, 859, and 860).

In order to generate the complete global radiation data series for the whole time span (1 February 2000–31 December 2018), we first apply a quality control check to the recorded data at the weather stations.

3.2 Data quality control

Numerous studies on quality control of measured solar radiation data can be found in the literature (Geiger et al., 2002; Younes et al., 2005; Moradi, 2009; Journée and Bertrand, 2011). Compared to other meteorological variables, solar radiation measurement is more prone to errors (Moradi, 2009). Younes et al. (2005) state two main sources of errors related to in situ measurement of solar radiation: those related to equipment and uncertainty and operational errors. Thus, prior to any computation two basic screenings were applied to recorded daily global radiation data to discard suspicious records associated with equipment and operational errors (Younes et al., 2005).

1. Observed daily global radiation (R_go) must be between the daily extraterrestrial radiation (R_ext) and a minimum 3 % of R_ext (Geiger et al., 2002; Moradi, 2009).

2. Observed daily global radiation (R_go) must be lower than the clear daily global radiation (R_gcs) observed under a highly transparent clear sky (Wu et al., 2007). R_gcs values were calculated with the model developed by Ineichen and Perez (2002) and the parameterization of Kasten and Young (1989) for the air mass. More detail regarding the equation as well as its parameters can be found in Aguilar et al. (2010).

The excluded values from these tests did not reach 1 % of the data at any weather station.

Figure 2 Data availability in the analyzed period (1 February 2000–31 December 2018) for each weather station. Stations are sorted by increasing altitude from the top to the bottom row.

A third quality control screening was applied following Younes et al. (2005) to detect erroneous data due to operational errors related with particularities of weather stations in high altitudes (e.g., shadows, impacts of snow, mechanical failures, etc.). They suggest a semi-automatic procedure that allows the creation of an expectancy envelope in the clearness index (CI) diffuse-to-global-irradiance ratio (k) domain to reject data too obviously erroneous. The CI data range is divided into bands of equal width, within which the mean and standard deviation of the k values, μ_k and σ_k, are calculated. The top and bottom boundary shapes are identified by fitting two polynomials through the points μ_k±bσ_k limited between 0 and 1 to respect the physical range of the CI. In this study b values between 2 and 3 were applied to limit both the rejection of good data and the acceptance of erroneous data to small percentages.

The CI was calculated with the observed data at each weather station. However, no measurements of daily diffuse radiation, R_d, were available. Thus, the model proposed by Aguilar et al. (2010) was applied to generate daily diffuse radiation (R_dp) at each weather station without considering the observed global data at such station. Obviously, this assumption depends on the validity of the model as well as on the quality of R_go datasets at the remaining weather stations. However, under the common lack of diffuse solar radiation measurements like the present one, modeling them can be an alternative (e.g., Yang et al., 2020) to reject erroneous R_g observations. This approach was proposed once the model had already been validated in a previous study (Aguilar et al., 2010) but keeping in mind the intrinsic limitations and assumptions previously stated.

After this quality test, the percentage of excluded values did not reach 10 % at any weather station, with a mean value close to 2 % when the whole set of stations was considered. Table 1 shows selected descriptors of the datasets at each station in this study after all the quality check process and Fig. 2 shows the chronogram of the final input data availability per station (N in Table 1) used in this study.

3.3 Generation of global radiation maps

The GIS-based solar radiation model proposed by Aguilar et al. (2010) that was previously implemented and validated in a small subwatershed located in the southwest of Sierra Nevada (Fig. 1) was extended to the whole area in this study. For validation purposes, data registered at weather stations are considered to represent the average values of the 30 m cell of the DEM on which they are located (Batllés et al., 2008; Martínez-Durbán et al., 2009).

The main equations and flowchart of the model are shown in Appendix A. The complete explanation of the algorithms as well as the justification of the assumptions of the model can be found in detail in Aguilar et al. (2010).

The model was developed to be run using limited data but considering the agents that constitute the main sources of the spatial and temporal variability of solar radiation. Results generated by the model include hourly maps of diffuse, beam, and reflected solar radiation values with minimum input data requirements as only topographic data, albedo estimations, and measured daily global radiation records (R_go) at least at one weather station are required. As for the daily global radiation registers, even when they are missing, their estimation from other more readily available meteorological data could always be a choice from the literature (Hargreaves and Samani, 1982; Bristow and Campbell, 1984; Allen, 1997; Bechini et al., 2000; Winslow et al., 2001; Donatelli et al., 2003, 2006; Yang and Koike, 2005; Diodato and Bellocchi, 2007; Wu et al., 2007; Ruiz-Arias et al., 2011; J. Liu et al., 2012; El Ouderni et al., 2013; Mullen et al., 2013).

The generation of global radiation maps with the model applied (Aguilar et al., 2010) requires a proper characterization of the spatiotemporal patterns of albedo in the study site. A total of 240 cloud-free Landsat images available for the study period from Landsat 5 TM (49 images), Landsat 7 ETM+ (141 images), and Landsat 8 OLI (50 images) with a 30 m spatial resolution were used. Figure 3 shows the specific dates and sensors of the 240 images analyzed in this study. All images were first properly corrected, and their reflectivity values were computed (Pimentel et al., 2014). Albedo was then derived for each image following the same procedure applied in Aguilar et al. (2010), which is based on the methodology described by Brest and Goward (1987), and linearly interpolated on a daily timescale for the whole study period.

Figure 3 Dates and sensors of each Landsat image analyzed in the study period (1 February 2000–31 December 2018).

Figure 4 Cross-validation analysis. Linear fits of daily predicted vs. observed R_g (MJ m⁻² d⁻¹) at each of the selected stations for the global data (black) and cloudy (CI < 0.3 – red), partly cloudy (0.3 < CI < 0.6 – blue), and clear-sky (CI > 0.6 – orange) days. Stations are sorted by increasing altitude from left to right and from the top to the bottom row.

3.4 Cross validation at weather stations

Once daily global radiation estimates were generated by the model a cross validation was applied at each weather station on the daily scale. This was carried out on a leave-one-out process; i.e., data from a weather station were removed from the input dataset to the model, and predicted values (R_gp) at that weather station were then compared to observed data (R_go).

Different indicators were computed to quantitatively evaluate the performance of the model (Muneer et al., 2007):

• The root mean square error (RMSE) (Eq. 1), where R_gp and R_go are the predicted and observed daily global radiation (MJ m⁻² d⁻¹), respectively, and N is the number of observed daily data. It measures the difference between values predicted by the model and those which were observed.

  $$\begin{array}{}\text{(1)}& \mathrm{RMSE}=\sqrt{{\displaystyle \frac{\sum {\left(R_{\mathrm{gp}}-R_{\mathrm{go}}\right)}^{\mathrm{2}}}{N}}}\end{array}$$

• The deviation from the 1:1 line of observed vs. predicted daily solar radiation values. Linear fits forced through the origin were obtained (Eq. 2), and the slopes (α in Eq. 2) are desired to be equal to 1. The coefficient of determination, R², as the ratio of the explained variation to the total variation, was also computed.

  $$\begin{array}{}\text{(2)}& R_{\mathrm{gp}}=\mathit{\alpha }\cdot R_{\mathrm{go}}\end{array}$$

The RMSE values and linear fits were obtained for the whole dataset at each weather station and also for different cloudiness levels to consider different atmospheric states that may condition the performance of the model according to previous studies (Batllés et al., 2008; Martínez-Durbán et al., 2009; Ruiz-Arias et al., 2009). Based on the cloudiness, three types of weather conditions were analyzed: cloudy days (CI < 0.3), partly cloudy days (0.3≤ CI < 0.6), and clear-sky days or cloudless days (CI ≥ 0.6).

The cross-validation analysis was also carried out with deseasonalized daily data to remove the expected intra-annual course of global radiation data. The deseasonalization of the daily series was carried out applying a stable seasonal filter (Brockwell and Davis, 2002) as already done in a previous study with other hydrometeorological datasets (Aguilar et al., 2017). In addition, as the reliability of solar radiation estimates is conditioned by the availability of recorded data, the cross-validation analysis for the whole study period was also computed with limited data. Thus, global radiation estimated was generated with only the four stations (601, 602, 604 and 608 in Fig. 1) with the longest records (Fig. 2) as inputs to the model. Results are shown in Appendix B.

Finally, to contrast the modeling scheme applied, another well-known GIS-based solar radiation model, Solar Analyst (SA) (Fu and Rich, 2000b), was also applied in the study site. Error values in the approximation to observed data and linear fits obtained in SN are shown in Appendix C.

The cross-validation assessment is summarized in Fig. 4. With the global datasets (in black in Fig. 4), a very close approximation of the model estimates to recorded data was obtained (mean α value of 0.98 and mean R² values of 0.91). RMSE values varied for the different stations and ranged from 1.81 (station 804) to 3.76 (station 860) with a mean value of 2.63 MJ m⁻² d⁻¹.

When the analysis was carried out in terms of the cloudiness level, a general overestimation by the model (e.g., a mean α value of 1.41) was always seen on cloudy days (CI ≤ 0.3). In contrast, on clear-sky days (CI > 0.6) slopes were very close to 1 with a mean α value of 0.96. An intermediate behavior was found on partly cloudy days (0.3 < CI ≤ 0.6) when the model slightly underpredicted (e.g., stations 854 and 608) or overpredicted depending on the weather station. As for RMSE values, the lowest values were always found for clear-sky days, when the cloud influence is minimal and the attenuation is mostly explained by changes in the atmospheric transmittance, followed by partly cloudy days with mean values of 2.07 and 3.07 MJ m⁻² d⁻¹, respectively. The highest RMSE values were always found on cloudy days with mean values of 3.70 MJ m⁻² d⁻¹. The high proportion of clear-sky days (65 %) and the low RMSE values on these days (2.07 MJ m⁻² d⁻¹) revealed the general good agreement of the model estimates with observed data. This is especially important in semiarid environments, where energy-limited hydrological processes (e.g., soil moisture depletion, evaporation, or snowmelt) are more relevant on clear-sky days, and they must be carefully computed in water and energy balance modeling, irrigation scheduling, etc. (Chen et al., 1999; Mamassis et al., 2012).

There is no clear pattern in the errors obtained with the elevation of the stations. The goodness of the model estimates was more affected by the interaction of the different characteristics of the weather station (e.g., slope, aspect, surrounding terrain configuration, orographic effects in the vertical development of clouds, etc.) than by the height of the station itself.

With the deseasonalized time series (Fig. B1), differences were reduced among the different cloudiness levels. The most remarkable change was a significant improvement in the estimates of cloudy days in every station when the range of RMSE values shifted from 2.54–7.52 (in red in Fig. 4) to 1.72–5.16 MJ m⁻² d⁻¹ (in red in Fig. B1). Also, the range of the slopes significantly narrowed from 1.18–1.74 (red α values in Fig. 4) to 0.92–1.09 (red α values in Fig. B1). Thus, the comparison with deseasonalized data showed a higher accuracy of the model than the one obtained with the original datasets (Fig. 4).

The comparison with limited input datasets shown in Fig. B2 confirmed the lower reliability of global radiation estimates in the first 5 years when datasets recorded at only four stations (601, 602, 604 and 608 in Fig. 1) were available in SN. Here, higher-elevation stations are subjected to a slightly greater overestimation of solar radiation (1.34–2.04 in red in Fig. B2), especially during cloudy conditions when the RMSE values increased to 3.62–8.45 MJ m⁻² d⁻¹ (in red in Fig. B2).

The application of Solar Analyst (Table C1) revealed in general worse approximations to observed data than those shown in Fig. 4, with mean R² values of 0.66 and RMSE values ranging from 3.59 (station 853) to 5.11 (station 859) with the global datasets.

The errors obtained in Fig. 4 were within the order of magnitude of those found in previous studies in other mountainous areas (Yang et al., 2006; 2010; Zhang et al., 2020) and slightly improved those previously obtained on a small subarea (10 × 5 km²) in the northeastern side of SN. Here, Tovar-Pescador et al. (2006) analyzed the application of SA in clear-sky days with 168 daily global radiation values in the dataset from 14 weather stations located between 1091 and 1659 m a.s.l. They obtained R² values of 0.75, similar to the value here obtained with SA estimates in the whole SN area (0.77 in Table B1) but lower than the R² equal to 0.99 obtained with the model (in orange in Fig. 3). Then, Batllés et al. (2008) in another application of SA in the same area with a 2-year daily dataset obtained the best performances for clear-sky days. RMSE values obtained in clear-sky days in the present study, of 11.1 % (2.07 MJ m⁻² d⁻¹), were the same as those obtained by Batllés et al. (2008) for clear-sky days (11 %). Later, Ruiz-Arias et al. (2009) evaluated the application of four different GIS-based solar radiation models with 523 daily global radiation values in the dataset at the same study site. RMSE values for the global dataset ranged between 1.99 and 7.28 MJ m⁻² d⁻¹ depending on the model.

Figure 5 Statistical distribution of daily R_g (MJ m⁻² d⁻¹) time series at each of the selected stations over the study area. The box shows 50 % of the data, delimited by Q1 (lower) and Q3 (upper), the solid line represents the median, and whiskers show 10th and 90th percentiles. Brown, orange, and yellow dots represent daily maximum, mean, and minimum time series values. Stations are sorted by increasing altitude from left to right.

The order of magnitude of the errors (Fig. 4) and its comparison with those obtained with more computationally and data-demanding GIS-based models in previous studies led us to conclude that the model is the best choice to generate global radiation data series in SN.

Therefore, once the model was validated in the study site, daily R_g maps were generated and aggregated at the monthly and annual scales.

4 Results

Daily, monthly, and annual R_g datasets in SN are analyzed in this section at two spatial scales. First, the results at the weather station scale are presented. Thus, possible relationships between altitude and/or location of the weather station with the different R_g statistics and how this relation changes with the temporal scale of analysis can be assessed. Then, the R_g maps that can be downloaded as specified in Sect. 5 are analyzed.

4.1 Daily time series of global radiation in Sierra Nevada

Figure 5 shows the statistical distribution of the daily R_g at each weather station ordered by increasing altitude and illustrates several questions. First, there is a very similar interquartile range among stations. Second, there are greater variations in the maximum daily R_g among the different stations with a mean value of 34.0 MJ m⁻² d⁻¹. Third, even though a slight increase with altitude can be shown in the high extreme statistics of the daily R_g values (e.g., in the maximum or in the 90th percentile), there is not a clear trend. Therefore, other factors such as orientation, proximity to the sea, or the terrain configuration in the surrounding terrain as suggested by Batllés et al. (2008) constitute relevant features in the study site.

Figure 6 shows an example of the spatial distribution in three representative days of cloudy, partially cloudy, and clear-sky conditions. Here the spatial distribution is clearly influenced by the topography of SN, especially in the clear-sky day.

Figure 6 Daily R_g (MJ m⁻² d⁻¹) in SN for three selected days that represent the three levels of cloudiness considered in this study: cloudy, partially cloudy, and clear sky.

Figure 7 Statistical distribution of monthly R_g (MJ m⁻² per month) time series at each of the selected stations over the study area. The box shows 50 % of the data, delimited by Q1 (lower) and Q3 (upper); the solid line represents the median; and whiskers show the 10th and 90th percentiles. Brown, orange, and yellow dots represent monthly maximum, mean, and minimum time series values. Stations are sorted by increasing altitude from left to right and from the top to the bottom row.

4.2 Monthly time series of global radiation in Sierra Nevada

The statistical distribution of monthly R_g per weather station (Fig. 7) shows that in every station (i) July and December constitute the months with the highest and the lowest values of R_g, respectively; (ii) there is a quite linear increase in the monthly R_g values from January to July and a sudden drop in August with a slightly convex evolution till December; and (iii) the interquartile range is significantly higher in the spring and fall than in the summer and winter months.

The increase in the high extreme statistics of radiation with the altitude of the weather station becomes more apparent at the monthly scale (Fig. 7) than at the daily scale (Fig. 5) previously analyzed. Thus, maximum values of around 1000 MJ m⁻² per month are reached in July in the highest stations (e.g., 1002, 802, 1001, and 860 in Fig. 7), whereas this value decreases to around 910 MJ m⁻² per month in the four lowest stations except for station 608.

Monthly R_g maps show significant spatial differences of up to 200 MJ m⁻² per month in both the mean monthly values (Fig. 8) that clearly follow the terrain configuration, with summits and valleys receiving high and low solar radiation values, respectively. For example, the area in the north of SN that is highly shadowed by the highest peaks in the Iberian Peninsula (Mulhacen and Veleta with 3482 and 3396 m a.s.l., respectively) is easily visible, with the lowest relative levels of insolation received within SN especially in the summer months (June, July, and August in Fig. 8).

Figure 8 Monthly average and standard deviation of R_g (MJ m⁻² per month) in the study period (2000–2018) in SN.

Both maps of the monthly mean and standard deviation of R_g (Fig.  8) and the statistical distribution of the monthly R_g in the study site (Fig. 9) show the same behavior as the one obtained at the weather stations regarding (i) July and December as the months with the highest and lowest values of R_g received in SN and (ii) the highest scatter in the monthly R_g values in the spring and fall months.

Figure 9 Statistical distribution of the monthly R_g (MJ m⁻² per month) values throughout the study area. Whisker boxes represent the 10th, 25th, 50th, 75th, and 90th percentiles of each monthly map per year.

For the study period (2000–2018), there is a great heterogeneity in the statistical distribution of the monthly R_g in the study site (Fig. 9), especially in the incoming radiation along the months of the wet season (October–May). In this way, in the most insolated years in the study period (2005 and 2012), significantly higher monthly radiation values were found in certain months of the springtime (March and May 2012 and April 2005). In those months, the higher-than-usual rate of clear-sky over cloudy days finally determines the annual differences in the incoming global radiation in SN.

When considering the temporal evolution of the distribution of R_g within the monthly maps in SN (Fig. 10), certain interannual differences can be observed along the study period, such as the existence of certain months in spring with unexpected low monthly radiation values (e.g., 2001, 2004, 2007 and 2008) or two relative maximum monthly R_g values (e.g., 2009, 2010, and 2014). Moreover, Fig. 10 shows a higher scatter in the monthly maximum (June–August) and minimum (November–January) R_g values in SN than when the analysis is carried out at each weather station (Fig. 7).

Figure 10 Evolution of the statistical distribution of monthly R_g (MJ m⁻² per month) in the study period (2001–2018) throughout the study area. Grayscale colors represent the following percentiles: 5th, 10th, 25th, 50th, 75th, 90th, and 95th.

4.3 Annual times series of global radiation in Sierra Nevada

Unlike at the daily scale (Fig. 5), a great variability among the different weather stations in terms of the global radiation received at the annual temporal scale is found (Fig. 11). Thus, we find minimum annual R_g values from 5920 MJ m⁻² yr⁻¹ in station 854 to around 6750 MJ m⁻² yr⁻¹ in station 1002. This difference is even bigger in the maximum annual R_g values from 6700 to 7720 MJ m⁻² yr⁻¹ in stations 854 and 802, respectively, and is also shown in the interquartile range.

Figure 11 Statistical distribution of annual R_g (MJ m⁻² yr⁻¹) time series at each of the selected stations over the study area. The box shows 50 % of the data, delimited by Q1 (lower) and Q3 (upper); the solid line represents the median; and whiskers show 10th and 90th percentiles. Brown, orange, and yellow dots represent annual maximum, mean, and minimum time series value. Stations are sorted by increasing altitude from left to right.

When analyzing the influence of altitude, the weather stations above 1500 m a.s.l (854, 857, 859, 804, 855, 858, 1002, 802, 1001, and 860 in Fig. 11) show their altitudinal gradient in all the statistics of the annual R_g values considered.

Annual R_g maps (Fig. 12) show the same spatial differences that follow the terrain configuration as those observed in the monthly time series (Fig. 8). For example, the area in the north of SN that is highly shadowed as previously mentioned corresponds to the area with the mean minimum annual values received in the study period, 4063 MJ m⁻² yr⁻¹, that only represents 63 % of the mean annual accumulated values in SN (6316 MJ m⁻² yr⁻¹).

Figure 12 Annual global radiation (MJ m⁻² yr⁻¹) in the study period (2001–2018) in SN.

Significant interannual differences can be easily shown with differences in the mean annual R_g value in the study area of up to 800 MJ m⁻² yr⁻¹ between 2005 and 2018. Such years with particularly high and low annual incoming radiation also presented higher (6800 MJ m⁻² yr⁻¹) and lower median annual R_g values (6200 MJ m⁻² yr⁻¹), respectively, than the annual median for the whole study period in SN (6456 MJ m⁻² yr⁻¹) (Fig. 13). These results agree with the annual irradiation map obtained by Batllés et al. (2008) in the northeastern part of SN. They reported maximum and minimum annual values of 7516 and 2342 MJ m⁻² yr⁻¹ on the summits and in deep valleys, respectively, and thus concluded that irradiation levels were more related to topographic characteristics than to altitude.

Figure 13 Evolution of the statistical distribution of annual R_g (MJ m⁻² yr⁻¹) in the study period (2001–2018) throughout the study area. Dashed lines represent the mean values of the percentiles analyzed.

5 Data availability

The daily, monthly, and annual global radiation maps derived in this study can be accessed and downloaded in .ncdf format from https://doi.org/10.1594/PANGAEA.921012 (Aguilar et al., 2021). In addition, a .txt file containing the availability (code 1) or gaps (code 0) in the daily R_g observations at each weather station has been added as a Supplement to this paper. Hourly datasets were also computed in this study, but due to their large storing capacity requirements they have not been included in the data repository specified above. Thus, hourly maps can be provided for certain dates upon request to the authors. However, a validation of these hourly datasets like the one applied in the daily estimates at the weather stations has not been specifically carried out in this study. Therefore, in case hourly maps are requested to the authors, these data should be taken with caution as the only available validation in SN was carried out at one weather station (802 in Fig. 1) and for a shorter period (2004–2010) in Aguilar et al. (2010).

6 Final remarks

This study presents 19 years (2000–2018) of daily, monthly, and annual global radiation maps of high spatial resolution (30 m) in a high-mountain Mediterranean site. In these areas the common lack of weather stations in high altitudes makes it difficult to accurately assess solar radiation spatial patterns.

A GIS-based modeling scheme based on measurements or estimations of incoming daily global radiation was applied and validated in the 16 weather stations available at this unique study site. Mean RMSE values ranged from 1.81 to 3.76 MJ m⁻² d⁻¹, depending on the weather station. The best estimations were always obtained on clear-sky days, when mean RMSE values decreased to 2.07 MJ m⁻² d⁻¹. The largest errors were obtained on cloudy days, which constitute on average 10 % of the daily datasets, and, therefore, future research should be conducted in order to improve the estimations in these situations, keeping the minimum input data requirement (daily global radiation data) advantage of the model. However, the high proportion (65 %) of clear-sky days, and the low RMSE values on those days, allow one to conclude that there is a good agreement between the model estimates and observed data in the study site.

Spatial differences of around 2000 MJ m⁻² yr⁻¹ were found within each year analyzed. In addition, significant differences were easily shown between the years in mean incoming values of up to 800 MJ m⁻² yr⁻¹. Those differences were mostly due to the variability in the incoming radiation at the wet season (October–May), with higher rates of clear-sky days in the most insolated years (e.g., 2005).

Thus, we can affirm that the modeling scheme here applied is an efficient option in semiarid mountainous areas, where daily global radiation datasets constitute the only source of solar radiation data.

Time series of these surface global radiation datasets can be used to analyze interannual and seasonal variation characteristics of the global radiation received in SN with high spatial detail (30 m). The availability of long global radiation datasets allows us to capture the annual variability within each cycle of the sun activity, as reported in the literature (Scaffetta and Wilson, 2013), and thus estimate its contribution to the annual variability of other climate variables in these semiarid mountainous areas.

Dense and properly maintained weather station networks in mountainous areas are rarely available. Thus, these datasets can also be used as cross-validation reference data for other global radiation distributed datasets generated in SN with different spatiotemporal interpolation techniques.

These results can also assess the order of magnitude of different sources of spatial variability (altitude/slope/aspect gradients) as well as the seasonal range of variation at different timescales and their annual variability. This estimation may provide a first estimate of the order of magnitude of uncertainty of average calculations or spatial interpolation from a scarce number of weather stations in Mediterranean and semiarid mountain areas.

The correct assessment of the solar radiation regime is crucial to correctly determine the temporal evolution of energy-limited hydrological processes such as the snow layer dynamics, soil moisture depletion, and evapotranspiration (Tomas-Burguera et al., 2019). Thus, as a key input parameter for the water and energy balance, these high-spatial-resolution solar radiation time series are useful not only for research on the snow domain and water planning in SN in the application of hydrological modeling, but also in many other applications – for example, for estimations of evapotranspiration for irrigation scheduling in the agricultural sector, for ecological and biodiversity-related studies, in the design and location of stand-alone solar energy facilities, and for outdoor recreational activities in the area that strongly rely on the hydrometeorological conditions of SN. Finally, this work contributes to feed research related to some key questions in hydrology, such as UPH 16 and UPH 5 identified by Blöschl et al. (2019).

Appendix A: Solar radiation equations

The sequence followed by the model is summarized in Fig. A1. Computations are classified at the point scale of weather stations (Point) and the distributed scale of grids of the digital elevation model (DEM) (Distributed). The complete explanation of the algorithms and assumptions of the model can be found in detail in Aguilar et al. (2010).

Figure A1 Flowchart of the solar radiation model.

Firstly, daily extraterrestrial radiation (R_ext in MJ m⁻² d⁻¹) is computed by integrating the extraterrestrial radiation incident upon a horizontal surface relative to the sun's beams from sunrise to sunset (Eq. A1).

$$\begin{array}{}\text{(A1)}& R_{\mathrm{ext}}=E_o\cdot I_{\mathrm{SC}}\cdot \cos \left({\mathit{\theta }}_{\mathrm{z}}\right),\end{array}$$

where I_SC is the solar constant (1367 W m⁻²), θ_z is the zenith angle, and E_o is the eccentricity factor. These variables were computed following the equations in Dozier et al. (1981).

Then, the daily clearness index (CI), as the ratio of observed daily global radiation (R_go in MJ m⁻² d⁻¹) to the daily extraterrestrial radiation, is computed at each weather station (Eq. A2).

$$\begin{array}{}\text{(A2)}& \mathrm{CI}={\displaystyle \frac{R_{\mathrm{go}}}{R_{\mathrm{ext}}}}\end{array}$$

CI is expressed in terms of two factors, CI_CS and fCI_cl. The first term represents the influence of atmosphere under clear-sky conditions over solar radiation, while the second term includes the cloudiness effects that decrease the final incoming solar radiation (Eq. A3). The approximation of Ineichen and Pérez (2002) is used to compute the global radiation under clear-sky conditions, R_gcs, and thus distributed hourly R_gcs values are obtained from the sun elevation angle, the height of the cell, the Linke turbidity factor (T_L), and the atmospheric mass obtained following the parameterization of Kasten and Young (1989). Thus, hourly CI_CS values can be computed cell by cell, and then the mean daily distributed values are generated. Once daily CI and CI_CS values are known, fCI_cl is obtained at each weather station from Eq. (A3) and spatially interpolated following the inverse distance weighted (IDW) method. From daily CI_CS and fCI_cl maps, daily interpolated CI and R_g values can be obtained at the cell scale from Eqs. (A3) and (A4).

$$\begin{array}{}\text{(A3)}& {\displaystyle }\mathrm{CI}={\mathrm{CI}}_{\mathrm{CS}}\cdot {\mathrm{fCI}}_{\mathrm{cl}}\text{(A4)}& {\displaystyle }{R}_{\mathrm{g}}=R_{\mathrm{ext}}\cdot {\mathrm{CI}}_{\mathrm{CS}}\cdot {\mathrm{fCI}}_{\mathrm{cl}}\end{array}$$

Topographic effects need to be evaluated for the different sun positions during the day, and thus hourly values of the different components need to be derived. Two different procedures are currently available in the model. The first one proposed in Aguilar et al. (2010) applies Jacovides et al. (1996) (Eq. A5a) to produce the daily diffuse (R_d in MJ m⁻² d⁻¹) and daily beam values (R_b in MJ m⁻² d⁻¹). The model finally computes hourly beam and diffuse values on horizontal surfaces (r_b and r_d, both in MJ m⁻² h⁻¹), from the daily amounts and following the temporal pattern of extraterrestrial hourly radiation during the day.

$$\begin{array}{}\text{(A5a)}& {\displaystyle \frac{R_{\mathrm{d}}}{R_{\mathrm{g}}}}=\left\{\begin{array}{ll}\mathrm{0.992}-\mathrm{0.0486}\mathrm{CI}& \mathrm{CI}\le \mathrm{0.1}\\ \mathrm{0.954}+\mathrm{0.734}\mathrm{CI}-\mathrm{3.806}{\mathrm{CI}}^{\mathrm{2}}& \\ +\mathrm{1.703}{\mathrm{CI}}^{\mathrm{3}}& \mathrm{0.1}<\mathrm{CI}\le \mathrm{0.71}\\ \mathrm{0.165}& \mathrm{CI}>\mathrm{0.71}\end{array}\right.\end{array}$$

The second approach uses the temporal pattern of extraterrestrial hourly radiation, r_ext, to generate hourly global values, R_g, according to previous studies (Chen et al., 1999; Ruiz-Arias et al., 2011). Then, the hourly regressive model developed by Ruiz-Arias et al. (2010) is applied to estimate the hourly diffuse values (Eq. A5b) from the hourly CI, CI_h, as the ratio of R_g to r_ext. This model was implemented as it has been validated over Europe and the USA using ground data from different sites, including some Spanish stations (Ruiz-Arias et al., 2010). Hourly beam values (r_b) are thus obtained on a cell basis as the difference between global and diffuse hourly radiation distributed values.

$$\begin{array}{}\text{(A5b)}& {\displaystyle \frac{r_{\mathrm{d}}}{r_{\mathrm{g}}}}=\mathrm{0.952}-\mathrm{1.041}{e}^{-\exp \left(\mathrm{2.3}-\mathrm{4.702}\cdot {\mathrm{CI}}_{\mathrm{h}}\right)}\end{array}$$

First applications at the study site have shown negligible differences between both partitioning schemes. The differences with daily recorded data were insignificant in the second decimal place of error values. Thus, the results presented in this study were obtained with the original scheme of Aguilar et al. (2010) (Eq. A5a) while the authors continue working on the improvement on the partitioning scheme of the model.

Then, the topographic correction is carried out, and, depending on the component, different procedures are applied.

Hourly beam radiation on a surface of slope β and orientation γ (r_b,βγ in MJ m⁻² h⁻¹) is calculated according to Eq. (A6) in terms of r_b, θ_z, and a new corrected zenith angle for the sloping surface, θ (Iqbal, 1983). Then, the model checks the shading effects. Self-shading will occur if the angle between the normal to the surface and the solar vector is greater than 90^∘. Finally, shading by nearby terrain takes place when the illumination angle is greater than the horizon angle in the same direction. The model previously obtained the horizons following the algorithms of Dozier et al. (1981) and Dozier and Frew (1990) by comparing the slopes between cells in the eight directions.

$$\begin{array}{}\text{(A6)}& r_{\mathrm{b},\mathit{\beta }\mathit{\gamma }}=r_{\mathrm{b}}\left({\displaystyle \frac{\cos \mathit{\theta }}{\cos {\mathit{\theta }}_{\mathrm{z}}}}\right)\end{array}$$

Hourly diffuse radiation on a surface of slope β and orientation γ (r_d,βγ in MJ m⁻² h⁻¹) is calculated according to Eq. (A7) in terms of r_d and SVF, the sky view factor that modifies the incoming radiation incident on a flat surface to consider possibly obstruction effects on a sloping surface (Dubayah, 1992). Dozier and Frew (1990) obtained an analytical expression for the estimation of the SVF in terms of the different horizons in each direction considered assuming an isotropic sky.

$$\begin{array}{}\text{(A7)}& r_{\mathrm{d},\mathit{\beta }\mathit{\gamma }}=r_{\mathrm{d}}\cdot \mathrm{SVF}\end{array}$$

Finally, hourly reflected radiation on a surface of slope β and orientation γ (r_r,βγ in MJ m⁻² h⁻¹) and albedo ρ is calculated according to Dozier and Frew (1990) as expressed in Eq. (A8).

$$\begin{array}{}\text{(A8)}& r_{\mathrm{r},\mathit{\beta }\mathit{\gamma }}=\mathit{\rho }\cdot \left[{\displaystyle \frac{\mathrm{1}+\cos \mathit{\beta }}{\mathrm{2}}}-\mathrm{SVF}\right]\cdot \left(r_{\mathrm{d}}+r_{\mathrm{b}}\right)\end{array}$$

Hourly global distributed radiation (r_gp in MJ m⁻² h⁻¹) is obtained by addition of the three hourly components at each cell according to Eq. (A9).

$$\begin{array}{}\text{(A9)}& r_{\mathrm{gp}}=r_{\mathrm{b},\mathit{\beta }\mathit{\gamma }}+r_{\mathrm{d},\mathit{\beta }\mathit{\gamma }}+r_{\mathrm{r},\mathit{\beta }\mathit{\gamma }}\end{array}$$

Finally, daily global distributed radiation (R_gp in MJ m⁻² d⁻¹) is obtained as the summation of hourly global distributed radiation values (Eq. A10).

$$\begin{array}{}\text{(A10)}& R_{\mathrm{gp}}=\sum _{\mathrm{24}h}{r}_{\mathrm{gp}}\end{array}$$

Appendix B: Cross validation with deseasonalized data series and limited input data

Figure B1 Linear fits of daily deseasonalized predicted vs. observed R_g (MJ m⁻² d⁻¹) at each one of the selected stations for the global data (black) and cloudy (CI < 0.3 – red), partly cloudy (0.3 < CI < 0.6 – blue), and clear-sky (CI > 0.6 – orange) days. Stations are sorted by increasing altitude from left to right and from the top to the bottom row.

Figure B2 Linear fits of daily predicted with limited input data vs. observed R_g (MJ m⁻² d⁻¹) at each one of the selected stations for the global data (black) and cloudy (CI < 0.3 – red), partly cloudy (0.3 < CI < 0.6 – blue), and clear-sky (CI > 0.6 – orange) days. Stations are sorted by increasing altitude from left to right and from the top to the bottom row.

Appendix C: Comparison with Solar Analyst estimates

Solar Analyst (SA) is one of the most used GIS-based solar radiation models. It calculates the insolation across a landscape or for specific locations, based on the methods developed by Fu and Rich (2000a, 2000b, 2002). The total amount of radiation is given as global radiation and depends on the latitude of the site, topography, shadow cast, and atmospheric attenuation. Global radiation is computed in SA as the sum of direct and diffuse radiation. The equations and modeling scheme can be found in detail in Fu and Rich (2000b).

Daily global radiation time series were generated in the study site with the Points Solar Radiation tool of SA at each weather station. Then a cross validation on a leave-one-out process was applied. The main inputs to the model, the diffuse fraction, k, and the atmospheric transmittivity, τ, were estimated in the study site from observed global radiation data following Batllés et al. (2008).

Table C1 Model performance with Solar Analyst: slope (α) and r² of the linear fit between daily predicted vs. observed R_g and RMSE (MJ m⁻² d⁻¹).

SA underestimated observed daily values with a mean α value of 0.78 in the global datasets (Table C1). In terms of the cloudiness level, a general underestimation by SA was always seen on cloudy (CI ≤ 0.3) and clear-sky days (CI > 0.6) with slopes of the fits significantly lower than 1 (mean α values of 0.42 and 0.74 respectively). In contrast, a slight overestimation with a mean α value of 1.03 was found on partly cloudy days (0.3 < CI ≤ 0.6). As for RMSE values, the lowest mean values were always found for cloudy days (1.59 MJ m⁻² d⁻¹), which are also lower than those obtained in Fig. 4 (3.70 MJ m⁻² d⁻¹). However, despite the lower RMSE values the deviation from the 1:1 linear fit in cloudy days with SA estimates was significant (mean α value of 0.42 and R² value of 0.39 in Table C1). The highest RMSE values with SA estimates were always found on partly cloudy days with a mean value of 5.35 MJ m⁻² d⁻¹ followed by clear-sky days with a mean RMSE value of 3.37 MJ m⁻² d⁻¹, both considerably higher than those obtained in Sect. 3.4 (3.07 and 2.07 MJ m⁻² d⁻¹ respectively).

Appendix D: Nomenclature

Symbols
CI:	daily clearness index
CICS:	daily clearness index in a cloudless atmosphere
CIh:	hourly clearness index
Eo:	eccentricity factor
fCIcl:	cloudiness effects factor
ISC:	solar constant
k:	diffuse-to-global-irradiance ratio
N CI < 0.3:	rate of days for cloudy conditions
N 0.3 < CI < 0.6:	rate of days for partially cloudy conditions
N CI > 0.6:	rate of days for clear-sky conditions
No:	number of initially available daily records in the study period
N:	number of available daily records after the quality check
Q1:	quartile 1
Q3:	quartile 3
Rb:	daily beam radiation
Rbp:	daily beam radiation predicted by the model
Rd:	daily diffuse radiation
Rdp:	daily diffuse radiation predicted by the model
Rext:	daily extraterrestrial radiation
Rg:	global radiation
Rgcs:	global radiation under clear-sky conditions
Rgo_max:	maximum daily global radiation observed value
Rgo_mean:	mean daily global radiation observed value
Rgo_min:	minimum daily global radiation observed value
Rgp:	daily global radiation predicted by the model
rb:	hourly beam radiation on horizontal surfaces
rb,βγ:	hourly beam radiation on a surface of slope β and orientation γ
rd:	hourly diffuse radiation on horizontal surfaces
rd,βγ:	hourly diffuse radiation on a surface of slope β and orientation γ
rext:	hourly extraterrestrial radiation
rr,βγ:	hourly reflected radiation on a surface of slope β and orientation γ
Rg:	hourly global radiation on horizontal surfaces
rgp:	hourly global radiation predicted by the model
R2:	coefficient of determination
TL:	Linke turbidity factor
z:	elevation
Abbreviations
DEM:	digital elevation model
IDW:	inverse distance weighted
RMSE:	root mean square error
SA:	Solar Analyst
SN:	Sierra Nevada mountain range
SVF:	sky view factor
UPH:	unsolved problems in hydrology
Greek symbols
α:	slope of the fit between Rgp and Rgo
β:	slope
γ:	orientation
μk:	mean of the diffuse-to-global-irradiance ratio
ρ:	albedo
σk:	standard deviation of the diffuse-to-global-irradiance ratio
θ:	corrected zenith angle for the sloping surface
θz:	zenith angle