Estimating reference crop evapotranspiration with elevation based on an improved HS model

Crop evapotranspiration (ET), the only term that appears both in water balance equations and surface energy balance equations (Xu & Singh 2005), is an essential component of the hydrologic cycle. Therefore, crop ET plays a key role in applying rainfall-runoff and ecological models, calculating regional water balance, and determining irrigation water demands. The first step of crop evapotranspiration estimation is estimating the reference crop evapotranspiration (ET₀). Inaccurate ET₀ estimates could lead to error propagation in the model and finally unreliable groundwater recharge estimates (Allen et al. 1998). Although there are many methods to estimate ET₀, each of them is generally limited by geographical and meteorological conditions (Liu & Pereira 2001; Irmak et al. 2003). The applicability of the Penman-Monteith model (PM) has been shown in many regions of different climate types (Jabloum & Sahli 2008; Tabari et al. 2013). Thus, the Food and Agriculture Organization (FAO) has adopted the PM equation as the standard method to estimate ET₀. The main shortcoming of the PM model is its requirement of some hard-to-obtain meteorological factors (Valiantzas 2015). It is very difficult to obtain continuous and comprehensive data in many cases. Furthermore, installation and maintenance of meteorological stations in high-elevation terrain is very expensive and complicated. Therefore, it is crucial to find a simple and accurate method to estimate ET₀ for high-elevation areas.

The existing empirical equations to estimate ET₀ with limited data were categorized as temperature-based, radiation-based, evaporation-based, material conversion-based and mixed-type (Tabari et al. 2013). The simplified temperature-based models, which only require easy-to-access data, have been widely used (Almorox et al. 2015). The temperature-based Hargreaves-Samani model (HS) (Hargreaves & Samani 1985) calculated solar radiation using the highest temperature, lowest temperature and external radiation. Almorox et al. (2015) found that the HS model was the simplest and most accurate temperature-based model to use in arid and semi-arid areas after studying data from 4,362 meteorological stations to estimate ET₀ using 11 temperature-based formulas. Er-Raki et al. (2010) also showed that the HS model was the most accurate temperature-based model in semi-arid regions of the Americas. Some scholars analyzed the applicability of commonly-used ET₀ empirical models in various climate zones throughout China. Hu et al. (2011) evaluated the applicability of the HS model at monthly scale in different climate zones of China. In addition, Wang et al. (2008) validated the applicability of the HS model in some arid areas and semi-arid areas of China.

Although the HS model has been popular for estimating ET₀ due to its simplicity, there is evidence that ET₀ is underestimated in arid regions and overestimated in humid areas by the HS model. Martinez & Thepadia (2010) found that ET₀ was overestimated by the HS model after an analysis of data from 72 stations in Florida. Yoder et al. (2005) estimated the daily and monthly ET₀ in the Cumberland Plateau in the southeastern humid region of the United States and found that the HS model was more suitable for estimating the monthly ET₀ compared with the PM equation.

Many ET₀ empirical models were only applied in certain areas (Jabloum & Sahli 2008), which means that the parameters in the empirical models had to be recalibrated for new areas. For example, Allen et al. (1998) pointed out that the value of the temperature coefficient in the HS model was not reasonable at altitudes over 1,500 m. Annandale et al. (2002) stated that the temperature coefficient and temperature index in the HS model should be modified to account for atmospheric pressures in different areas, and the atmospheric pressure decreased with altitude (Yang & Zhang 2009).

However, elevation factors, which were overlooked before, can be easily obtained and are closely related to the hygrometer constant (γ), the zenith radiation (R_a, solar radiation absorbed by the top of the earth's atmosphere) and the net solar radiation (R_n, the energy absorbed by the Earth's surface) (Allen et al. 1998). Existing studies on simplified models mostly focus on temperature, relative humidity, wind speed, sunshine hours, rainfall and other meteorological factors. The main objective of this study is to provide a simple and accurate model to estimate ET₀ based on the HS model, using elevation information.

Tibet, located in the southwest of China and known as the ‘roof of the world’, has an average altitude of 4,000 m covering an area of 1.23 million km². The high altitude leads to unique climate characteristics: low oxygen, low pressure (less than 2/3 atmospheric pressure than at sea level), long sunshine duration (more than 3,000 hours per year) and strong radiation (annual solar radiation ranging from 6,000 to 8,000 MJ m⁻²). The humidity and temperature vary widely, and there are distinctive dry and wet seasons in the high-altitude areas of Tibet (above 2,000 m) due to the frequent cold-heat exchanges in the near surface layer. According to the characteristics of Tibetan geomorphology and climate, this study fully considered the altitude variations of Tibet (the overall trend is higher in the west and lower in the east). Twenty representative Tibetan sites with altitude ranging from 2,000 m to 5,000 m were selected (Figure 1) for the study. Their geographical data are shown in Table 1.

The daily meteorological data including the minimum and maximum temperatures, humidity, sunshine duration, wind speed, rainfall and elevation from the 20 sites over 30 years were obtained from the National Meteorological Information Center (NMIC, http://data.cma.cn). The daily and monthly ET₀ were calculated using the Penman-Monteith equation. Then, the 30-year meteorological data were divided into two sub sets: the training dataset – the daily weather data of the 20 sites from 1981 to 2000 (n = 135144, no interpolation); and the testing dataset – the daily weather data of the 20 sites from 2001 to 2010 (n = 73,040, no interpolation). An elevation-adjusted HS model was developed using the training dataset, and validated using the testing dataset.

The HS equation was tested based on lysimeter measurements over an 8-year period (Hargreaves & Samani 1985). The HS equation can utilize atmospheric temperature differences (T_max–T_min, indicating how much of the zenith radiation can reach the earth's surface) to estimate solar radiation. Many studies showed that the HS equation gave a reasonable ET₀ throughout the world in the condition of some missing meteorological data. However, the empirical coefficients, such as the temperature coefficient (0.0023), temperature constant (−17.8) and temperature index (0.5), were obtained through regression analysis of many observations collected from around the world over many years. Thus, the parameters in the HS equation should be modified before it is used in a new area.

The analysis of the variation of ET₀ with elevation and temperature showed that the annual ET₀ increased gradually with increasing elevation (decreasing temperature) in areas below 4,000 m, but decreased with increasing elevation (decreasing temperature) in areas above 4,000 m (Figure 2 and Table 1). The annual highest cumulative value of ET₀ (about 1,280 mm/year on average) was found in the Lazi area (about 4,000 m ASL), and the lowest cumulative value of ET₀ (about 950 mm/year on average) was found in Chayu (the lowest altitude site selected in this study, about 2,000 m ASL) and the Ando area (the highest altitude site selected in this study, about 5,000 m ASL). In summary, elevation and temperature were very important in ET₀ estimation, and this was why we modified the temperature constant and the elevation function to improve the HS model.

The ET₀ and the annual average temperature were abnormal at the Damxung area (4,200 m above sea level), which has a lower annual ET₀ compared to those of Tingri (4,300 m above the sea level) and Lazi (4,000 m above sea level) (Figure 2). This was due to the unique terrain and landscape of Damxung, which led to lower annual average temperature (1.81 °C) in Damxung than those of Tingri (3.12 °C) and Lazi (6.54 °C), and further caused the lower annual cumulative evapotranspiration in the Damxung area.

Principal component analysis (PCA) is a statistical method that is used to extract the most important information (i.e. variance) in a dataset (Wang 2015; Wang et al. 2018). The new variables obtained from the PCA are called principal components (PCs), which are linear combinations of the original variables, and each PC is orthogonal to other PCs.

When calculating ET₀ using the PM equation, seven meteorological parameters were used, including the daily maximum temperature, daily minimum temperature, daily average temperature, daily average wind speed, daily average relative humidity and daily sunshine duration, all of which were obtained from the NMIC. The non-continuous data series contained geospatial data in three dimensions: longitude, latitude and altitude. Due to the large number of continuous meteorological parameters needed for this complex calculation equation, a simpler alternative is often necessary when data are sparse. Statistical Product and Service Solutions (SPSS) software (IBM, USA) was used to analyze the PCs based on the geospatial data and 30-year meteorological data from the 20 study sites in Tibet. The original data were automatically standardized to eliminate the influence of the different variance scales of the variables.

The first three PCs (Table 2) had eigenvalues larger than 1, and the total variance explained by them was nearly 80%. Therefore, the first three PCs could be used to illustrate the main information contained in the original nine variables. The correlations between the first three PCs and the original variables were used to find the relationship between the PCs and the original meteorological indexes or geographic information (Table 3). There were significant correlations between the first principal component and four of the initial variables: the daily minimum temperature, daily mean temperature, daily maximum temperature and elevation. The second and third PCs only had significant correlation with sunshine duration and latitude, respectively.

The PCA analysis further confirmed that temperature and elevation were important indicators of ET₀.

FAO recommended that the temperature constant should be −17.8, when the unit of R_a was measured in mm d⁻¹. However, T_mean is often lower than −17.8 °C in high-altitude areas, especially in areas above 3,500 m in January and December. This led to negative ET₀ values calculated by the HS model. This result is contrary to common sense (actual evapotranspiration cannot be negative). Therefore, the method to determine the temperature constant needs to be improved.

The lowest daily T_mean of the 20 representative sites (Table 1, elevation ranging from 2,000 to 5,000 m) from 1960 to 2015 was −36.6 °C, which was observed at the Gerze region (4,414.9 m altitude) on December 26th, 1987. Therefore, the temperature constant was set to −36.6 to avoid negative values of ET₀ calculated using the new model.

The new HS-E model was validated using the testing dataset (the meteorological data from the 20 typical sites from 2001 to 2010) at different time scales.

The MRE is a percentage and is used to characterize the accuracy and error scope of the prediction; the RMSE is used to measure the deviation between the observed values and true values. The closer the MRE and RMSE are to zero, the higher reliability of the model. The closer the NES value is to 1, the higher the reliability of the model.

Daily ET₀ (n = 3,653 × 20) was estimated by the HS model and the HS-E model using the testing dataset. Then RMSE, MRE and NSE were calculated, using the results from PM equation as the true values, to evaluate the performance of the HS model and the HS-E model (Table 4).

First, the NSE value of the HS-E model was between 0.60 and 0.88 with an average of 0.79; while the NSE value of the HS model was between 0.35 and 0.84 with an average of 0.69 (Table 4). Thus, with NSE closer to 1, the simulation results of the HS-E model were more accurate than those of the HS model. Second, the RMSE of the HS-E model was between 0.40 and 0.80 mm day⁻¹, with an average of 0.56 mm day⁻¹; while the RMSE of the HS model was between 0.51 and 0.89 mm day⁻¹, with an average of 0.65 mm day⁻¹ (Table 4). Third, the MRE of the HS-E model was between 1.13% and 18.72% with an average (the average MRE absolute value) of 9.80%, while the MRE of the HS-E model was between −31.64% and 18.41% (the average absolute value was 15.13%) (Table 4). Therefore, with the RMSE and MRE closer to 0, the values of ET₀ calculated by the HS-E model had a smaller deviation compared with the true values. Fourth, when plotting against the ET_{0 PM}, the slope of ET_{0 HS-E} ranged from 1.03 to 1.24 with an average of 1.09; while the ET_{0 HS} ranged from 0.63 to 1.21. Therefore, the HS-E model was more stable and its estimation was closer to the true value. In summary, when estimating daily ET₀, the HS-E model avoided the occurrence of negative values and provided estimations closer to true values calculated by the PM equation, compared with the HS model.

Monthly ET₀ (n = 120 × 20) was also estimated by the HS model and the HS-E model using the testing dataset (Figure 4), and then the RMSE, MRE and NSE were calculated, using the results from PM equation as the true values, to evaluate the performance of the HS model and the HS-E model (Table 5). First, the NSE value of HS-E model for monthly ET₀ was between 0.77 and 0.98 with an average of 0.89; while the NSE of HS model for monthly ET₀ ranged from 0.06 to 0.93, with an average of 0.69 (Table 5). Therefore, the estimated results of the HS-E model were of higher quality than those of the HS model. Second, the average RMSE of 9.52 mm month⁻¹ for ET₀ obtained by the HS-E model was less than the average RMSE of 15.51 mm month⁻¹ for ET₀ obtained by the HS model. Third, the MRE of ET₀ according to HS-E model varied from −4.26% (Pulan) to 22.10% (Changdu), with an average (the average of MRE absolute values) of 7.21%; while the average MRE of the HS model was 15.51%. Therefore, with the RMSE and MRE closer to 0, the values of ET₀ calculated by the HS-E model had a smaller deviation compared with the true values.

On the other hand, first, the NSE of 0.89 for the monthly ET_{0 HSE} was greater than 0.79 at the daily scale. Second, the RMSE of 9.52 mm month⁻¹ for the monthly ET_{0 HSE} was lower than the 0.53 mm day⁻¹ at the daily scale. Third, the MRE of 7.21% for the monthly ET_{0 HSE} was less than the 7.44% at the daily scale in 20 different stations, meaning the monthly ET₀ calculated by the HS-E model was more accurate than that at daily scale. Therefore, the HS-E model has smaller errors with the increasing calculation period and more obvious advantages over the HS model.

So, the simulation quality of the HS-E model was much better than that of the HS model at monthly time scale.

Seasonal ET₀ (n = 40 × 20) was estimated by the HS model and the HS-E model using the testing dataset. Then RMSE, MRE and NSE were calculated, using the results from PM equation as the true values, to evaluate the performance of the HS model and the HS-E model (Figure 5). First, The NSE of HS-E model for the seasonal ET₀ was between 0.74 and 0.99 with an average of 0.90; while the NSE of HS model for seasonal ET₀ ranged from −0.09 to 0.94, with an average of 0.69 (Figure 5). Second, the average RMSE of HS-E model with 25.47 mm season⁻¹ was less than the average RMSE of HS model with 42.39 mm season⁻¹. Third, the MRE of HS-E model varied from −4.74% (Lasa) to 20.44% (Changdu), with an average (the average of MRE absolute values) of 7.0%; while the HS model had an average MRE of 15.43%. Therefore, the HS-E model was superior to the HS model.

The ET₀ values calculated by the HS-E model for seasonal data produced an average NSE of 0.90 at the seasonal scale, slightly greater than the 0.89 at the monthly scale and obviously greater than the 0.79 at the daily scale. An average RMSE of 25.47 mm season⁻¹ was found for the seasonal scale, slightly lower than the 9.52 mm season⁻¹ for the monthly scale and significantly lower than the 0.53 mm day⁻¹ for the daily scale. The seasonal ET₀ calculated by the HS-E model was slightly more reliable than calculated at the monthly scale, and both were much better than that at the daily scale. Therefore, the HS-E model has smaller errors with the increasing calculation period and with more obvious advantages over the HS model.

So, the simulation quality of the HS-E model was much better than that of the HS model at larger time scales.

Due to the difficulty in obtaining continuous climate data for ET₀ calculation by PM method, many studies were conducted to simplify the calculation method of ET₀ (Hargreaves & Allen 2003; Trajkovic 2005; Traore et al. 2010; Marti et al. 2011). It should be noted that the HS model used in this study is simpler compared with other existing studies. Many previous studies added continuous meteorological parameters, such as the relative humidity and sun hours, into the HS model to enhance its computational accuracy (Valiantzas 2015). However, increasing the complexity of the HS model also leads to higher demands for resources, including obtaining basic meteorological data and running a more complicated simulation (Ishak et al. 2010), which undermined the essence of HS model: using a simple temperature algorithm that only requires a few parameters.

The atmospheric temperature difference (T_max–T_min) was utilized as an indicator to express the amount of zenith radiation reaching the surface of the earth in the HS model. Ren et al. (2016a, 2016b) estimated the daily ET₀ for arid to sub-humid areas in Inner Mongolia, China base on HS model. However, due to the unique nature of Tibet, the temperature at high elevation was not as high as expected under strong radiation (China Tibet Information Center 2005). The lower than expected atmospheric temperature resulted from the thinner atmosphere and lower absorption of long-wave radiation from the ground at higher elevation. Meanwhile, there is less cloud cover at higher elevation areas. And the less cloud cover means less long-wave radiation absorption from the ground in the day time and weaker counter radiation of the cloud against the ground at night, both leading to poor insulation (Allen et al. 1998). All of the above reasons weakened the sensitivity of temperature differences (T_max–T_min) to the radiation. Therefore, the elevation factor was included to enhance the ability of the model to estimate radiation and improve the accuracy of the HS model.

This paper mainly focused on the high-elevation areas of Tibet, so validation of the HS-E equation was also limited to that region. However, adjusting the temperature constant and using the elevation factor, the proposed HS-E model has the potential to easily estimate ET₀ in a large variety of areas with easy to obtain parameters.

The standard PM method (adopted by FAO) requires a lot of continuous meteorological data that are difficult to obtain in areas of complex terrain and high elevation. And the popular HS model could not avoid negative ET₀ values due to the inappropriate temperature constant. We adjusted the temperature constant in HS model by taking the minimum of the annual average temperature in 50 years; and we included an elevation factor in the HS model to better relate temperature difference to radiation. The elevation factor was derived by principle component analysis with geospatial data in Tibet. The improved HS model, HS-E model, was calibrated using meteorological data from 1981 to 2000 at 20 representative sites in Tibet, and validated using meteorological data from 2001 to 2010 at the same sites (using the results from PM equation as the true values).

The HS-E model not only avoided the occurrence of negative values in the ET₀ estimation, but also enhanced the accuracy of the ET₀ estimation. The accuracy of the HS-E model was shown to be superior to that of the HS model under different time scales (daily , monthly and seasonal) using our testing dataset. The NSE of the HS-E model reached 0.78 (daily), 0.89 (monthly) and 0.90 (seasonal), showing the high reliability of the model. Meanwhile, by analyzing the results of the NSE, RMSE and MRE at different time scales, the monthly ET₀ was found to be closer to the PM standard ET₀ compared to the daily ET₀. The larger the time scale, the better the HS-E model performed. Overall, the HS-E model is feasible and convenient for use at high elevation areas. However, this study only used data from sites with elevations ranging from 2,000 to 5,000 m. In the future, it is recommended to consider the effect of elevation using more meteorological data at sites covering wider ranges of elevations.

This study was supported by the National Science Foundation of China (51579158, 51609154), State Key Laboratory of Simulation and Regulation of Water Cycle in River Basin (SKL2018ZY01) and scientific researching fund projects of China Institute of Water Resources and Hydropower Research (MK2017J02). The contributions of the editor and anonymous reviewers whose comments and suggestions significantly improved this article are also appreciated.