Estimation of evapotranspiration by the Food and Agricultural Organization of the United Nations (FAO) Penman–Monteith temperature (PMT) and Hargreaves–Samani (HS) models under temporal and spatial criteria – a case study in Duero basin (Spain)

The evapotranspiration-based scheduling method is the most common method for irrigation programming in agriculture. There is no doubt that the estimation of the reference evapotranspiration (ET_o) is a key factor in irrigated agriculture. However, the high cost and maintenance of agrometeorological stations and high number of sensors required to estimate it make it non-plausible, especially in rural areas. For this reason, the estimation of ET_o using air temperature, in places where wind speed, solar radiation and air humidity data are not readily available, is particularly attractive. A daily data record of 49 stations distributed over Duero basin (Spain), for the period 2000–2018, was used for estimation of ET_o based on seven models against Penman–Monteith (PM) FAO 56 (FAO – Food and Agricultural Organization of the United Nations) from a temporal (annual or seasonal) and spatial perspective. Two Hargreaves–Samani (HS) models, with and without calibration, and five Penman–Monteith temperature (PMT) models were used in this study. The results show that the models' performance changes considerably, depending on whether the scale is annual or seasonal. The performance of the seven models was acceptable from an annual perspective (R²>0.91, NSE > 0.88, MAE < 0.52 and RMSE < 0.69 mm d⁻¹; NSE – Nash–Sutcliffe model efficiency; MAE – mean absolute error; RMSE – root-mean-square error). For winter, no model showed good performance. In the rest of the seasons, the models with the best performance were the following three models: PMT_CUH (Penman–Monteith temperature with calibration of Hargreaves empirical coefficient – k_RS, average monthly value of wind speed, and average monthly value of maximum and minimum relative humidity), HS_C (Hargreaves–Samani with calibration of k_RS) and PMT_OUH (Penman–Monteith temperature without calibration of k_RS, average monthly value of wind speed and average monthly value of maximum and minimum relative humidity). The HS_C model presents a calibration of the Hargreaves empirical coefficient (k_RS). In the PMT_CUH model, k_RS was calibrated and average monthly values were used for wind speed and maximum and minimum relative humidity. Finally, the PMT_OUH model is like the PMT_CUH model except that k_RS was not calibrated. These results are very useful for adopting appropriate measures for efficient water management, especially in the intensive agriculture in semi-arid zones, under the limitation of agrometeorological data.

1 Introduction

A growing population and its need for food increasingly demand natural resources such as water. This, linked with the uncertainty of climate change, makes water management a key consideration for future food security. The main challenge is to produce enough food for a growing population that is directly affected by the challenges created by the management of agricultural water, mainly by irrigation management (Pereira, 2017).

Evapotranspiration (ET) is the water lost from the soil surface and surface leaves by evaporation and, by transpiration, from vegetation. ET is one of the major components of the hydrologic cycle and represented a loss of water from the drainage basin. ET information is key to understanding and managing water resource systems (Allen et al., 2011). ET is normally modeled using weather data and algorithms that describe aerodynamic characteristics of the vegetation and surface energy.

In agriculture, irrigation water is usually applied based on the water balance method in the soil water balance equation, which allows the calculation of the decrease in soil water content as the difference between outputs and inputs of water to the field. In arid areas where rainfall is negligible during the irrigation season, an average irrigation calendar may be defined a priori using mean ET values (Villalobos et al., 2016). The Food and Agricultural Organization of the United Nations (FAO) improved and upgraded the methodologies for reference evapotranspiration (ET_o) estimation by introducing the reference crop (grass) concept, described by the FAO Penman–Monteith (PM-ET_o) equation (Allen et al., 1998). This approach was tested well under different climates and time step calculations and is currently adopted worldwide (Allen et al., 1998; Todorovic et al., 2013; Almorox et al., 2015). Estimated crop evapotranspiration (ET_c) is obtained by a function of two factors (ET${}_{\mathrm{c}}=K_{\mathrm{c}}\cdot {\mathrm{ET}}_{\mathrm{o}}$): reference crop evapotranspiration (ET_o) and the crop coefficient (K_c; Allen et al., 1998). ET_o was introduced to study the evaporative demand of the atmosphere independently of crop type, crop stage development and management practices. ET_o is only affected by climatic parameters and is computed from weather data. Crop influences are accounted for by using a specific crop coefficient (K_c). However, K_c varies predominately with the specific crop characteristics and only to a limited extent with climate (Allen et al., 1998).

The ET is very variable locally and temporally because of the climate differences. Because the ET component is relatively large in water hydrology balances, any small error in its estimate or measurement represents large volumes of water (Allen et al., 2011). Small deviations in ET_o estimations affect irrigation and water management in rural areas in which crop extension is significant. For example, in 2017 there was a water shortage at the beginning of the cultivation period (March) at the Duero basin (Spain). The classical irrigated crops, i.e., corn, were replaced by others with lower water needs, such as sunflower.

Wind speed (u), solar radiation (R_s), relative humidity (RH) and temperature (T) of the air are required to estimate ET_o. Additionally, the vapor pressure deficit (VPD), soil heat flux (G) and net radiation (R_n) measurements or estimates are necessary. The PM-ET_o methodology presents the disadvantage that required climate or weather data are normally unavailable or of low quality (Martinez and Thepadia, 2010) in rural areas. In this case, where data are missing, Allen et al. (1998) in the guidelines for PM-ET_o recommend two approaches: (a) using the equation of Hargreaves–Samani (Hargreaves and Samani, 1985) and (b) using the Penman–Monteith temperature (PMT) method that requires data of temperature to estimate R_n (net radiation) and VPD for obtaining ET_o. In these situations, temperature-based evapotranspiration (TET) methods are very useful (Mendicino and Senatore, 2013). Air temperature is the most available meteorological data, which are available from most climatic weather stations. Therefore, TET methods and temperature databases are a solid base for ET estimation all over the world, including areas with limited data resources (Droogers and Allen, 2002). The first reference of the use of PMT for limited meteorological data was Allen (1995); subsequently, studies like those of Annandale et al. (2002) were carried out, having similar behavior to the Hargreaves–Samani (HS) method and FAO PM, although there was the disadvantage of requiring greater preparation and computation of the data than the HS method. Regarding this point, it should be noted that the researchers do not favor using the PMT formulation and adopt the HS equation, which is simpler and easier to use (Paredes et al., 2018). Authors like Pandey et al. (2014) performed calibrations based on solar radiation coefficients in Hargreaves–Samani equations. Today, the PMT calculation process is easily implemented with the new computers (Pandey and Pandey, 2016; Quej et al., 2019).

Todorovic et al. (2013) reported that, in Mediterranean hyper-arid and arid climates, PMT and HS show similar behavior and performance, while for moist sub-humid areas, the best performance was obtained with the PMT method. This behavior was reported for moist sub-humid areas in Serbia (Trajovic, 2005). Several studies confirm this performance in a range of climates (Martinez and Thepadia, 2010; Raziei and Pereira, 2013; Almorox et al., 2015; Ren et al., 2016). Both models (HS and PMT) improved when local calibrations were performed (Gavilán et al., 2006; Paredes et al., 2018). These reduce the problem when wind speed and solar radiation are the major driving variables.

Studies in Spain comparing HS and PMT methodologies were studied in moist sub-humid climate zones (northern Spain), showing a better fit in PMT than in HS. (López-Moreno et al., 2009). Tomas-Burguera et al. (2017) reported, for the Iberian Peninsula, a better adjustment of PMT than HS, provided that the lost values were filled by interpolation and not by estimation in the model of PMT.

Normally the calibration of models for ET_o estimation is done from a spatial approach, calibrating models in the locations studied. Very few studies have been carried out to test models from the seasonal point of view, with the annual calibration being the most studied. Meanwhile spatial and annual approaches are of great interest to climatology and meteorology, as agriculture, seasonal or even monthly calibrations are relevant to crops (Nouri and Homaee, 2018). To improve the accuracy of ET_o estimations, Paredes et al. (2018) used the values of the calibration constant values in the models that were derived for the October–March and April–September periods.

The aim of this study was to evaluate the performance of temperature models for the estimation of reference evapotranspiration against the FAO 56 Penman–Monteith model, with a temporal (annual or seasonal) and spatial perspective in the Duero basin (Spain). The models evaluated were two HS models, with calibration and without calibration, and five PMT models analyzing the contribution of wind speed, humidity and solar radiation in a situation of limited agrometeorological data.

Figure 1 Location of study area. The point with the number indicates the location of the agrometeorological stations according to Table 1.

2 Materials and method

2.1 Description of the study area

The study focuses on the Spanish part of the Duero hydrographic basin. The international hydrographic Duero basin is the most extensive of the Iberian Peninsula; with an area 98 073 km², it includes the territory of the Duero River basin as well as the transitional waters of the Porto estuary and the associated Atlantic coastal waters (CHD, 2019). It is a shared territory between Portugal, with 19 214 km² (19.6 % of the total area), and Spain, with 78 859 km² (80.4 % of the total area). The Duero River basin is located in Spain, between 43^∘5^′ and 40^∘10^′ N and 7^∘4^′ and 1^∘50^′ W (Fig. 1). This basin aligns almost exactly with the so-called Submeseta Norte, an area with an average altitude of 700 m, delimited by mountain ranges with a much drier central zone that contains large aquifers, which is the most important area of agricultural production; 98.4 % of the Duero basin belongs to the autonomous community of Castilla and Léon, and 70 % of the average annual precipitation is used directly by the vegetation or evaporated from surface – this represents 35 000 h m³. The remaining (30 %) is the total natural runoff. The Mediterranean climate is the predominant climate; 90 % of the surface is affected by summer drought conditions. The average annual values are a temperature of 12 ^∘C and precipitation of 612 mm. However, precipitation ranges from minimum values of 400 mm (south–central area of the basin) to a maximum of 1800 mm in the northeast of the basin (CHD, 2019). According to Lautensach (1967), 30 mm is the threshold definition of a dry month. Therefore, between two and five dry periods can be found in the basin (Ceballos et al., 2004). Moreover, the climate variability, especially precipitation, exhibited in the last decade has decreased the water availability for irrigation in this basin (Segovia-Cardozo et al., 2019).

The Duero basin has 4×10⁶ ha of rainfed crops and some 500 000 ha of irrigated crops that consume 75 % of the basin's water resources. Barley (Hordeum vulgare L.) is the most important rainfed crop in the basin, occupying 36 % of the national crop surface, followed by wheat (Triticum aestivum L.), with 30 % (MAPAMA, 2019). Sunflower (Helianthus annuus L.) represents 30 % of the national crop surface. This crop is mainly unirrigated (90 %). Maize (Zea mays L.), alfalfa (Medicago sativa L.) and sugar beet (Beta vulgaris L. var. sacharifera) are the main irrigated crops. These crops represents 29 %, 30 % and 68 % of the national crop area, respectively. Finally, vines (Vitis vinifera L.) fill 72 000 ha, being less than 10 % irrigated. For the irrigated crops of the basin there are water allocations that fluctuate depending on the availability of water during the agricultural year and the type of crop. These values fluctuate from 1200–1400 m³ ha⁻¹ for vines to 6400–7000 m³ ha⁻¹ for maize and alfalfa. The use rates of the irrigation systems used in the basin are as follows: 25 %, 68 % and 7 % for surface, sprinkler and drip irrigation, respectively (Plan Hidrológico, 2019).

2.2 Meteorological data

The daily climate data were downloaded from 49 stations (Fig. 1b) from the agrometeorological network SIAR (Agroclimatic Information System for Irrigation; SIAR in Spanish), which is managed by the Spanish Ministry of Agriculture, Fisheries and Food (SIAR, 2018), providing the basic meteorological data from weather stations distributed throughout the Duero basin (Table 1). Each station incorporates measurements of air temperature (T) and relative humidity (RH; Vaisala HMP155), precipitation (ARG100 rain gauge), solar global radiation (pyranometer Skye SP1110), and wind direction and wind speed (u; wind vane and R.M. Young 05103 anemometer). Sensors were periodically maintained and calibrated, and all data were recorded and averaged hourly on a data logger (Campbell CR10X and CR1000). Characteristics of the agrometeorological stations were described by (Moratiel et al., 2011, 2013a). For quality control, all parameters were checked, and the sensors were periodically maintained and calibrated, with all data being recorded and averaged hourly on a data logger. The database calibration and maintenance are carried out by the Ministry of Agriculture. Transfer of data from stations to the main center is accomplished by modems; the main center incorporates a server which sequentially connects to each station to download the information collected during the day. Once the data from the stations are downloaded, they are processed and transferred to a database. The main center is responsible for quality control procedures that comprise the routine maintenance program of the network, including sensor calibration, checking for validity values and data validation. Moreover, the database was analyzed to find incorrect or missing values. To ensure that high-quality data were used, we used quality control procedures to identify erroneous and suspect data. The quality control procedures applied are the range and limit test, step test, and internal consistency test (Estevez et al., 2016).

The period studied was from 2000 to 2018, although the start date may fluctuate depending on the availability of data. Table 1 shows the coordinates of the agrometeorological stations used in the Duero basin and the aridity index based on UNEP (1997). Table 1 shows the predominance of the semi-arid climate zone, with 42 stations of the 49 being semi-arid, 2 being arid, 4 being dry sub-humid and 1 being moist sub-humid.

2.3 Estimates of reference evapotranspiration

2.3.1 FAO Penman–Monteith (FAO PM)

The FAO recommends the PM method for computing ET_o and evaluating other ET_o models like the Penman–Monteith model using only temperature data (PMT) and other temperature-based models (Allen et al., 1998). The method estimates the potential evapotranspiration from a hypothetical crop with an assumed height of 0.12 m, having an aerodynamic resistance of (r_a) 208∕u₂ (u₂ is the mean daily wind speed measured at a 2 m height over the grass), a surface resistance (r_s) of 70 s m⁻¹ and an albedo of 0.23, closely resembling the evaporation of an extension surface of green grass with a uniform height that is actively growing and adequately watered. The ET_o (mm d⁻¹) was estimated following FAO 56 (Allen et al., 1998):

$$\begin{array}{}\text{(1)}& {\mathrm{ET}}_{\mathrm{o}}={\displaystyle \frac{\mathrm{0.408}\mathrm{\Delta }\left(R_{\mathrm{n}}-G\right)+\mathit{\gamma }\frac{\mathrm{900}}{T+\mathrm{273}}{u}_{\mathrm{2}}\left(e_{\mathrm{s}}-e_{\mathrm{a}}\right)}{\mathrm{\Delta }+\mathit{\gamma }\left(\mathrm{1}+\mathrm{0.34}{u}_{\mathrm{2}}\right)}}.\end{array}$$

In Eq. (1), R_n is net radiation at the surface (MJ m⁻² d⁻¹), G is ground heat flux density (MJ m⁻² d⁻¹), γ is the psychrometric constant (kPa ^∘C⁻¹), T is mean daily air temperature at 2 m height (^∘C), u₂ is wind speed at 2 m height (m s⁻¹), e_s is the saturation vapor pressure (kPa), e_a is the actual vapor pressure (kPa) and Δ is the slope of the saturation vapor pressure curve (kPa ^∘C⁻¹). According to Allen et al. (1998) in Eq. (1), G can be considered to be zero.

Table 1 Agrometeorological station used in the study. Coordinates and aridity index are shown.

2.3.2 Hargreaves–Samani (HS)

The scarcity of available agrometeorological data (mainly global solar radiation, air humidity and wind speed) limit the use of the FAO PM method in many locations. Allen et al. (1998) recommended applying the Hargreaves–Samani expression in situations where only the air temperature is available. The HS formulation is an empirical method that requires empirical coefficients for calibration (Hargreaves and Samani, 1982, 1985). The Hargreaves–Samani (Hargreaves and Samani, 1982, 1985) method is given by the following equation:

$$\begin{array}{}\text{(2)}& {\mathrm{ET}}_{\mathrm{o}}=\mathrm{0.0135}\cdot k_{\mathrm{RS}}\cdot \mathrm{0.408}\cdot H_{\mathrm{o}}\cdot \left(T_m+\mathrm{17.8}\right)\cdot {\left(T_x-T_n\right)}^{\mathrm{0.5}},\end{array}$$

where ET_o is the reference evapotranspiration (mm d⁻¹); H_o is extraterrestrial radiation (MJ m⁻² d⁻¹); k_RS is the Hargreaves empirical coefficient; and T_m, T_x and T_n are the daily mean, maximum air temperature and minimum air temperature (^∘C), respectively. The value k_RS was initially set to 0.17 for arid and semi-arid regions (Hargreaves and Samani, 1985). Hargreaves (1994) later recommended using the value of 0.16 for interior regions and 0.19 for coastal regions. Daily temperature variations can occur due to other factors such as topography, vegetation, humidity, etc.; thus using a fixed coefficient may lead to errors. In this study, we use 0.17 as the original coefficient (HS_O) and the calibrated coefficient k_RS (HS_C). k_RS reduces the inaccuracy, thus improving the estimation of ET_o. This calibration was done for each station.

2.3.3 Penman–Monteith temperature (PMT)

The FAO PM, when applied using only measured temperature data, is referred to as PMT and retains many of the dynamics of the full-data FAO PM (Pereira et al., 2015; Hargreaves and Allen, 2003). Humidity and solar radiation are estimated in the PMT model using only air temperature as input for the calculation of ET_o. Wind speed in the PMT model is set to the constant value of 2 m s⁻¹ (Allen et al., 1998). In this model, where global solar radiation (or sunshine data) is lacking, the difference between the maximum and minimum temperature can be used, as an indicator of cloudiness and atmospheric transmittance, for the estimation of solar radiation (Eq. 3; Hargreaves and Samani, 1982). Net solar shortwave and longwave radiation estimates are obtained as indicated by Allen et al. (1998), in Eqs. (4) and (5), respectively. The expression of PMT is obtained as indicated in the following:

$$\begin{array}{}\text{(3)}& {\displaystyle }& {\displaystyle }{R}_{\mathrm{s}}=H_{\mathrm{o}}\cdot k_{\mathrm{RS}}\cdot {\left(T_x-T_n\right)}^{\mathrm{0.5}},\text{(4)}& {\displaystyle }& {\displaystyle }{R}_{\mathrm{ns}}=\mathrm{0.77}\cdot H_{\mathrm{o}}\cdot k_{\mathrm{RS}}\cdot {\left(T_x-T_n\right)}^{\mathrm{0.5}},\end{array}$$

where R_s is solar radiation (MJ m⁻² d⁻¹); R_ns is net solar shortwave radiation (MJ m⁻² d⁻¹); H_o is extraterrestrial radiation (MJ m⁻² d⁻¹); and H_o was computed as a function of site latitude, solar angle and the day of the year (Allen et al., 1998). T_x is the daily maximum air temperature (^∘C), and T_n is the daily minimum air temperature (^∘C). For k_RS Hargreaves (1994) recommended using k_RS=0.16 for interior regions and k_RS=0.19 for coastal regions. For better accuracy the coefficient k_RS can be adjusted locally (Hargreaves and Allen, 2003). In this study two assumptions of k_RS were made: one where a value of 0.17 was fixed and another where it was calibrated for each station,

$$\begin{array}{}\text{(5)}& \begin{array}{rl}{R}_{\mathrm{nl}}& =\left(\mathrm{1.35}\cdot \left({\displaystyle \frac{k_{\mathrm{RS}}\cdot {\left(T_x-T_n\right)}^{\mathrm{0.5}}}{\mathrm{0.75}-\mathrm{2}z{\mathrm{10}}^{-\mathrm{5}}}}\right)-\mathrm{0.35}\right)\\ & \cdot \left(\mathrm{0.34}-\mathrm{0.14}{\left(\mathrm{0.6108}\cdot \exp \left({\displaystyle \frac{\mathrm{17.27}\cdot T_{\mathrm{d}}}{T_{\mathrm{d}}-\mathrm{237.3}}}\right)\right)}^{\mathrm{0.5}}\right)\\ & \cdot \mathit{\sigma }\cdot \left({\displaystyle \frac{{\left(T_x+\mathrm{273.15}\right)}^{\mathrm{4}}+{\left(T_n+\mathrm{273.15}\right)}^{\mathrm{4}}}{\mathrm{2}}}\right),\end{array}\end{array}$$

where R_nl is net longwave radiation (MJ m⁻²  d⁻¹), T_x is daily maximum air temperature (^∘C), T_n is daily minimum air temperature (^∘C), T_d is dew point temperature (^∘C) calculated with the T_n according to Todorovic et al. (2013), σ is the Stefan–Boltzmann constant for a day ($\mathrm{4.903}\times {\mathrm{10}}^{-\mathrm{9}}$ MJ K⁻⁴ m⁻² d⁻¹) and z is the altitude (m):

$$\begin{array}{}\text{(6)}& {\mathrm{PMT}}_{\mathrm{rad}}=\left({\displaystyle \frac{\mathrm{0.408}\mathrm{\Delta }}{\mathrm{\Delta }+\mathit{\gamma }\left(\mathrm{1}+\mathrm{0.34}{u}_{\mathrm{2}}\right)}}\right)\cdot \left(R_{\mathrm{ns}}-R_{\mathrm{nl}}-G\right),\end{array}$$

$$\begin{array}{}\text{(7)}& {\displaystyle }& {\displaystyle }{\mathrm{PMT}}_{\mathrm{aero}}={\displaystyle \frac{\mathit{\gamma }\cdot \frac{\mathrm{900}\cdot u_{\mathrm{2}}}{T_{\mathrm{m}}+\mathrm{273}}\cdot \left(\left(\frac{e_{\mathrm{s}}\left(T_x\right)+e_{\mathrm{s}}\left(T_n\right)}{\mathrm{2}}\right)-e_{\mathrm{s}}\left(T_{\mathrm{d}}\right)\right)}{\mathrm{\Delta }+\mathit{\gamma }\left(\mathrm{1}+\mathrm{0.34}{u}_{\mathrm{2}}\right)}},\text{(8)}& {\displaystyle }& {\displaystyle }\mathrm{PMT}={\mathrm{PMT}}_{\mathrm{rad}}+{\mathrm{PMT}}_{\mathrm{aero}},\end{array}$$

where PMT is the reference evapotranspiration estimate by the Penman–Monteith temperature method (mm d⁻¹), PMT_rad is the radiative component of PMT (mm d¹), PMT_aero is the aerodynamic component of PMT (mm d⁻¹), Δ is the slope of the saturation vapor pressure curve (kPa ^∘C⁻¹), γ is the psychrometric constant (kPa ^∘C⁻¹), R_ns is net solar shortwave radiation (MJ m⁻² d⁻¹), R_nl is net longwave radiation (MJ m⁻² d⁻¹), G is ground heat flux density (MJ m⁻² d⁻¹), considered to be zero according to Allen et al. (1998), T_m is mean daily air temperature (^∘C), T_x is maximum daily air temperature, T_n is mean daily air temperature, T_d is dew point temperature (^∘C) calculated with the T_n according to Todorovic et al. (2013), u₂ is wind speed at 2 m height (m s⁻¹) and e_s is the saturation vapor pressure (kPa). In this model two assumptions of k_RS were made: one where a value of 0.17 was fixed and another where it was calibrated for each station.

2.3.4 Calibration and models

We studied two methods to estimate the ET_o: the HS method and the reference evapotranspiration estimate by PMT. Within these methods, different adjustments are proposed based on the adjustment coefficients of the methods and the missing data. The parametric calibration for the 49 stations was applied in this study. In order to decrease the errors of the evapotranspiration estimates, local calibration was used. The seven methods used with the coefficient (k_RS) of the calibrated and characteristics in the different locations studied are shown in Table 2. The calibration of the model coefficients was achieved by the nonlinear least-squares fitting technique. The analysis was made on a yearly and seasonal basis. The seasons were the following: (1) winter (December, January and February, or DJF), (2) spring (March, April and May, or MAM), (3) summer (June, July and August, or JJA) and (4) autumn (September, October and November, or SON).

Table 2 Characteristics of the models used in this study.

¹ Dew point temperature obtained according to Todorovic et al. (2013). ² Average monthly value of wind speed. ³ Average monthly value of maximum and minimum relative humidity.

2.4 Performance assessment

The model's suitability, accuracy and performance were evaluated using the coefficient of determination (R²; Eq. 9) of the n pairs of observed (O_i) and predicted (P_i) values. Also, the mean absolute error (MAE; mm d⁻¹; Eq. 10), root-mean-square error (RMSE; Eq. 11) and the Nash–Sutcliffe model efficiency (NSE; Eq. 12; Nash and Sutcliffe, 1970) coefficient were used. The coefficient-of-regression line (b), forced through the origin, is obtained by predicted values divided by observed values (ET_model∕ET_FAO56).

The results were represented in a map applying the Kriging method with the Surfer^® 8 program:

$$\begin{array}{}\text{(9)}& {\displaystyle }& {\displaystyle }{R}^{\mathrm{2}}={\left\{{\displaystyle \frac{\sum _{i=\mathrm{1}}^n\left(O_i-\stackrel{\mathrm{‾}}{O}\right)\cdot \left(P_i-\stackrel{\mathrm{‾}}{P}\right)}{{\left[\sum _{i=\mathrm{1}}^n{\left(O_i-\stackrel{\mathrm{‾}}{O}\right)}^{\mathrm{2}}\right]}^{\mathrm{0.5}}\cdot {\left[\sum _{i=\mathrm{1}}^n{\left(P_i-\stackrel{\mathrm{‾}}{P}\right)}^{\mathrm{2}}\right]}^{\mathrm{0.5}}}}\right\}}^{\mathrm{2}},\text{(10)}& {\displaystyle }& {\displaystyle }\mathrm{MAE}={\displaystyle \frac{\mathrm{1}}{n}}\sum _{i=\mathrm{1}}^n\left(\left|O_i-P_i\right|\right)\left(\mathrm{mm}{\mathrm{d}}^{-\mathrm{1}}\right),\text{(11)}& {\displaystyle }& {\displaystyle }\mathrm{RMSE}={\left[{\displaystyle \frac{\sum _{i=\mathrm{1}}^n{\left(O_i-P_i\right)}^{\mathrm{2}}}{n}}\right]}^{\mathrm{0.5}}\left(\mathrm{mm}{\mathrm{d}}^{-\mathrm{1}}\right),\text{(12)}& {\displaystyle }& {\displaystyle }\mathrm{NSE}=\mathrm{1}-\left[{\displaystyle \frac{\sum _{i=\mathrm{1}}^n{\left(O_i-P_i\right)}^{\mathrm{2}}}{\sum _{i=\mathrm{1}}^n{\left(O_i-\stackrel{\mathrm{‾}}{O}\right)}^{\mathrm{2}}}}\right].\end{array}$$

3 Results and discussion

In the study period the data indicated that the Duero basin is characterized by being a semi-arid climate zone (94 % of the stations), where the P∕ET_o ratio is between 0.2 and 0.5 (Todorovic et al., 2013). The mean annual rainfall is 428 mm, while the average annual ET_o for Duero basin is 1079 mm, reaching the maximum values in the center–south zone, with values that slightly surpass 1200 mm (Fig. 2). The great temporal heterogeneity is observed in the Duero basin, with values of 7 % of the ET_o during the winter months (DJF), while during the summer months (JJA) they represent 47 % of the annual ET_o. In addition, the months from May to September represent 68 % of the annual ET_o, with similar values to those reported by Moratiel et al. (2011).

Figure 2 Mean values season of ET_o (mm) during the study period 2000–2018. (a) Annual, (b) winter (December, January and February, or DJF), (c) spring (March, April and May, or MAM), (d) summer (June, July and August, or JJA) and (e) autumn (September, October and November, or SON).

Table 3 shows the different statistics analyzed in the seven models studied as a function of the season of the year and annually. From an annual point of view all models show R² values higher than 0.91, an NSE higher than 0.88, a MAE less than 0.52, a RMSE lower than 0.69, and underestimates or overestimates of the models by ±4 %. The best behavior is shown by the PMT_CHU model, with a MAE and RMSE of 0.39 and 0.52 mm d⁻¹, respectively. PMT_CUH shows no tendency to overestimate or underestimate the values in which a coefficient of regression b of 1.0 is observed. This model shows values of the NSE and R² of 0.93. The models HS_C and PMT_OUH have similar behavior, with the same MAE (0.41 mm d⁻¹), NSE (0.92) and R² (0.91). The RMSE is 0.55 mm d⁻¹ for the PMT_OUH model and 0.54 mm d⁻¹ for the HS_C model. The models PMT_OUT and HS_O showed slightly better performance than PMT_O2T and PMT_C2T, given that the last two models showed the worst behavior (Fig. 3). The performance of the models (PMT_O2T, PMT_OUT and PMT_OUH) improves as the averages of wind speed (u) and dew temperature (T_d) values are incorporated. The same pattern is shown between the PMT_CUH models, where the mean u values and T_d are incorporated, and PMT_C2T, where u is 2 m s⁻¹ and dew temperature is calculated with the approximation of Todorovic et al. (2013). These adjustments are supported because the adiabatic component of evapotranspiration in the PMT equation is very influential in the Mediterranean climate, especially wind speed (Moratiel et al., 2010).

Figure 3 Performance of the models with an annual focus. (a) Average annual values of ET_o (mm d⁻¹). Mean values of MAE (mm d⁻¹): (b) PMT_O2T model, (c) H_O model, (d) H_C model, (e) PMT_C2T model, (f) PMT_OUT model, (g) PMT_OUH model and (h) PMT_CUH model.

Figure 4 Performance of the models with a winter focus (December, January and February). (a) Average values of ET_o (mm d⁻¹) in winter. Mean values of MAE (mm d⁻¹): (b) PMT_O2T model, (c) H_O model, (d) H_C model, (e) PMT_C2T model, (f) PMT_OUT model, (g) PMT_OUH model and (h) PMT_CUH model.

Table 3 Statistical indicators for ET_o estimation in the seven models studied for different seasons. Average data for the 49 stations studied.

From a spatial perspective, it is observed in Fig. 3 that the areas where the values of the MAE are higher are to the east and southwest of the basin. This is due to the fact that the average wind speed in the eastern zone is higher than 2.5 m s⁻¹; for example, the Hinojosa del Campo station shows average annual values of 3.5 m s⁻¹. The southwestern area shows values of wind speeds below 1.5 m s⁻¹, such as at the Ciudad Rodrigo station, with annual average values of 1.19 m s⁻¹.

These MAE differences are more pronounced in the models in which the average wind speed is not taken, such as the PMT_C2T and PMT_O2T models. Most of the basin has values of wind speeds between 1.5 and 2.5 m s⁻¹. The lower MAE values in the northern zone of the basin are due to the lower average values of the VPD than the central area, with values of 0.7 kPa in the northern zone and 0.95 kPa in the central zone. The same trends in the effect of wind on the ET_o estimates were detected by Nouri and Homaee (2018), who indicated that values outside the range of 1.5–2.5 m s⁻¹ in models where the default u was set at 2 m s⁻¹ increased the error of the ET_o. Even for models such as HS, where the influence of the wind speed values is not directly indicated outside of the ranges previously mentioned, their performance is not good, and some authors have proposed HS calibrations based on wind speeds in Spanish basins such as the Ebro Basin (Martínez-Cob and Tejero-Juste, 2004). In our study, the HS_C model showed good performance, with MAE values similar to PMT_CUH and PMT_OUH (Fig. 3). The performance of the models by season of the year changes considerably, obtaining lower adjustments, with values of R²=0.53 for winter (DJF) in the models HS_O and HS_C and for summer (JJA) in the models PMT_O2T and PMT_C2T. All models during spring and autumn show R² values above 0.8. The NSE for models HS_O, PMT_C2T, PMT_O2T and PMT_OUT in summer and winter is at unsatisfactory values below 0.5 (Moriasi et al., 2007). The mean values (49 stations) of the MAE (Fig. 4) and RMSE for the models in the winter were 0.24–0.30 and 0.3–0.37 mm d⁻¹, respectively. For spring, the ranges were between 0.42 and 0.52 mm d⁻¹ for the MAE (Fig. 5) and 0.55 and 0.65 mm d⁻¹ for RMSE. In summer, the MAE (Fig. 6) fluctuated between 0.53 and 0.72 mm d⁻¹, and the RMSE fluctuated between 0.68 and 0.91 mm d⁻¹. Finally, in autumn, the values of the MAE (Fig. 7) and RMSE were 0.38–0.58 and 0.49–0.70 mm d⁻¹, respectively (Table 3).

Figure 5 Performance of the models with a spring focus (March, April and May). (a) Average annual values of ET_o (mm d⁻¹) in spring. Mean values of MAE (mm d⁻¹): (b) PMT_O2T model, (c) H_O model, (d) H_C model, (e) PMT_C2T model, (f) PMT_OUT model, (g) PMT_OUH model and (h) PMT_CUH model.

Figure 6 Performance of the models with a summer focus (June, July and August). (a) Average values of ET_o (mm d⁻¹) in summer. Mean values of MAE (mm d⁻¹): (b) PMT_O2T model, (c) H_O model, (d) H_C model, (e) PMT_C2T model, (f) PMT_OUT model, (g) PMT_OUH model and (h) PMT_CUH model.

Figure 7 Performance of the models with an autumn focus (September, October and November). (a) Average values of ET_o (mm d⁻¹) in autumn. Mean values of MAE (mm d⁻¹): (b) PMT_O2T model, (c) H_O model, (d) H_C model, (e) PMT_C2T model, (f) PMT_OUT model, (g) PMT_OUH model and (h) PMT_CUH model.

The model that shows the best performance independently of the season is PMT_CUH. The models that can be considered in a second level are HS_C and PMT_OUH. During the months of more solar radiation (summer and spring) the performance of the HS_C model is slightly better than the PMT_OUH model. The HS_O, PMT_O2T, PMT_C2T and PMT_OUT models have a much poorer performance than the previous models (PMT_OUH and HS_C). The model that has the worst performance is PMT_O2T.

The northern area of the basin is the area in which a lower MAE shows in most models and for all seasons. This is due in part to the fact that the lower values of ET_o (mm d⁻¹) are located in the northern zone. On the other hand, the eastern zone of the basin shows the highest values of the MAE due to the strong winds that are located in that area.

During the winter the seven models tested show no great differences between them, although PMT_CUH is the model with the best performance. It is important to indicate that during this season the RMSE (%) is placed in all the models above 30 %, so they can be considered to be very weak models. According to Jamieson et al. (1991) and Bannayan and Hoogenboom (2009) the model is considered excellent with a normalized RMSE (%) less than 10 %, good if the normalized RMSE (%) is greater than 10 and less than 20 %, fair if the normalized RMSE (%) is greater than 20 % and less than 30 %, and poor if the normalized RMSE (%) is greater than 30 %. All models that are made during the spring season (MAM) can be considered to be good or fair, since their RMSE (%) fluctuates between 17 % and 20 %. The seven models that are made during summer season (JJA) can be considered to be good, since their RMSE varies from 12 % to 16 %. Finally, the models that are made during autumn (SON) are considered to be fair or poor, fluctuating between 22 % and 32 %. The models that reached values greater than 30 % during autumn were PMT_C2T (31 %) and PMT_O2T (32 %), which also had a clear tendency to overestimate (Table 3). In the use of temperature models for estimating ET_o, it is necessary to know the objective that is set. For the management of irrigation in crops, it is better to test the models in the period in which the species require the contribution of additional water. In many cases, applying the models with an annual perspective with good performance can lead to more accentuated errors in the period of greater water needs. The studies of different temporal and spatial scales of the temperature models for ET_o estimation can give valuable information that allows for managing the water planning in zones where the economic development does not allow the implementation of agrometeorological stations due to its high cost.

4 Discussion

In annual seasons our data of the RMSE fluctuate from 0.69 mm d⁻¹ (PMT_O2T) to 0.52 mm d⁻¹ (PMT_O2T). These data are in accordance with the values cited by other authors in the same climatic zone. Jabloun and Sahli (2008) cited a RMSE of 0.41–0.80 mm d⁻¹ for Tunisia. The authors showed the PMT model performance to be better than that for the Hargreaves non-calibrated model. Raziei and Pereira (2013) reported data of the RMSE for a semi-arid zone in Iran to be between 0.27 and 0.81 mm d⁻¹ for the HS_C model and 0.30 and 0.79 mm d⁻¹ for PMT_C2T, although these authors use monthly averages in their models. Ren et al. (2016) reported values of RMSE to be in the range of 0.51 to 0.90 mm d⁻¹ for PMT_C2T and in the range of 0.81 to 0.94 mm d⁻¹ for HS_C in semi-arid locations in Inner Mongolia (China). Todorovic et al. (2013) found the PMT_O2T method to have better performance than the uncalibrated HS method (HS_O), with a RMSE average of 0.47 mm d⁻¹ for PMT_O2T and 0.52 HS_O. At this point, we should highlight that in our study daily-value data were used.

The original Hargreaves equation was developed by regressing cool-season grass ET in Davis, California; the k_RS coefficient is a calibration coefficient. The aridity index for Davis is semi-arid ($P/\mathrm{ET}=\mathrm{0.33}$; Hargreaves and Allen, 2003; Moratiel et al., 2013b), like 94 % of the stations studied, which explains why the behavior of the HS_O model is often very similar to HS_C. Even so, the calibration coefficient needs to be adjusted for other climates. Numerous studies in the literature have demonstrated the relevance of the k_RS calibration model to estimating FAO 56 (Todorovic et al., 2013; Raziei and Pereira, 2013; Paredes et al., 2018).

PMT models have improved, considering the average wind speed. In addition, trends and fluctuations of u have been reported as the factor that most influences ET_o trends (Nouri et al., 2017; McVicar et al., 2012; Moratiel et al., 2011). Numerous authors have recommended including, as much as possible, average data of local wind speeds for the improvement of the models, like Nouri and Homaee (2018) and Raziei and Pereira (2013) in Iran, Paredes et al. (2018) in the Azores (Portugal), Djaman et al. (2017) in Uganda, Rojas and Sheffield (2013) in Louisiana (USA), Jabloun and Shali (2008) in Tunisia, and Martinez-Cob and Tejero-Juste (2004) in Spain, among others. In addition, even ET_o prediction models based on PMT focus their behavior based on the wind speed variable (Yang et al., 2019). It is important to note that PMT_OUT generally has better performance than PMT_C2T except for in spring. The difference between both models is that in PMT_C2T, k_RS is calibrated with wind speed set to 2 m s⁻¹, and in PMT_OUT, k_RS is not calibrated and has an average wind speed. In this case the wind speed variable has less of an effect than the calibration of k_RS, since the average values of wind during spring (2.3 m s⁻¹) are very close to 2 m s⁻¹ and there is no great variation between both settings. In this way, k_RS calibration shows a greater contribution than the average of the wind speed to improve the model (Fig. 5e and f). In addition, although u is not directly considered for HS, this model is more robust in regions with speed averages around 2 m s⁻¹ (Allen et al., 1998; Nouri and Homaee, 2018). On the other hand errors in the estimation of relative humidity cause substantial changes in the estimation of ET_o, as reported by Nouri and Homaee (2018) and Landeras et al. (2008).

The results of RMSE values (%) of the different models change considerably by season; values are between 16.6 % and 12.3 % for summer and between 41.2 % and 33.5 % for winter. Similar results were obtained in Iran by Nouri and Homaee (2018), where in the months of December–January and February the performance of the PMT and HS models tested had RMSE (%) values above 30 %. Very few studies, as far as we know, have been carried out with adjustments of evapotranspiration models from a temporal point of view, and generally the models are calibrated and adjusted from an annual point of view. Some authors, such as Aguilar and Polo (2011), differentiate seasons between wet and dry, and others, such as Paredes et al. (2018), divide them into summer and winter; Vangelis et al. (2013) take two periods into account, and Nouri and Homaee (2018) do it from a monthly point of view. In most cases, the results obtained in these studies are not comparable with those presented in this study, since the timescales are different. However, it can be noted that the results of the models according to the timescale season differ greatly with respect to the annual scale.

5 Conclusions

The performance of seven temperature-based models (PMT and HS) was evaluated in the Duero basin (Spain) for a total of 49 agrometeorological stations. Our studies revealed that the models tested on an annual or seasonal basis provide different performance. The values of R² are higher when they are performed annually, with values between 0.91 and 0.93 for the seven models, but when performed from a seasonal perspective, there are values that fluctuate between 0.5 and 0.6 for summer or winter and 0.86 and 0.81 for spring and autumn. The NSE values are high for models tested from an annual perspective, but for the seasons of spring and summer they are below 0.5 for the models HS_O, PMT_O2T, PMT_C2T and PMT_OUT. The fluctuations between models with an annual perspective of the RMSE and MAE were greater than if those models were compared from a seasonal perspective. During the winter none of the models showed good performance, with values of R²>0.59, NSE > 0.58 and RMSE (%) > 30 %. From a practical point of view, in the management of irrigated crops, winter is a season where crop water needs are minimal, with daily average values of ET_o around 1 mm due to low temperatures, radiation and VPD. The model that showed the best performance was PMT_CUH, followed by PMT_OUH and HS_C for annual and seasonal criteria. PMT_OUH is slightly less robust than PMT_CUH during the maximum radiation periods of spring and summer, since PMT_CHU performs the k_RS calibration. The performance of the HS_C model is better in the spring period, which is similar to PMT_CHU. The spatial distribution of MAEs in the basin shows that it is highly dependent on wind speeds, obtaining greater errors in areas with winds greater than 2.8 m s⁻¹ (east of the basin) and lower than 1.3 m s⁻¹ (south–southwest of the basin). This information of the tested models at different temporal and spatial scales can be very useful for adopting appropriate measures for efficient water management under the limitation of agrometeorological data and under the recent increments of dry periods in this basin. It is necessary to consider that these studies are carried out at a local scale, and in many cases the extrapolation of the results at a global scale is complicated. Future studies should be carried out in this way from a monthly point of view, since there may be high variability within the seasons.