Regionalization of the Hargreaves coefficient to estimate long-term reference evapotranspiration series in SE Spain

This study employs a methodological approach for estimating long-term series of monthly reference evapotranspiration (ET o ) from historical data. To carry it out a regionally calibrated version of the Hargreaves equation was applied at old ordinary weather stations which only provide data of air temperature and precipitation. The proposed approach was based on the analysis of: (1) the Hargreaves coefficient obtained by local calibration from data of 66 modern automatic weather stations; (2) the regional characterization of the spatial variability of that coefficient by means of a “regional function”; and (3) the final application of this function to the old ordinary weather stations. This approach was assessed under the semiarid conditions of the Segura River Basin (south-eastern Spain) by comparing ET o estimates against those obtained with the Penman-Monteith method, which was used as reference. Spatial variability of the Hargreaves coefficient was well correlated with the annual and monthly means of daily temperature range, so they were selected as explanatory variables for the regionalization of the Hargreaves coefficient following two approaches: a global regional function and monthly regional functions. The regionally calibrated version of the Hargreaves equation by monthly functions clearly improved the performance of its original parameterization (average relative error decreased from 19.8% to 10.1%) although, as expected, estimates were not as good as those obtained with the local calibration (average relative error = 7.7%).

. It is a physicallybased method that requires data for a large number of meteorological variables which are not often all available. In those situations, the Hargreaves (HG) equation is one of the most popular temperature-based methods that provides reasonable ET o estimates with a global validity (Allen et al., 1998).
The knowledge of long-term series of ET o is currently of great interest for studying and modelling hydrological and agricultural systems under dynamic scenarios related to land use and/or climate change (Elgaali et al., 2007;Sánchez-Toribio et al., 2010;Bae et al., 2011;Espadafor et al., 2011;Rey et al., 2011). However, determining long-term series of ET o from historical data presents a major drawback: the lack of reliable records for long periods due to the progressive changes in measuring devices and sitting of weather stations throughout the 20 th century. In addition, applying the PM method is even more limited since weather stations with the required data for long periods are still very scarce (Droogers & Allen, 2002;Pereira & Pruitt, 2004;Simolo et al., 2010). For instance, in the Segura River Basin (SRB), a semiarid region located in southeastern Spain, there are only a handful of meteorological stations with more than 30 years of records with the required data for applying the PM method. Besides, those records are often incomplete and not very reliable. However, long records (> 50 years) of precipitation and air temperature are available at about 50 old ordinary weather stations (OWSs), some of which even have records going back more than 90 years. The international standardisation of the PM method has promoted the implementation of new agrometeorological networks in the SRB since 2000, consisting of modern automatic weather stations (AWSs) specifically equipped for applying the PM method. However, they were mainly located in irrigation districts and far from the site of the old OWSs, not making it possible to relate the historical records of the old OWSs with the new readings at the AWSs.
Given the difficulty of applying the PM method for calculating long-term ET o series in the SRB (similar to other Spanish and worldwide regions), other alternative low data demanding methods must be used. The HG equation (Hargreaves & Samani, 1985) is an appealing method for estimating ET o at meteorological stations when ordinary weather data are available, as it only requires air temperature data. This equation can be extremely useful for both determining long-term series of ET o from historical records, as well as for generating projected ET o series from monthly temperature projections provided by climate change models (Milzow et al., 2010).
Previous works have demonstrated that, in general, the HG equation can provide precise estimations for weekly or longer time predictions (Hargreaves, 1994;Hargreaves & Allen, 2003;López-Urrea et al., 2006a). Moreover, other scientists such as Shuttleworth (1993) recommend that the HG method should not be used for shorter periods than 1 month. The original parameterization of the HG equation (Hargreaves & Samani, 1985) usually overestimates ET o in humid regions, and underestimates it in dry areas (Saeed, 1986;Jensen et al., 1990;Amatya et al., 1995;Temesgen et al., 2005). It has also been reported that the HG equation overestimates ET o at low evapotranspiration rates, and vice versa (Droogers & Allen, 2002;Xu & Singh, 2002;Itenfisu et al., 2003). Those reports make it clear that the HG equation performance is strongly influenced by the climatic conditions where it was parameterised. Therefore, such equation must be evaluated and, if it is necessary, calibrated for accurate use in other zones (Jensen et al., 1997;Gavilán et al., 2006;Shahidian et al., 2013).
An appropriate re-parameterization of the HG equation must be carried out by local calibration of the HG coefficient. The adjustment of the HG coefficient has usually been carried out twofold: by comparison with weighing lysimeter measurements of the reference crop (Jensen et al., 1997;Martínez-Cob & Tejero-Juste, 2004;López-Urrea et al., 2006b), or more frequently by comparison against the ET o estimations provided by the application of the PM method at the same weather station (Itenfisu et al., 2003;Martínez-Cob & Tejero-Juste, 2004;Gavilán et al., 2006;Jabloun & Sahli, 2008;Sentelhas et al., 2010;Espadafor et al, 2011;Mendicino & Senatore, 2013). However, when the HG equation is applied at weather stations where it cannot be locally calibrated, as occurs at the old OWSs of the SRB, other approaches for its regional calibration should be considered.
Several solutions of different complexity were proposed for the regionalization of the HG equation, most of them adopting a regression based calibration of the HG coeff icient using auxiliary parameters such as temperature range (Samani, 2000;Vanderlinden et al., 2004;Lee, 2010;Mendicino & Senatore, 2013), relative humidity (Hargreaves & Allen, 2003), wind speed (Jensen et al., 1997;Martínez-Cob & Tejero-Juste, 2004), presence of large water bodies or distance to coast (Vanderlinden et al., 2004;Mendicino & Senatore, 2013), and rainfall (Droogers & Allen, 2002). Shahidian et al. (2013) carried out an in-depth analysis of the parameters previously used for the spatial and seasonal calibration of the HG method, concluding that it is possible to improve the precision of the estimates for new sites where no reliable records of climatic data exist by using regional averages of such parameters.
The purpose of this study was to develop a methodological approach for obtaining long-term series of monthly ET o by applying the HG equation, with a regionally calibrated HG coefficient, in old OWSs, where historical air temperature data are available. In contrast to previous analysis of the HG equation, monthly-averaged weather data were used to estimate monthly-averaged ET o. This approach was selected since a monthly time step is usual in long-term hydrological and agricultural modelling, especially if future projections of weather data are used. The proposed methodology was based on the estimation of locally calibrated HG coefficients at a set of modern AWSs; the regionalization of the HG coefficient based on the formulation of a suitable function linking it to available information at the OWSs; and the final application of this function at the OWSs. The approach was evaluated under the semiarid conditions of the SRB.

Study area and weather data
The Segura River Basin (SRB) is characterised by a Mediterranean semi-arid climate with warm, dry summers and mild winter conditions. The average annual temperature is 17.5°C, reaching maximum temperatures of 38ºC in summer, and minimum temperatures of 1°C in winter. Annual rainfall averaged 350 mm during the study period, with high seasonal and inter-annual variability. Most precipitation occurred during the fall and winter months, but inter-annual droughts were also common.
Data for an eight-year study period (2001-2008) from two sets of weather stations were used in the study. On the one hand, 66 modern AWSs with available daily data of air temperature, relative humidity, wind speed, and global solar radiation over the study period were selected. These AWSs belong to three different weather and agro-meteorological services: 38 stations are managed by the Agricultural Information Service of Murcia Region (SIAM, http://siam. imida.es); 16 stations are part of the National Agroclimatic Information System for Irrigation (SIAR, http://www.mapa.es/siar); and the remaining 12 stations pertain to the National Meteorology Agency (AEMET, http://www.aemet.es). Most stations in this set were remotely monitored and specifically equipped for calculating ET o with the PM method. Air temperature and relative humidity were measured from 1.5 to 2.0 m above soil surface, whereas wind speed was usually recorded at 2.0 m height. Some of the AEMET stations measured wind speed at higher heights and then their data were adjusted to 2.0 m height by means of a logarithm wind profile equation. The reader is referred to the aforementioned web pages for detailed information about the sensor type and model for recording each meteorological variable at the AWSs. On the other hand, a set of 77 old OWSs were available. These provided historical long-term series of daily maximum and minimum air temperature, and precipitation. The oldest of them have records from 1913 and is even today still in operation. The OWSs belong to the National Meteorological Agency and are graphed as triangles in Fig. 1. They are equipped with traditional analogical thermometers and rain gauges. The three networks (SIAM, SIAR and AEMET) are responsible for quality control procedures, including sensor calibration and data validation. Since showing and discussing the results for all the stations would be unfeasible in a scientific paper, twelve stations of the AWSs were chosen for displaying the results of the local calibration and for assessing the regional calibration performance. A summary of the main features of these stations is provided in Table 1.

ET o estimates
For both sets of weather stations, monthly-averaged data throughout the study period were derived from the daily observations, assuming that a monthly average value was valid if at least 25 daily observations were available. These monthly values of the meteorological variables were used to compute monthlyaveraged ET o using the PM and HG equations. The value of the ET o calculated with monthly-averaged weather data is indeed very similar to the average of the daily ET o values calculated with daily weather data for that month (Allen et al., 1998).
The FAO-56 version of the PM method is considered the most precise and standard method to estimate ET o (Allen et al., 1998), such as it was corroborated in the study region (López-Urrea et al., 2006a). It refers the ET o concept to the rate of evapotranspiration from an extensive area with an ideal 0.12 m high crop with a fixed surface resistance of 70 s m -1 and an albedo of 0.23, that is well provided with water and nutrients. The PM equation is as follows: where ET o is the reference evapotranspiration (mm day -1 ); Δ is the slope of the saturated vapour pressure curve (kPa °C -1 ); R n is the net radiation (MJ m -2 day -1 ); G is the soil heat flux density (MJ m -2 day -1 ); T m is the mean air temperature (°C) at 2.0 m; U 2 is the average wind speed at 2.0 m height (m s -1 ); e s is the saturation vapour pressure (kPa) at temperature T m ; e a is the actual vapour pressure (kPa); (e s -e a ) is the vapour pressure deficit (kPa); and γ is the psychrometric constant (kPa°C -1 ). The soil heat flux (G) was assumed to be negligible over the calculation time step period (1 month). Eq. [1] was used in this study, adopting the procedure suggested by Allen et al. (1998) for calculating monthly-averaged ET o , starting from the monthly-averaged values of temperature, solar radiation, air relative humidity, and wind speed.
As aforementioned, the Hargreaves (HG) equation requires only maximum and minimum air temperature, as well as extraterrestrial radiation. The monthly-averaged ET o (mm day -1 ) were calculated using the following equation (Hargreaves, 1994): where C refers to the HG coefficient, which value is 0.0023 according to the original parameterization proposed by Hargreaves & Samani (1985); R a is the water equivalent of the monthly-averaged daily extraterrestrial radiation (mm day -1 ), calculated according to Allen et al. (1998); T max , and T min are the monthly-averaged maximum and minimum values of daily air temperature (°C); and T is the monthlyaveraged daily temperature, calculated as the average of T max and T min .

Calibrations of the Hargreaves equation
Monthly-averaged ET o for the study period was calculated with the PM and HG methods (ET o,PM and ET o,HG , respectively) at the AWSs. ET o,HG estimations were assessed by comparison against ET o,PM , which was selected as a standard reference for ET o . The assessment at each station entailed analysing the relationship between the reference and the estimated values of ET o by means of linear regressions y = a + b·x, where y = ET o,PM and x = ET o,HG . Additionally, the performance of the HG equation was evaluated using the mean bias error (MBE), the root mean square error (RMSE) and the relative error (RE) in accordance with the following expressions (Willmott, 1982): where n is the sample size (96 months); and ET - After assessing the performance of ET o,HG , HG equation was re-parameterized by local calibration of the Hargreaves coefficient C at each AWS. This local calibration was performed at global and monthly scales.

Local calibration of the Hargreaves equation
The first approach for locally calibrating the HG equation consisted of the estimation of a globally calibrated HG coefficient (C g ) for each single available AWS, which minimizes the differences with the estimates provided by the PM method (ET o,PM ). This was achieved by using a simple optimization technique aimed at minimizing the MBE at each AWS. The monthly-averaged values of ET o obtained with the global local calibration of HG coefficient are referred to as ET o,HGg. As a result, a set of 66 C g values for the SRB were obtained.
The other approach was the fitting of twelve monthly calibrated HG coefficients (C m,j ) for each single available AWS. C m,j minimizes the MBE with the estimates provided by the PM method for month j at each AWS: where n is the sample size and j denotes the month of the year (j=1,…, 12). The monthly-averaged values of ET o obtained with the monthly local calibration of HG coefficient are referred to as ET o,HGm . As a result, 792 C m,j values, one for each month of the year at each AWS were obtained.
The performance assessment for both calibrations proposals was carried out by comparison with ET o,PM , in the same way as the original parameterization of the HG equation.

Regional estimation of the Hargreaves coefficient
In the SRB, the OWSs were located in different sites than the modern AWSs considered for the local analysis and hence, the local calibration of the HG coefficient alone was not enough for this study purpose. Nevertheless, it would be easy to solve this problem if a relationship between the calibrated HG coefficients (C g or C m,j ) and other data also available at the OWSs could be established, such as elevation, the mean daily temperature, the mean daily temperature range, or the mean precipitation. Therefore, the regionalization of the HG coefficient was based on the formulation of a suitable mathematical equation (regional function) linking it to the aforementioned data. This approach was similar to that followed in other regions for previous regionalization of the HG coefficient attempts (Samani, 2000;Vanderlinden et al., 2004;Lee, 2010;Mendicino & Senatore, 2013;Shahidian et al., 2013).

Estimation of long-term ET o series with Hargreaves equation in SE Spain
The regional function of the HG coefficient was incorporated to Eq. [2] for obtaining the regionally calibrated version of the HG equation, which was validated in terms of ET o at the twelve selected AWSs using cross-validation. This procedure consisted of (i) recalculating the regional function removing each time one of the AWS from the analysis, and (ii) assessing for the removed AWSs the ET o estimates achieved with that regional function by comparison with ET o,PM . Finally, the regionally calibrated version of the HG equation was applied to the set of OWSs, from which historical time series of monthly-averaged ET o were obtained from ancillary records of air temperature.

Performance of the original parameterization of the Hargreaves equation
The results of the linear regressions between the values of ET o,PM and ET o,HG at the 12 AWSs chosen for displaying the results are shown in Table 2. This table also shows the average values for all the AWSs. The coefficients of determination were always over 0.96, thus representing a high and steady correlation between both ET o estimation methods. The slope of the adjusted linear functions, b, was > 1 for all the stations except for Almoradí, where it was 0.948. On average, b was 13.5% higher than the unity, indicating a systematic ET o underestimation of the original parameterization of the HG equation. It should be noted that the regression slopes were statistically different than 1 (Ttest, α = 0.95) for most AWSs. The analysis of the performance statistics (MBE, RMSE, and RE) showed substantial errors and a great variability of results among stations. Considering all locations, the MBE values ranged from 0.397 to 1.132 mm day -1 , with a mean of 0.741 ± 0.268 mm day -1 , RMSE ranged from 0.471 to 1.210 mm day -1 , with a mean of 0.862 ± 0.264 mm day -1 and RE ranged from 12.03 to 26.98%, with a mean of 19.78 ± 5.05%. The positive values of MBE again clearly showed the systematic underestimation of ET o,HG . The comparison of these statistics with the reported by Gavilán et al. (2006) in the nearby region of Andalusia indicated that although the average RMSE and RE were almost the same (similar scatter with respect to the reference method), the MBE was higher and with opposite sign (very different bias with respect to the reference method). This MBE behaviour was expected since the SRB is considerably arider than Andalusia.   Fig. 3 for Murcia AWS. The underestimation of the HG equation was systematic, but higher in summer than in winter, in agreement with previous reports indicating the increasing underestimation of the HG equation with increasing ET o rates (Amatya et al., 1995;Xu & Singh, 2002).

Performance of the Hargreaves equation locally calibrated by a global coefficient
The results of the linear regressions between the monthly values of ET o,HGg and ET o,PM are shown in Table 3.  In this case, the regression slopes were not statistically different than 1 (T-test, α = 0.95) for most AWSs. Moreover, the underestimation in the monthly evolution of the ET o,HG disappeared in winter and was reduced in summer, as can be observed when depicting the monthly evolution of ET o,HGg (Fig. 3). The comparison between ET o,HGg and ET o,PM for Murcia and Alcantarilla AWSs (Fig. 2) evidenced an improvement in the estimations with respect to ET o,HG , although a slight underestimation for high ET o rates (summer months) was still observed. The F-test (α = 95%) stated that the regression slopes and intercepts were statistically significantly different than the provided with the original HG equation for both AWSs. Considering the 66 AWSs, the MBE values ranged from -0.031 to 0.089 mm day -1 , with a mean of 0.007 ± 0.028 mm day -1 , RMSE ranged from 0.285 to 0.615 mm day -1 , with a mean of 0.412 ± 0.095 mm day -1 and RE ranged from 7.02 to 14.14%, with a mean of 9.61 ± 2.13 %. It should be noted that the RE decreased on average about 51% with respect to ET o,HG .
The globally calibrated HG coefficient (C g ) values were higher than the proposed in the original HG equation (C = 0.0023) for all the stations with the exception of Almoradí, where it was 0.00227. C g ranged between 0.00227 and 0.00362, with an average value of 0.00285 ± 0.00035. Fig. 4 shows the spatial variation of C g in the SRB, calculated from its values at the 66 AWSs and by applying an inverse distance weighting interpolation method (Watson & Philip, 1985). An important regional variation of C g was observed, with a marked decreasing gradient from the coast to inland, inversely related to the altitude va-

Performance of the Hargreaves equation locally calibrated by monthly coefficients
The results of the linear regressions between ET o,HGm and ET o,PM , as well as the performance statistics MBE, RMSE, and RE are included in Table 4. As expected, the ET o,HGm values fitted with ET o,PM even better than ET o,HGg in all studied stations. The slope of the regression ranged from 0.99 to 1.01, thus indicating that the slight underestimation observed with ET o,HGg was eventually corrected by the monthly local calibration. The negligible MBE value at all stations confirmed this circumstance. The monthly trend of ET o,HGm and ET o,PM (Fig. 3) also depicted very small differences between them, without any clear systematic error for both winter and summer months. Additionally, the ET o,HGm versus ET o,PM relationship for Murcia and Alcantarilla stations (Fig. 2) was very near the straight line 1:1. The F-test (α = 95%) also stated that the regression slopes and intercepts were significantly different than the provided by ET o,HGg for both stations.
The RMSE of ET o,HGm ranged from 0.217 to 0.443 mm day -1 , with a mean of 0.338 ± 0.071 mm day -1 , and RE ranged from 5.42 to 10.68%, with a mean of 7.71 ± 1.67%. These values entailed an average RE decrease of 61.24% and 19.77% with respect to ET o,HG and ET o,HGg , respectively. These results evidence than the monthly local calibration, which was not considered in previous similar studies, performs better than the global local calibration.
The monthly evolution of the average value and the standard deviation for the 66 AWSs of the monthly calibrated HG coefficient (C m,j ) is displayed in Fig. 5. The average C m,j presented an increasing trend during the spring, reaching its maximum value in July (C m,July = 0.00312). Subsequently, it decreased steadily throughout the autumn and showed its minimum value in January (C m,January = 0.00273). C m,j pattern was similar to that observed for the average monthly temperature (T m , secondary axis in Fig. 5), in correspondence with the reported overestimation and underestimation of the HG equation for low and high ET o rates (Xu & Singh, 2002), since higher C m,j were observed in the hottest months to correct the higher underestimation of the HG equation. It should be noted that all average C m,j values were above 0.0023, indicating again that the original parameterization of the HG equation underestimated it throughout the whole year (Fig. 3). The link between C m,j and T m also led to higher variations of C m,j at inland stations than at those located near the coast, due to the moderating effect of the sea. Moreover, the standard deviation of C m,j was higher in winter than in summer, similar again to T m behaviour.

Regional calibration of the Hargreaves equation
Extrapolation of locally calibrated HG coefficients to the old OWSs location would be possible if a relationship between C g or C m,j and another parameters measured at those OWSs could be established. Most recent HG coefficient regionalization attempts were proposed by Samani (2000) in USA; Vanderlinden et al. (2004) and Gavilán et al. (2004) in Andalusia; Lee (2010) in Korea; and Mendicino & Senatore (2013) in southern Italy. These authors found reliable formulations for C g based on several parameters such as the annual mean of daily temperature (T -), the annual mean of daily temperature range (ΔT -) or the ratio , all of them being available at the OWSs. Bearing in mind the results of these precedent studies, an overall correlation analysis including C g , T -,

ΔT
-, , and other available parameters at the OWSs, such as elevation and the mean annual precipitation (P -), is shown in Table 5. It can be seen that C g was highly correlated with some temperature related parameters, such as ΔT and , whereas the correlation with the other parameters was not statistically significant (pvalue > 0.05, α = 95%). An equivalent correlation analysis including C m,j ; the monthly mean of daily temperature (T j -); the monthly mean of daily temperature range (ΔT i -); ; the mean monthly precipitation (P j -); and elevation also confirmed the high   correlation with ΔT j and at monthly scale. Specifically, the correlations of C g and C m,j with ΔT and ΔT j -, were high in both cases: -0.78 for C g , and ranging from -0.62 in January to -0.82 in June for C m,j . These high correlation values indicated that linear equations would be sufficient to explain most of the regional variability of the HG coefficient. Consequently, the following relationships were found between C g and ΔT -, and between C m,j and ΔT j -, at each AWS: [7] [8] where a, b, a j and b j are adjustment parameters and j denotes the month of the year (j = 1,…, 12). Opposite to other authors (Samani, 2000;Vanderlinden et al., 2004;Mendicino & Senatore, 2013) quadratic or power functions between C g and ΔT -, or between C m,j and ΔT j -, did not explain the variability of the HG coefficient significantly better than the linear functions. Fig. 6 graphically represents the relationship between C g and ΔT -, along with the regression curves proposed by Samani (2000), Vanderlinden et al. (2004) and Mendicino & Senatore (2013). It can be observed that none of the previous proposals fitted properly with C g in the study region. The relationship proposed by Vanderlinden et al. (2004) in Andalusia was the best approximation, as it was expected due to the vicinity of both studied regions, but in any case it showed a satisfactory estimation. This result highlights regional functions for the HG coefficient cannot be extrapolated to other regions, even in their vicinity. Table 6 shows the parameters value, their standard errors, the correlation coefficient and the p-value for Eqs.
[7] and [8]. The coefficients a and b for the global regional function were equal to -0.000132 and 0.00442, respectively. Taking into account that only one parameter was considered in Eq. [7], the value of the coefficient of determination (R 2 = 0.61) was slightly  Eq.
[7] and [8]. SE: standard error of the function. SE i : standard error of parameter i.
high in comparison with the aforementioned precedent regionalization attempts, probably due to the smaller size and uniform climatology in the study area. The slopes of the monthly regional functions, a j , decreased from January (a January = -0.000267) to July (a July = -0.000095), following an opposite behavior than the solar radiation. The interception, b j , showed a similar trend, with the maximum value in December (b December = 0.00674) and the minimum in July (b July 0.00396). The coefficient of determination presented similar values than for C g throughout the year, with the exception of the winter months (January to March), when it presented lower values (0.41 to 0.47).
The regional estimation of the HG coefficient by means of a global regional function or monthly regional functions were validated by comparison with ET o,PM at the twelve selected AWSs using crossvalidation. Table 7 shows the statistics (MBE, RMSE, and RE) for both the global (Eq. [7]) and the monthly (Eq. [8]) regional functions. Fig. 7 compares the RE value after both regionalization approaches with the original parameterization of the HG equation (ET o,HG ) and with the local calibrations (ET o,HGg and ET o,HGm ). The global regional function for the HG coefficient performed better than ET o,HG , reaching lower values of MBE, RMSE and RE at the twelve AWSs. The average RE value achieved was 12.5%, leading to a reduction of the error with respect the ET o,HG of 36.9%. The average RE increased only a 30.2% with respect to the global local calibrations (ET o,HGg ). The monthly regional functions for the HG coefficient performed even better than the global regional function, improving the results for the twelve AWSs. The average RE value achieved was 10.1%, leading to a reduction of 48.9% with respect to ET o,HG , and to a rise of 31.2% with respect to ET o,HGm . Moreover, the monthly regional functions performed better than the global local calibration of the HG coefficient at two location (Abanilla and Beniel), showing similar average RMSE and RE values. This circumstance is specially relevant for our study, as this implies that the regionalization by monthly regional functions can perform almost as well as a global local calibration, which usually was the followed approach in previous regionalization attempts.
Finally, long-term series of ET o were estimated at the OWSs by combining Eqs.
[2] and [8], since it clearly performed better than Eq. [7]. As an illustration for this practical application, Fig. 8 shows the historical series of monthly-averaged ET o at the Mula-Embalse de la Cierva OWS. The series extends throughout the period 1933-2008 with the only noteworthy discontinuity being caused by the Spanish Civil War (1936)(1937)(1938)(1939).   (Saeed, 1986;Hargreaves & Allen, 2003). Therefore, calibration of the HG coefficient is widely required in the study area.

Discussion
The local calibration at each AWS by a global coefficient allowed estimating ET o values (ET o,HGg ) very close to ET o,PM (average RE = 9.61%). This statistic represents a higher improvement than the ob-served by Vanderlinden et al. (2004) in Andalusia, where a similar calibration process was followed. The spatial variation of C g in the SRB is in agreement with the results for other Mediterranean regions like Andalusia (Vanderlinden et al., 2004) and southern Italy (Mendicino & Senatore, 2013), where significantly higher C g values were found at coastal than at inland stations. This coastal effect was related with the habitual windier conditions at coastal locations (Vanderlinden et al., 2004;Gavilán et al., 2006), which usually produces systematic underestimations of HG equation (Martínez-Cob & Tejero-Juste, 2004;Shahidian et al., 2013).
The local calibration with monthly coeff icients (ET o,HGa ) resulted in even better estimations than those A b a n il la Á g u il a s A lc a n t a r il la A lm o r a d i B e n ie l C a r a v a c a C e h e g in H e ll ín M u r c ia S a n C a y e t a n o Y e c la T o t a n a A v e r a g e v a lu e [7] to the global regional function for the HG coefficient; and Eq.
[8] to the monthly regional functions for the HG coefficient. The results are at the considered 12 modern automatic weather stations (AWSs). obtained by the global coefficient (average RE = 7.71%), which meant an excellent behaviour, bearing in mind that only air temperature data were used. This satisfactory result was in accordance with previous attempts at locally calibrating the HG equation (e.g., ASCE, 1996;Jensen et al., 1997;Itenfisu et al., 2003;Jabloun & Sahli, 2008;Sentelhas et al., 2010;Cobaner, 2011;Espadafor et al., 2011). The annual pattern of C m,j agrees with that reported by Shahidian et al. (2013) in California, a location with similar climate to the SRB, and where a parabolic function was fitted to the C m,j annual trend.
In spite of those ET o accurate estimates, the local calibration were not enough for the purpose of this study, which entailed the estimation of ET o long-term series from historical temperature data in a set of old OWSs, which were placed in different locations that the AWSs. In that sense, two approaches for the regional estimation of the HG coefficient based on global and monthly regional functions were tested. Both methods were established in relation with the annual and the monthly means of daily temperature range, respectively, since they were highly correlated parameters with C g and C m,j . Contrary to other authors (Droogers & Allen, 2002) a good fitting with rainfall was not found, probably due to its high seasonal and inter-annual variability within the studied area. It should be noted that, in agreement with the results by Martínez-Cob & Tejero-Juste (2004) and Gavilán et al. (2006) in other semi-arid Spanish regions, a significant correlation was also found with wind speed, but it was not considered in the regionalization process because of the OWSs did not record such information.
Linear functions for global (Eq. [7]) and monthly (Eq. [8]) regional functions were used since, opposite to other authors (Samani, 2000;Vanderlinden et al., 2004;Mendicino & Senatore, 2013), quadratic or power functions did not provide significant correlation improvement. Other regional functions proposed in previous regionalization attempts systematically underestimated HG coefficient, even though some of them had been previously proposed for very nearby regions. This result could be ascribed to two main reasons: (i) the higher aridity of the study area with respect to the regions where the other regionalization attempts were developed, and (ii) the data treatment, since our study managed monthly-averaged weather data instead of daily data. However, the monthly data treatment did not affect the quality of the estimations after the regionalization of the HG coefficient, which was similar to that obtained in the aforementioned attempts.
The assessment of the regionalization process by cross validation clearly provided better ET o estimations for the monthly (Eq. [8]) than for the global (Eq. [7]) regional function (average REs were 10.1 and 12.5% respectively). Moreover, the quality of the ET o estimations provided by the monthly regional functions was only slightly lower than the obtained with the global local calibration, which can be considered a very satisfactory performance. Therefore, although both approaches for the regionalization of the HG coeff icient allowed estimating ET o values close to those obtained with the PM method (Table 7), the use of the monthly approach is highly recommended in the SRB as it clearly performed better than the global one.
In conclusion, this study demonstrates that the HG coefficient presents a great spatial variability throughout the SRB, which mainly depends on temperature related parameters. Eqs.
[7] and [8] can be used as an approximation to estimate global or monthly HG coefficient, respectively, at stations in the SRB where only minimum and maximum temperature data are available, which hence can be incorporated into the HG equation to estimate long-term series of monthly ET o . The proposed methodology appears the most straightforward and can be extended to other regions and climates, provided that the regional functions for calculating the HG coefficient are available. Therefore, the consideration of a regionally calibrated version of the HG equation is encouraged for the calculation of longterm series of ET o when a full dataset of air temperature, relative humidity, wind speed and solar radiation is not provided.