Topographic effects on solar radiation distribution in mountainous watersheds and their influence on reference evapotranspiration estimates at watershed scale

Distributed energy and water balance models require time-series surfaces of the climatological variables involved in hydrological processes. Among them, solar radiation constitutes a key variable to the circulation of water in the atmosphere. Most of the hydrological GIS-based models apply simple interpolation techniques to data measured at few weather stations disregarding topographic effects. Here, a topographic solar radiation algorithm has been included for the generation of detailed time-series solar radiation surfaces using limited data and simple methods in a mountainous watershed in southern Spain. The results show the major role of topography in local values and differences between the topographic approximation and the direct interpolation to measured data (IDW) of up to +42% and −1800% in the estimated daily values. Also, the comparison of the predicted values with experimental data proves the usefulness of the algorithm for the estimation of spatiallydistributed radiation values in a complex terrain, with a good fit for daily values ( R2 = 0.93) and the best fits under cloudless skies at hourly time steps. Finally, evapotranspiration fields estimated through the ASCE-Penman-Monteith equation using both corrected and non-corrected radiation values address the hydrologic importance of using topographicallycorrected solar radiation fields as inputs to the equation over uniform values with mean differences in the watershed of 61 mm/year and 142 mm/year of standard deviation. High speed computations in a 1300 km 2 watershed in the south of Spain with up to a one-hour time scale in 30 × 30 m2 cells can be easily carried out on a desktop PC. Correspondence to: C. Aguilar (caguilar@uco.es)


Introduction
There are several methods available for the development of digital elevation models for hydrological studies but regular grid structures provide the best compromise between accuracy and computational efficiency (Moore et al., 1991).For this, all the inputs to distributed hydrological modeling must be available at this spatial scale.Among such inputs to hydrological models, solar radiation plays an important role in most of the processes involved, as it is a key variable in the circulation of water from the earth's surface to the atmosphere, especially in Mediterranean regions.At a global scale, latitudinal gradients caused by the earth's rotation and translation movements are well-known.However, at a smaller scale, apart from cloudiness and other atmospheric heterogeneities, topography determines the distribution of the incoming solar radiation; variability in slope angle and slope orientation, as well as the shadows cast by topographic agents, can lead to strong local gradients in solar radiation (Dozier, 1980;Dubayah, 1992;Dubayah and van Katwijk, 1992), with the corresponding influence on the energy-mass balance of the snow cover and its evolution (Dubayah and van Katwijk, 1992;Li et al., 2002;Herrero et al., 2009), the vegetation canopy (Dubayah, 1994), the surface soil layer, surface water bodies (Annear and Wells, 2007), etc.
The regional climate in Mediterranean areas is characterized by great inland-coast, valley-hill contrasts, and is subject to cyclical fluctuations in cloud cover, precipitation and drought, thus exhibiting considerable spatial and temporal variations (Diodato and Bellocchi, 2007).At these latitudes, during periods of a lack of rainfall -a common event at different spatial and temporal scales -radiation is the main force in the system which causes both snowmelt C. Aguilar et al.: Topographic effects on solar radiation distribution in mountainous watersheds and evapotranspiration.Here, an accurate estimation of timeseries solar radiation surfaces is required for distributed energy and water balance modeling (Ranzi and Rosso, 1995;Herrero et al., 2007).
One of the main drawbacks in the assessment of solar radiation is the lack of reliable data.In mountainous areas, where the monitoring network ineffectively covers the complex heterogeneity of the terrain, simple geostatistical methods for spatial interpolation are not always representative enough, and algorithms that explicitly or implicitly account for the features creating strong local gradients in the incoming radiation must be applied (Susong et al., 1999;Garen and Marks, 2005;Chen et al., 2007;Batllés et al., 2008;Yang et al., 2010).Thus, the implementation of the spatial variability in the incoming radiation at the cell scale for distributed hydrological modeling is of major concern, especially in mountainous areas (Allen et al., 2006).Here, the combination of extreme gradients in the spatial distribution of solar radiation, together with the lack of measurements at detailed spatial and temporal scales, calls for the integration of algorithms simple enough to be run with common measurements but at the same time able to capture the agents that constitute the main sources of the spatial and temporal variability of solar radiation.
At the local scale, the amount of solar radiation reaching a given location is called global solar radiation and it depends mainly on the cloud cover, the turbidity of the clean air, the time of the year, latitude, and surface geometry (Iqbal, 1983;Essery and Marks, 2007).As radiation penetrates the atmosphere, it is depleted by absorption and scattering.Not all of the scattered radiation is lost, since part of it eventually arrives at the surface of the earth in the form of diffuse radiation (Liu and Jordan, 1960).Global radiation is the sum of direct or beam radiation from the sun, diffuse radiation from the sky, where a portion of the overlying hemisphere may be obstructed by the terrain, and direct and diffuse radiation reflected by nearby terrains (Dubayah, 1994).Therefore, global radiation received on a surface with a random slope and aspect is largely controlled by atmospheric and topographic conditions (Flint and Childs, 1987;Tian et al., 2001;Diodato and Bellocchi, 2007).In a very rough terrain, some areas may not receive any direct radiation during the whole year -even if facing south -because of high peaks surrounding them.Under such conditions, a GIS-based solar radiation model that considers the impact of terrain shading should be applied (Allen et al., 2006).
In any case, for the estimation of radiation incident on tilted surfaces, the partition of global horizontal radiation into its beam and diffuse components is of major concern, as the topographic effects are different for each one and therefore have to be modeled separately (Iqbal, 1980;Antonic, 1998;González-Dugo et al., 2003).Thus, diffuse radiation is affected by the unobstructed portion of the overlying hemisphere, while reflected radiation is affected by terrain slope and the portion of the overlying hemisphere obstructed by the terrain (Dubayah, 1994).As for the beam component, self-shadowing and shadows cast by the surrounding terrain have to be considered for each sun position in the sky during the day.
For the quantification of the diffuse component, many parameters related to the atmospheric properties are required in order to express the scattering properties of the atmosphere.However, these parameters are not easily available or their computation from common measurements may be timeconsuming (Dubayah and van Katwijk, 1992), and so simpler procedures need to be applied, especially at watershed scale and in rough terrain.Thus, in these situations, the basic procedure for the partition of global radiation into its components is the calculation of correlations between the daily global radiation and its diffuse component from measured values of both quantities, and then applying these correlations at locations where diffuse radiation data are not available (Iqbal, 1980).In the literature, there are several reviews about the different correlations available, depending on the averaging procedure and on the time scale of the radiation data (e.g.hourly or daily) (Iqbal, 1980;Spencer, 1982;Kambezidis et al., 1994;Jacovides et al., 1996).Liu and Jordan (1960) were the first authors to develop a model for the estimation of diffuse radiation from global data, establishing the basis for later empirical analyses of global radiation from daily data.Ruth and Chant (1976) obtained a very similar figure and demonstrated a latitudinal dependence in the models.Other authors developed hourly correlations (Orgill and Hollands, 1977;Bugler, 1977;Erbs et al., 1982).In 1979, Collares-Pereira and Rabl, maintaining the assumption of isotropic approximation for the diffuse radiation previously proposed by Liu and Jordan, improved some aspects of the model (correction for the shade ring effect and use of daily values for extraterrestrial radiation instead of single monthly values), and defined the daily clearness index as the ratio of global radiation to extraterrestrial radiation.
Despite the availability of topographically corrected models for the estimation of solar radiation fields as Dozier and Frew (1990), Dubayah (1992Dubayah ( , 1994), etc), etc., these approaches are not commonly included in GIS-based hydrological models, which usually adopt simple approaches to estimate the incident radiation throughout the watershed, as explained next.In AnnAGNPS, a distributed-parameter, physically-based, continuous-simulation, watershed-scale, nonpoint-source pollutant model (Cronshey and Theurer, 1998), correction factors to take account of the effect of dust, water vapor, path length, and reflection and rescattering are applied to the extraterrestrial radiation in order to obtain the short wave radiation received at the ground surface.However, these corrections are simplified into two multiplicative factors, one reflecting the effects of the atmosphere as a function of the elevation and another as the influence of clouds, which depends on sunshine fraction for each day (Bingner and Theurer, 2003).Land area representations of a watershed are used to provide landscape spatial variability, so climate variables remain constant at a subwatershed scale, and, therefore, do not involve topographic factors.SWAT (Arnold et al., 1998) is a lumped model, in which each subwatershed is associated with a unique radiation gauge.Here, topographic corrections are not considered and the measured solar radiation data, when available, are directly applied on the whole region of influence by means of estimated extraterrestrial radiation.MIKE-SHE (Refsgaard and Storm, 1995), is a comprehensive, deterministic, distributed and physicallybased modeling system capable of simulating all major hydrological processes in the land phase of the hydrological cycle (Singh et al., 1999).However, up to now, atmospheric processes have not generally been modeled explicitly, and, whereas precipitation is a direct input in MIKE SHE, radiation and water vapor transport in the atmosphere are typically bound up in evapotranspiration models (Graham and Butts, 2005), and simple methods, such as Thiessen polygons or other areal methods, are usually applied to extrapolate the point scale values for the referred stations at a watershed scale (Singh et al., 1999;Vázquez et al., 2002;Vázquez and Feyen, 2003).
The aim of this study was to address the importance of incorporating in hydrological models the effects of topography on the spatial distribution of global solar radiation at watershed scale.For this purpose, different topographic algorithms have been coupled in order to estimate series of distributed solar radiation values and calculations have been made to quantify this influence on evapotranspiration estimates in mountainous areas in Mediterranean locations.Thus, an algorithm was derived from Dozier (1980), Jacovides et al. (1996) and Ineichen and Pérez (2002) to take into account the lack of weather stations at high altitudes.To be exact, this algorithm should estimate hourly global values as well as the separation between its beam and diffuse components from the common measurements obtained on horizontal surfaces.The resulting algorithm was implemented on a GIS-based routine and applied to data from a mountainous watershed on the south coast of Spain.The distributed results were compared to those obtained from simpler interpolating methods and experiment data.Finally, in order to address the hydrologic importance of using topographically corrected solar radiation fields over uniform values, a simple evaluation in terms of their influence in the computation of reference evapotranspiration fields has been carried out.

Study area and data sources
The study area is the Guadalfeo river watershed, Southern Spain (Fig. 1), where the highest altitudes in Spain can be found (3482 m) with the coastline only 40 km away, in a 1300 km 2 area which results in the interaction between semiarid Mediterranean and alpine climate conditions, with the regular presence of snow (Díaz, 2007;Herrero, 2007;Aguilar, 2008;Millares, 2008).The combination of these altitudinal gradients, together with the large number of vegetation, landforms and soil types, produces a complex mountainous terrain with variable hydrological behavior.The main part of the watershed, in terms of hydrology, is comprised of the southern hillside of Sierra Nevada, where global radiation is high throughout the year due to its aspect and lack of cloud cover, even during winter, despite the cold temperatures and the presence of snow.However, the deep valleys with a characteristic south-facing orientation lead to important differences in the instantaneous global radiation between the east-and west-facing mountain slopes, especially after sunrise and before sunset when these valleys are mainly in the shade.
The meteorological data used in this study consisted of hourly and daily datasets provided by the three stations (Fig. 1) of the Agroclimate Information Network of Andalusia (RIA) available in the watershed: 601, 602 and 603, whose UTM coordinates are shown in Table 1.Measurements were made with a Skye Llandrindod Wells SP1110 pyranometer, with a characteristic range of 0.35 ∼ 1.1 µm.
The topographic input data are represented by a digital elevation model (DEM) with a horizontal resolution of 30 × 30 m and 1 m of vertical precision (Fig. 1).Surface slope and aspect were calculated for each point in the DEM, using the regression plane through the 3 × 3 neighborhood of a given point after Dozier and Frew (1990).
For the evaluation of the algorithm performance, the daily datasets applied were provided by one station of the Andalusian Alert and Phytosanitary Information Network (RAIF) (referred to as 702 in Fig. 1), which measures the variable with a Skye Llandrindod Wells SP1110 pyranometer, with a characteristic range of 0.35 ∼ 1.1 µm, as well as hourly data recorded at a new weather station installed in 2004 in Sierra Nevada by the University of Granada Environmental flow dynamics Research Group at an elevation of 2510 m (referred to as 802 in Fig. 1).Measurements of global radiation at station 802 were made with a Kipp and Zonen SP-Lite pyranometer, with a characteristic range of 0.4 ∼ 1.1 µm.The sensor was placed on a horizontal surface, partially surrounded by higher ground to the north but completely exposed to the south.Even though the measurements at weather stations constitute a point source of information, for the purposes of these studies, those records can be assumed to constitute representative average values for the cell on which they are located (Batllés et al., 2008;Martínez-Durbán et al., 2009).

Calculation of solar radiation components including topographic effects
The topographic effects on solar radiation are mainly variations in the illumination angle and shadowing from local horizons, the apparent intersection of the earth and the sky as seen by an observer in a certain direction.The local horizon information from the gridded data allows us to ascertain whether a given location at a certain sun position is shaded from direct sunlight by surrounding terrain and determines, at any location, the portion of the overlying hemisphere which is obscured by the terrain (Dozier et al., 1981;Dubayah, 1992).Thus, each hourly component, beam, diffuse and reflected radiation has been calculated separately to account for the topographic effects (González-Dugo et al., 2003) from the daily global radiation data measured at weather stations.
According to Essery and Marks (2007), despite the fact that since the availability of gridded data and powerful computers many efficient algorithms for calculating distributions of solar radiation over topographic grids have been developed, all of them implement the same basic geometric principles.Thus, the calculation of horizons in this study was made following the modification to the method by Dozier (1980), made more computationally efficient by Dozier et al. (1981) and Dozier and Frew (1990).They developed a simple and fast algorithm for the extraction of horizons from DEMs by comparing slopes between cells in a certain direction, and formulated the problem by determining the coordinates of the points which constitute the horizons in each cell.Then, by rotating the matrix and solving it in a one-dimensional way along each row as many times as the directions consid-ered, they derived the horizons in the whole hemisphere for each cell in the DEM.In this study, eight directions were considered in the calculation: the four cardinal points and their mid-way points.

Beam and diffuse component estimation on horizontal surfaces
The extraterrestrial radiation at any moment of time incident upon a horizontal surface located at an angle relative to the sun's beams (R o ), is a function of the solar coordinates (the zenith angle (θ z ) or its complementary angle, the sun elevation angle (h s ), and the solar azimuth (ψ)) and the solar constant (I CS ) corrected with the eccentricity factor (E o ) to account for the changes in distance from the earth to the sun along the elliptical trajectory: I CS was fixed as 1367 Wm −2 (Frölich and Brusa, 1981) and solar coordinates, which were calculated following the equations in Iqbal (1983), are functions of geographical latitude (φ L ), and solar declination (δ) in radians.For the solar declination, the value at noon can be used as a constant daily value (Iqbal, 1983) and, as a medium-sized watershed, unique values for latitude and longitude were considered, and, therefore, a constant daily value of extraterrestrial solar radiation was obtained for the whole watershed.
In order to obtain the total amount of global radiation during one day (MJ m −2 day −1 ), extraterrestrial radiation (Eq. 1) must be integrated from sunrise to sunset.Thus, by assuming that the solar beam angle originates from the center of the solar disk, Eq. ( 1). was integrated following the expressions in Iqbal (1983) between the beginning and ending sun-hour angles when the sun's beam first and last strikes the surface (Allen et al., 2006).
The basic procedure for the partition of daily global radiation into its components is the application of empirical correlations between the daily global radiation (R g ) and its diffuse component (R d ), and the ratio of daily global radiation to extraterrestrial radiation (R o ), known as the daily clearness index (CI) (Eq.2).The correlation applied in the present study, shown in Eq. ( 3) was obtained by Jacovides et al. (1996), who investigated the accuracy of some of the previously available correlations when applied locally for high-quality data recorded in Cyprus, and found that these correlations are location-independent.However, they developed a specific correlation, which is more suitable for applying in Mediterranean areas.
Beam daily solar radiation (R b ) can be obtained as the difference between global and diffuse radiation.Therefore, by applying Eq. ( 2), daily R d and R b values (MJ m −2 day −1 ) can be estimated at each weather station with daily R g (MJ m −2 day −1 ) available datasets.However, for the spatial interpolation of daily beam and diffuse solar radiation, the lower scattering and absorption effects in high elevations, due to the lower air mass and vice versa, must be taken into account, especially with a quickly changed altitude (Yang et al., 2010).In this way, an atmospheric correction has to be performed pixel by pixel (Li et al., 2002).
The estimation of clear-sky solar radiation helps to evaluate the effects of the atmospheric properties in a cloudless atmosphere.In this way, incident global daily radiation could be defined as being the product of the incident daily global radiation under clear sky conditions (R gcs ) and a cloudiness factor (f CI cl ) that attenuates R gcs .
Thus, a theoretical clearness index for clear sky conditions (CI cs ) could be defined as the ratio of global radiation under clear sky conditions to extraterrestrial radiation (R gcs /R o ), In the literature there are numerous reviews and evaluations of the different models for the estimation of clear-sky solar radiation (Ineichen, 2006;Annear and Wells, 2007;Tham et al., 2009).In this study, the approximation of Ineichen and Pérez (2002) was applied based on its implementation simplicity.Thus, R gcs can be obtained for each hour of the day from the sun elevation angle, the altitude in meters (z), the Linke turbidity factor (T L ) and the air mass (AM) calculated in terms of the height and solar elevation according to Kasten and Young (1989) and T L for air mass 2. However, even though T L can vary within a year (Mavromatakis and Franghiadakis, 2007), a constant value of 2 was considered in this study as a first approximation.
a 1 , a 2 , f h1 and f h2 are empirical expressions that relate the altitude of the station to the altitude of the atmospheric interactions (Rayleigh and aerosols).A detailed description can be found in Ineichen and Pérez (2002).
Replacing R gcs in Eq. ( 4) by CI cs leads to Eq. ( 6), similar to Suehrcke's formulation (2000): Therefore, the clearness index (Eq.7) can also be expressed as the combination of two factors, CI cs that shows the influence of the atmosphere over solar radiation in a cloudless atmosphere, and f CI cl , that decreases the incoming solar radiation due to cloudiness effects.Eq. ( 7) constitutes the basis for the distributed computation of the CI in order to calculate the spatial fields of the different daily solar radiation components.CI cs can be easily obtained at each cell and at an hourly scale as it varies with the sun elevation angle and the height of the cell (Eq.5), and, then, mean daily distributed values are calculated cell by cell.In addition, the daily CI values are calculated at each weather station from Eq. ( 1) with daily measured data, and from Eq. ( 7) with the mean daily values of CI cs at the cells where the weather stations are located, the mean daily values of f CI cl at each station can be obtained.Thus, f CI cl values are spatially interpolated following the inverse distance weighed (IDW) method, in order to distribute it throughout the watershed.This may appear to be an unrealistic simplification, but it is justified by the lack of more spatially distributed registering sites, which would allow a better assessment of the spatial variation of cloudiness in the watershed.From f CI cl and CI cs fields, daily global radiation values (MJ m −2 day −1 ) can be calculated (Eq.6) at cell scale and finally from Eq. ( 7) daily CI fields become available.Once daily global radiation and daily CI fields are obtained, the application of Eq. ( 3) is straightforward at the cell scale, and so distributed values of the daily diffuse and beam radiation components (MJ m −2 day −1 ) can be computed.
For the temporal downscaling of the solar radiation components, the application of hourly relations between hourly CI and hourly diffuse radiation values was initially considered following previous works in the literature (Orgill and Hollands, 1977;Bugler, 1977;Erbs et al., 1982).However, the aim of this work was to provide a feasible method to include topographic effects on radiation at watershed scale, and the size and heterogeneity of the study site together with the lack of weather stations, which, unfortunately, are usual circumstances in many locations, make it unreasonable to spatially interpolate hourly CI values as the spatial distribution of diffuse and direct radiation shows a better correlation at a daily scale.Thus, a simpler approach is proposed so that once the daily values of each component are obtained for each cell, the hourly values (r b and r d , both in (MJ m −2 h −1 )), are computed by distributing the daily amounts throughout the day following the temporal pattern of extraterrestrial hourly radiation during the day.
Finally, hourly values of beam and diffuse radiation on horizontal surfaces can then be transposed to give hourly radiation on tilted surfaces and the assessment of topographic effects as explained next.

Beam radiation
Under the isotropic assumption, hourly beam radiation on a surface of slope β and orientation γ (r b,βγ (MJ m −2 h −1 )) for a certain hour angle (ω) can be expressed in terms of the hourly beam radiation on an horizontal surface (r b ), the zenith angle (θ z ) and a new corrected zenith angle for the sloping surface (θ) (Iqbal, 1983): Therefore, for the calculation of hourly beam solar radiation on tilted surfaces, a correction in the solar coordinates is necessary, so that the cosine of the zenith angle includes the effect of slope and orientation.This corrected zenith angle or illumination angle (θ), function of the sun-earth-tilted surface geometrical relationship, can be obtained as (Iqbal, 1983;Allen et al., 2006): Therefore, the main factors conditioning the fraction of beam radiation are not only the slope and aspect of the location relative to its neighbors, but also the location of the sun relative to the slope at each time step.A certain location is receiving direct sunlight if none of the following situations are taking place: -Self-shadowing due to its own slope: this takes place if the vector normal to the surface forms an angle greater than 90 • with the solar vector (González-Dugo et al., 2003) (e.g.north-facing hill slope and 45 • slope, sun in the south at 30 • over the horizon).This situation is easy to calculate, as Eq. ( 7) yields a negative value.
-Shading cast by the nearby terrain: in this case, the sun is hidden by a local horizon.This case is more complex, since, unlike slope and orientation, it cannot be calculated with information restricted to the immediate neighborhood of a given point (Dozier et al., 1981).In order to express it mathematically, the term known as horizon angle in a certain direction φ, H , is introduced as the angle between the normal to the surface and the line joining this point or grid in the DEM with another point in the same direction high enough to block solar radiation.Thus, shading by the surrounding terrain will occur for each time step if the illumination angle is greater than the horizon angle in that direction.

Diffuse radiation
Topography influences the diffuse component by modifying the portion of the overlying hemisphere visible at a certain point.The computation of scattered and reflected radiation amounts from the atmosphere to the slopes is rather complicated, owing to their non-isotropic nature.A common assumption made is that the diffuse component of solar radiation (sky light) has an isotropic distribution over the hemispherical sky.However, the non-isotropic character of diffuse radiation fields (maximum intensities near the sun and the horizons, minimum intensities in the direction normal to that of the sun, etc.) makes the simplified assumption sufficiently unrealistic to introduce errors into calculations of the energy incident on sloping surfaces (Temps and Coulson, 1977).Nevertheless, following the ideas of Kondratyev and Manolova (1960), who concluded that the isotropic approximation is sufficient for practical purposes (Klutcher, 1979), the isotropic assumption will prevail in this study with the portion of overlying hemisphere visible at each cell as the main factor controlling this component.Thus, the hourly diffuse radiation (r d,βγ (MJ m −2 h −1 )) on a surface of slope β and orientation γ , is: where the sky view factor, SVF, is the ratio between the diffuse component at one given point and that on an unobstructed horizontal surface, so that it corrects the incoming radiation incident on a flat surface to radiation over a sloping and possibly obstructed surface (Dubayah, 1992).Under the assumption of an isotropic sky, a constant value for the SVF can be expressed analytically in terms of the different horizons in each direction considered, as (Dozier and Frew, 1990):

Reflected radiation
Albedo refers to the global reflectance of the surface to solar radiation.Both albedo and topography can vary over short distances, and their interaction can lead to a wide variability in global solar radiation on a scale of meters (Dubayah, 1992).Reflected radiation (r r,βγ (MJ m −2 h −1 )) can be computed following the ideas of Dozier and Frew (1990) from: where the term in brackets represents the terrain configuration factor for isotropic conditions and infinitely long slope, and ρ is the albedo of the surface.The spatial average of albedo is a factor which is difficult to estimate (Tasumi et al., 2006).In this study, the albedo was estimated by Díaz et al. (2007) from the remote sensing data available from Landsat-5 and Landsat-7 satellites during the study period.After the images had been properly corrected and their reflectivity values extracted, albedo values were obtained at the cell scale through the method proposed by Brest and Goward (1987) and interpolated for the whole time lapse on a daily basis.

Global radiation
Global radiation at an hourly scale at each cell (r g (MJ m −2 h −1 )) was obtained as the sum of each component at an hourly scale once: (1) direct radiation had been corrected by self-shadowing and shadows cast by nearby terrain; (2) diffuse sky radiation included the portion of the overlying hemisphere that could be obstructed by nearby terrain, and (3) direct and diffuse radiation reflected by nearby terrain towards the location of interest had been calculated from both corrected components (Dubayah, 1994).And finally, daily global radiation (R g (MJ m −2 day −1 )) was calculated by integration of the hourly global radiation estimates cell by cell.
For the purposes of computer simulation programs, which handle vast amounts of data, the algorithm was implemented in Matlab during the trials, and, finally, in C++ to obtain a sufficiently fast computation, considering all the processes involved at a cell scale.Also, some of the assumptions that could appear to be quite unrealistic due mainly to the scarcity of data, are intended to achieve a compromise between a sufficiently representative distributed approximation and a highspeed processing algorithm that can be run on a desktop PC, from the comparison with measured data and simpler interpolation techniques.Finally, through the use of daily samples, the availability of data is enhanced as not many hourly records are needed; this allows the use of the algorithm in mountainous areas which lack a high frequency monitoring network, so common in many other areas.

Evaluation of the topographic effects on solar radiation fields on reference evapotranspiration estimates
The choice of a method for the calculation of ET 0 depends on numerous factors.The energy available at the soil surface is the first control of the process, so the estimation of this factor from accessible data sometimes conditions the method (Shuttleworth, 1993).In this study, the ASCE-Penman Monteith equation (Eq.13) was applied for the estimation of evapotranspiration over a reference surface (Allen et al., 1998): where ET ASCE 0 is the reference evapotranspiration during a certain time step (mm/ t); the slope of the vapor pressuretemperature-curve saturation calculated at mean air temperature (kPa/ • C); γ the psychrometric constant (kPa/ • C); R n and G the net radiation (combination of net short-wave and net long-wave radiation) and soil heat fluxes, respectively, both in mm/ t water equivalent; e a and e s the actual and saturation vapor pressure (kPa), respectively; T the daily mean air temperature ( • C) and u 2 the wind speed, both measured at a height of 2 m above the soil surface (m/s).Finally, C d and C n are resistance coefficients which vary with the reference crop, temporal time-step and, in the case of hourly time-steps, with daytime and night time.Here, the reference surface defined in the FAO PM equation was considered so that for daily time steps the values of C d and C n were 900 and 0.34, respectively (Allen et al., 1998;Gavilán et al., 2008).
The calculation of some of the variables involved in the ASCE-PM equation can be found in detail depending on the input data available in Allen et al. (1998).Saxton (1975) found out that the variable to which the equation is most sensitive is net radiation, here expressed as the sum of net short-wave radiation, R ns , and net longwave radiation, R nl , at the Earth's surface.Net daily shortwave radiation (Eq.14) on the soil surface, as the difference between incident and reflected radiation, can be expressed in terms of the albedo of the surface, α (0.23 for the reference surface) and the predicted incoming solar global radiation, R g (MJ m −2 day −1 ).In the same way, the net long-wave radiation was approximated by Eq. ( 15) where ε atm is the atmospheric emissivity, T (K) the mean air temperature and σ Stefan-Boltzmann's constant (4.903 × 10 −9 MJ K −4 m −2 day −1 ).The atmospheric emissivity was calculated through a parametric expression by Herrero et al. (2009) based on near-surface measurements of solar radiation and relative humidity, valid for the local conditions of the study area.Therefore, the expression for net long-wave radiation remains as previously done by other authors (Doorenbos and Pruitt, 1977;Allen, 1986;Allen et al., 1998) as a modification to Stefan-Boltzmann's law due to the absorption and downward radiation from the sky.Thus, the product of the Stefan-Boltzmann's constant and the mean air temperature to the fourth power is multiplied by a cloudiness and an air humidity factor (Allen et al., 1998;Donatelli et al., 2006), both factors included in this study in the term ε atm and so, together with the mean air temperature, they constitute the only inputs to the equation.
Finally, as the algorithm is able to derive global radiation values at the cell scale and once the rest of the inputs to the ASCE-PM equation are also available at the cell scale (Herrero et al., 2007), the influence of topographic effects is evaluated in a distributed manner in ET ASCE 0 , estimated after using the topographically corrected values in comparison with distributed estimates by IDW of the solar radiation data recorded at the weather stations.

Results and discussion
In order to run the proposed set of algorithms at the watershed scale, hourly global radiation was calculated from each 30 × 30 m 2 cell of the DEM in the study area for the period comprised between 4 November 2004 and 2 July 2010.
The results are organized into three sections.Firstly, comparisons of the results obtained through the topographic radiation algorithm previously exposed with those derived from a classical interpolation technique are shown.Secondly, the suitability of the results at different temporal scales is presented through its comparison with field measurements, proving the accuracy of the estimated values for hydrological distributed modeling.Finally, in order to address the hydrologic importance of using topographically corrected solar radiation fields over uniform values obtained through classical interpolation techniques, the influence of both estimations as inputs to reference evapotranspiration computations is evaluated.

Topographic corrections v. classic interpolation techniques on solar radiation estimates
At the first stage, topographic information was derived from the DEM (Fig. 1) in the study area.Slope and orientation maps were obtained and the horizons for each cell were calculated as stated before.Once these parameters are available for a certain area they can be used in subsequent executions as they are considered to be independent of the time of year.
In order to compare the results obtained through the topographic algorithm with those of classic interpolation techniques at an hourly scale, a reference day was selected.This date, 20 November 2004, was chosen as it was cloudless and it had not rained for several days.This condition is very important for the albedo estimation from remote sensing images, as the presence of moisture in the environment influences the quality of the estimates, and, therefore, consecutively dry days are most suitable for an accurate performance.Moreover, remote sensing images were available for this date, and, therefore, the errors due to the temporal interpolation of albedo values were minimized.
The hourly sequence of global radiation, as the sum of each component at an hourly scale once each component has been properly corrected, is shown in Fig. 2a, where the spatial gradient in hourly global radiation is evident.On the whole, it can be seen that the locations receiving more radiation are those in the highest part of the watershed, with a south-facing orientation that remains unobstructed during most of the hours of daylight.
In order to assess the potential of the topographic corrections, a simple interpolation technique was applied for the same date.For this, the inverse distance weighed was applied to the hourly values of global radiation measured in the stations with available hourly datasets (Fig. 2b) and the basic statistics were computed.Thus, absolute minimum and maximum values in the watershed, as well as the mean and standard deviation of the hourly estimates through both distributed computations, are shown in Table 2.It can be seen that even though the mean values in the watershed present a slightly similar order of magnitude, differences in the standard deviation, maximum and minimum values between both methodologies were remarkable in every hour of the day.From contrasting results between Fig. 2a and b, not only was a huge difference visible in the distributed values of the variable cell by cell, but also noticeable was the wider range of global values in the watershed, when topographic factors are taken into account, as shown by the maximum, minimum and standard deviation values in Table 2.In this latter case, extreme values, far exceeding the measured values, represent extreme conditions, such as high areas remaining unobstructed most of the daytime and sometimes receiving almost double the values obtained through interpolation of the data recorded at the stations or, at the other end of the scale, valleys that receive minimal or even zero quantities of solar radiation, due to the configuration of the surrounding terrain.The higher maximum hourly estimates and standard deviation values, as well as the lower minimum hourly estimates obtained by the topographic approximation (Table 2), quantify this fact.As a result, processes such as evaporation or snowmelt, which rely heavily on solar radiation, can be miscalculated under a wide range of conditions, such as overestimations in areas obstructed by nearby terrain or underestimations in the upper and exposed regions of the watershed, among others.
Hourly values can be aggregated in each cell at the required temporal scale.In this way, Fig. 3 represents the spatial distribution of daily global radiation on 20 November 2004 estimated through the topographic algorithm (Fig. 3a) and from IDW (Fig. 3b), respectively.Again, the same conclusions can be drawn as at an hourly time step, since maximum and minimum values found in the watershed considering topographic effects are quite different to those obtained through IDW and the daily values recorded at the weather stations (Table 1).Thus, we found differences of as much as an extra 42% in the estimated daily values compared with those obtained through spatial interpolation without consideration of topography, predominantly on the southfacing hillsides in the northern part of the watershed, and underestimations of up to −1800 % in certain cells obstructed most of the daytime, with a mean difference of −4% in the whole watershed.In addition, the same computation was applied for different cloudless days when there were remote sensing images available (1 January 2005, 30 March 2005, 2 June 2005 and 20 July 2005).Table 3 again shows the same general trend in the basic statistics (absolute minimum and maximum values and the mean and standard deviation of the daily estimates in the watershed) for all the days considered: proximity in the mean values obtained by both methodologies, but considerable differences in the extreme values and in the standard deviation, as stated before.
Finally, the same comparison for the hydrological year 2004, from 1 September 2004 to 31 August 2005 (Fig. 3c Table 4. Maximum, minimum, mean and standard deviation (MJ m −2 year −1 ) of annual solar radiation estimates in the watershed obtained by the topographic approximation and IDW to the annual accumulated values measured at stations with daily available datasets (601, 602, 603, 702 and 802).and d) resulted in a mean excess of 323.75 MJ m −2 year −1 with a standard deviation of 860 MJ m −2 year −1 when applying IDW over the topographic computation and extreme differences of the same order of magnitude as at the daily time step.The same comparison for the rest of the hydrological years of the study period can be seen in Table 4, where the last year is incomplete as data from July and August 2010 were not yet available.Once again, higher extremes and standard deviation values through the topographic approximation were obtained.However, unlike at hourly and daily time scales, the accumulation of differences at annual scale leads to a certain variability in the mean values of annual radiation estimates in the watershed between both methodologies for every hydrological year (200-300 MJ m −2 year −1 ).

Validation of topographic corrections
The radiation values generated from datasets of stations 601, 602 and 603 were compared to the radiation measurements at stations 802 and 702 and the agreement between generated and measured data was evaluated through 1:1 lines.The period considered in the evaluation was determined by the availability of data in the weather station 802, which included almost six hydrological years, from 4 November 2004 to 2 July 2010.However, even though station 702 presents sparse measurements from the beginning 2008 until nowadays, the adjustment for the rest of the period is shown.
As regards daily values, a close agreement between generated (R gp ) and measured (R go ) values can be observed in Fig. 4 as points tend to line up around the 1:1 line.A very slight underestimation in the generated values can also be appreciated, with the topographic approximation underestimating the values measured at both stations 702 and 802 by 3 and 2%, respectively.At station 802, these underestimations take place especially in summer periods of years receiving more radiation than the average (2006 and 2007) as depicted in Fig. 5.However, this slight disagreement between both measured and generated values appears to be more related to the direct interpolation of the cloudiness factor and to the fact that in those periods the availability of remote sensing images for an accurate estimation of albedo was more limited than to the assumptions of the topographic approximation applied in this study.Also, the low RMSE values as well as the proximity of R 2 estimates to 1 at both stations demonstrate the ability of the topographic approximation to estimate solar radiation in high elevations even when the only available datasets come from weather stations that are located at a relatively low elevation compared to the mean height of the watershed.To conclude, the accuracy of predicted hourly values was assessed in station 802, where measurements at this time scale were available.Despite the scattering effect observed in Fig. 6, which shows the agreement between predicted (r gp ) and measured (r go ) hourly values for the evaluation period, we can say that the algorithm reasonably predicts the observed data with a R 2 of 0.87, especially considering the time scale and some of the assumptions of the algorithm, which at this time step might appear to be rather simplistic.Therefore, the installation of a denser monitoring network provided with solar devices recording hourly direct and diffuse radiation data may improve the results.Firstly, it would provide the spatial scheme required for the spatial interpolation of hourly values.Secondly, it would allow the inclusion of more factors for the spatial distribution of the f CI cl as previously suggested.Finally, the derivation of hourly correlations between the hourly diffuse radiation and the CI would be reasonable in the study area.To sum up, the possibility of working at finer scales would be the ideal as the geometrical relationships involved in the calculation of extraterrestrial radiation are continuous in time.However, independently of this continuous nature of extraterrestrial radiation, the time scale of the computation of the incoming solar radiation is determined by the temporal frequency of the monitoring network.Nevertheless, the calculations with aggregated hourly values at higher temporal scales such as the daily time step showed the same degree of detail as at the hourly time scale (R 2 around 0.9).
Finally, since cloudless skies are required for an accurate characterization of the albedo from remote sensing images, the results at an hourly time step were analyzed considering this effect.Thus, two different atmospheric situations in terms of the occurrence of rainfall are defined as an indicator of the cloudiness in the atmosphere: events (when it rains somewhere in the watershed) and non-events (periods between events).Figure 7 represents hourly values for event days on the left-hand side (a, b, c) and non-events on the right (d, e, f).
Table 5 shows different linear fits for each day represented in Fig. 7 and its calculated R 2 values.As was expected, the predicted values were much better for cloudless skies or non-events when acceptable R 2 values were obtained even when forcing the adjustment to reach the origin.However, as with daily values, the algorithm slightly underestimated the  observed hourly values, which could be improved with the consideration of the variation in cloudiness at an hourly time step and a more detailed study of the temporal variation in the albedo in the area.In any case, the worst fits obtained in situations when events occur are expected to be more closely related to the separation of the different components in the global radiation value than to the topographic interpolation process.
Nevertheless, these results are considered to be acceptable in the framework of the present study as the estimation of global radiation in semiarid environments is especially important for cloudless days, when evaporation processes and snowmelt must be considered in water and energy balance modeling.This is especially true considering that cloudless  days constitute a higher rate than cloudy days associated with situations when events occur in a Mediterranean area like the present study site: in this case, around 75% of clear sky days for the evaluation period.

Influence of the inclusion of topographic corrections on hydrological variables: ET 0
Finally, a distributed computation of ET 0 was applied for the same reference day of Sect.3.1 (Fig. 8a) and the hydrological year 2004 (Fig. 8c) (1 September 2004 to 31 August 2005) once all the variables involved in the ASCE-PM equation had been spatially derived including topographic effects.Also, the ASCE-PM equation was computed with solar radiation surfaces obtained through IDW of the data recorded at the stations as inputs to the equation (Fig. 8b and d) in order to prove the importance of solar radiation fields which include topographic corrections.Again, in Fig. 8a and c, not only the apparent spatial variability of ET 0 estimates cell by cell which follows the topographic gradient can be seen, but also a wider range of values in the watershed than with IDWinterpolated solar radiation fields as inputs (Fig. 8b and d).
Besides, in this latter case, ET 0 estimates in the watershed appear to be more influenced by the spatial distribution of other variables than by solar radiation (e.g.temperature in Fig. 8b).This spatial variability is quantified in terms of the basic statistics of the distributed estimates through both methodologies for the rest of the reference days and hydrological years, respectively (Table 6 and 7).Again, mean values of the same order of magnitude at daily scale but very extreme values for the maximum and minimum estimates when using the topographic corrected solar radiation, and, therefore, higher standard deviation values than with IDW-interpolated solar radiation fields as inputs (Table 6) can be noted.However, as annual estimates of ET 0 are obtained by integration of daily estimates of ET 0 , the  accumulation of the differences between both methodologies at finer scales becomes apparent in the mean annual values (around 60 mm/year) at watershed scale from both methodologies (Table 7).Finally, Fig. 9 represents the mean daily difference at watershed scale for each hydrological year of ET 0 estimations when using IDW-interpolated and topographically-corrected solar radiation fields.It can be seen that there is a general overestimation with IDW-interpolated solar radiation fields every day of the year.Also, these overestimations follow the temporal pattern of solar radiation, and so the biggest differences (0.4 mm/day) take place in summer, when the importance of ET 0 in the water balance is greater.Considering the mean statistics of the difference between both computations on an annual basis, a mean excess of 61 mm/year and a standard deviation of 142 mm/year in ET 0 estimations when using IDW-interpolated solar radiation fields were obtained for the study period.These differences in an area where the mean annual rainfall varies from 450 mm/year on the coast to 800 mm/year on the highest peaks may constitute a considerable source of error in the water balance when applying distributed hydrological models for the management and planning of water resources.

Conclusions
Difficulties are sometimes encountered in utilizing available solar radiation data, since they consist primarily of total (direct plus diffuse) radiation only, and a knowledge of the values for each component is often required, especially for the consideration of topographic effects as they affect each component differently.
Thus, detailed time-series radiation surfaces have been developed, using limited data and relatively simple methods, to run distributed energy and water balance models in www.hydrol-earth-syst-sci.net/14/2479/2010/ Hydrol.Earth Syst.Sci., 14, 2479-2494, 2010  4 2006-2007 17.9 1317.6 992.7 137.1 917.3 1221.8 1053.6 49.9 2007-2008 51.4 1360.1 1036.2 143.2 981.1 1245.6 1096.8 46.8 2008-2009 78.8 1325.7 1032.3 131.5 939.4 1236.2 1091.6 44.2 2009-2010 42.8 943.6 702.1 103.4 632.1 864.3 747.6 32.9 mountainous Mediterranean environments.The interpolation is managed through linear interpolation of the cloudiness effects and the variation of the CI with elevation as a clue to mean daily radiation, plus topographic properties geometrically related to the sun's position at hourly intervals.These calculations are easy to reproduce from standard weather station datasets.The significant incidence of topography on the values of global solar radiation throughout the watershed has been demonstrated by the results of the topographic solar radiation algorithm proposed.In this sense, differences of as much as an extra 42% in the estimated daily values compared with those obtained through spatial interpolation without consideration of topography, and underestimations of up to −1800% in certain cells obstructed most of the daytime were found.This affects the modeling of the slow but thorough drying out of the watershed during periods between events and the modeling of the snowmelt in the highest areas, among other processes.
The simulated results fit the measured values of global radiation at the 2510 m high monitoring point established for this work really well, with a correlation of 0.93 for daily values, and a RMSE of 2 MJ m −2 day −1 , that confirms the validity of the assumptions made in the algorithm, as the previous paragraphs have justified.However, the simulated results constitute a further approach to the accurate characterization of the spatial distribution of hourly global radiation values in mountainous areas with scant data recording sites.On-going work will develop a further approach, and test the inclusion of additional corrective terms in the cloudiness factor through the establishment of two additional weather stations equipped with pyranometers at points with an increasing height above sea-level and distance from the sea.
The importance of considering the topographic gradients in the spatial distribution of solar radiation for the study of hydrological processes in which this variable plays a crucial role became evident against ET 0 estimates with solar radiation fields obtained through classical interpolation techniques of data recorded at weather stations.Thus, a mean excess of 61 mm/year was found with IDW-interpolated solar radiation fields as inputs to the ASCE-PM equation with the highest overestimations taking place in summer periods.

Table 1 .
UTM coordinates of the weather stations, measured daily global radiation and clearness index on the 20 November 2004.

Table 2 .
Maximum, minimum, mean and standard deviation (MJ m −2 h −1 ) of hourly solar radiation estimates in the watershed on the 20 November 2004 obtained by the topographic approximation and IDW to the values measured at stations with hourly available datasets(601, 602, 603 and 802).

Table 5 .
Linear fits of observed (r go ), and predicted (r gp ) global hourly radiation (MJ m −2 h −1 ) at station 802 for certain dates.

Table 6 .
Maximum , minimum, mean and standard deviation (mm/day) of daily evapotranspiration estimates in the watershed on selected days with global radiation topographically corrected and IDW interpolated.Topographic (mm/day) IDW (mm/day)

Table 7 .
Maximum, minimum, mean and standard deviation (mm/year) of annual evapotranspiration estimates in the watershed on selected days with global radiation topographically corrected and IDW interpolated