Variations and Sensitivity Analysis of Reference Evapotranspiration During 2010-2019 in the Zhangye Farmland Oasis, Northwest China

Reference evapotranspiration (ET 0 ) is an important parameter for agricultural water management in the arid Zhangye farmland oasis. However, the ET 0 variations in this oasis over the last decade and meteorological forcings of these variations are unknown. This study investigated the ET 0 variations during 2010-2019 in this oasis using the FAO-56 Penman-Monteith (PM) and Hargreaves equations. Results showed that the ET 0 features daily and monthly variations with peak values in mid-July and an annual cycle. Although the estimated ET 0 series based on the two equations have high correlations in the time domain, the Hargreaves equation always underestimates the ET 0 compared to the PM equation. The yearly ET 0 showed statistically significant increasing trends (90% significance level) during 2010-2019, while statistically significant increasing trends in monthly ET 0 are found only in March and November. Increasing trends reflected in monthly and yearly ET 0 are mainly attributed to the increasing maximum temperature and sunshine duration and decreasing relative humidity. Sensitivity analysis demonstrated that the meteorological factor to which the ET 0 is most sensitive varies with time scale and equation. Moreover, regression equations used to correct the underestimation associated with the Hargreaves equation for estimating ET 0 in the Zhangye farmland oasis also were constructed.


Introduction
Reference evapotranspiration (ET0) is an important component of the regional hydrological cycle and is closely related to the biomass and productivity of terrestrial ecosystems [1][2][3] . ET0 is also one of the most important hydrological parameters for agricultural water resources management and estimating actual evapotranspiration 4,5 . Under the background of rapid global warming, studying ET0 variability at different spatio-temporal scales can help to maintain regional agricultural production and protect the regional eco-environment 3 . In addition, studies indicated that ET0 variations feature high spatial heterogeneity, especially in China's complex geographical and climatic conditions. For example, during the last decade, significant increasing trends in annual ET0 have occurred in the hilly regions of Southern China 6 , the upper Yangtze River basin 7 , the Pearl River Delta 8 , the non-monsoon region of China 9 , the Huang-Huai-Hai River Basin 10 and the Indo-China Peninsula 11 . In other regions, such as in Heilongjiang Province 12 and the central arid zone of Ningxia Province 13 , annual ET0 variations have experienced significant decreasing trends. In the Tibetan Plateau and the Loess Plateau, ET0 variations have been found to have more complex trends than other regions of China, namely both increasing and decreasing trends in the ET0 were observed in these regions 3,14 . The Zhangye farmland oasis is located in the Hexi Corridor (102°07'E-103°46'E, 36°31'N-37°55'N), which is a typical arid region in Northwest China. In addition, the Zhangye farmland oasis is the largest corn seed production base in China and uses a large amount of water of the Heihe River during the growing season of corn, giving rise to a conflict between the water requirements of natural ecosystems and agricultural irrigation. Therefore, an understanding of the ET0 variations during the last decade in this area is important for local agricultural production and oasis eco-environment protection. However, it is still unclear how ET0 has varied in the Zhangye farmland oasis during the last decade (i.e., [2010][2011][2012][2013][2014][2015][2016][2017][2018][2019] and what are the main meteorological factors driving these variations. Accordingly, the first aim of this study is to investigate the ET0 variations in the Zhangye farmland oasis during the last decade. Direct measurements of ET0 are time-consuming, cumbersome, and expensive [15][16][17] . Thus, ET0 usually is estimated using empirical or hybrid (energy and mass transfer) methods. Commonly-used mathematical equations for estimating ET0 include the Thornthwaite equation 18 , the Turc equation 19 , the Jensen-Haise equation 20 , the Priestley-Taylor equation 21 23 , as well as several variations of the PM equation. These equations have been widely compared with each other and applied in different climatic regions of the world [24][25][26][27][28] . Previous studies have shown that the PM equation is closest to actual evapotranspiration and is better than other ET0 estimation equations in many areas with different climatic conditions if all required meteorological factors can be fully satisfied. Usually, the PM equation is accepted as the standard method for estimating the ET0 and hence it has been adopted as a reference in many studies evaluating the performances of other ET0 estimation equations 17,29 . However, the PM equation requires the full availability of meteorological data including temperature, wind speed, relative humidity and sunshine duration. Globally, there are not many meteorological stations observing all these meteorological parameters. Consequently, the applicability of the PM equation is limited over areas with limited or unavailable meteorological data. The Hargreaves equation is one of the simplest equations used to estimate ET0, previous studies showed that the Hargreaves equation was on par with the PM equation in terms of accuracy for estimating ET0 in the arid and semi-arid climatic conditions 24,30 . In addition, the Thornthwaite equation usually provides very poor ET0 estimates compared to other equations 29 . The Turc equation is one of the most accurate empirical equations but it usually is used to estimate ET0 under humid conditions 31 . The Priestley-Taylor equation is a simplified version of the Penman equation that can be used when for wet surface areas 21,29 . Based on these considerations, the Hargreaves equation is a more applicable ET0 estimation equation in the arid Zhangye farmland oasis. Accordingly, the second aim of this study is to compare the ET0 variations estimated by the Hargreaves equation to those estimated by the PM equation.
ET0 is jointly affected by several meteorological factors, such as air temperature, wind speed, air humidity, solar radiation and sunshine duration. Different meteorological factors have different effects on ET0 variations [32][33][34] . Previous studies also showed that the sensitivities of ET0 variations to meteorological factors varied in space and time [34][35][36][37][38][39] . In order to understand the mechanisms of the regional ET0 variations caused by variability in meteorological factors, researchers have studied the sensitivities of ET0 to main meteorological factors using a variety of sensitivity analysis methods in different regions 16,36,38,[40][41][42][43] . However, at present, there is still a lack of systematic research on which meteorological factors have the greatest impact on ET0 in the Zhangye farmland during the last decade at different time scales. Therefore, the third aim of this study is to investigate the sensitivities of the ET0 variations estimated by the PM equation and the Hargreaves equation to main meteorological factors (e.g., maximum temperature, Tmax; minimum temperature, Tmin; mean temperature, Tmean; mean relative humidity, RHU; mean wind speed, WSP; mean sunshine duration, SSD) in the Zhangye farmland oasis during the last decade on daily and monthly time scales.

Variations in ET0
The daily, monthly and yearly ET0 series, which were estimated by the PM and Hargreaves equations based on six meteorological data (see Methods), are presented in Figs Fig. 1 demonstrates that the ET0 values estimated by both of the equations feature apparently daily and monthly variations and an annual cycle. The correlation coefficient between the two daily ET0 series estimated by the two equations is 0.95 (p < 0.01) and correlation coefficients between the monthly ET0 series estimated by the two equations from January to December are 0.85, 0.70, 0.94, 0.91, 0.69, 0.71, 0.89, 0.92, 0.92, 0.71, 0.56 and 0.89 (p < 0.01), respectively, indicating high correlations in the time domain between the two daily and monthly ET0 series respectively estimated by the two equations in the Zhangye farmland oasis. In addition, the ET0 increases from January to July, with peak values in mid-July, and then decreases from July to December (see Fig. 2a). Note that the estimated ET0 values using the PM equation are always greater than those estimated by the Hargreaves equation on daily, monthly and yearly time scales (see Figs. 1-3). Large differences between the two estimated ET0 series are found in March, April and May. Compared with the PM equation, the main reason for the Hargreaves equation underestimating ET0 may be that the PM equation considers the effects of WSP, RHU and SSD on the evapotranspiration. For example, the water vapor carried by winds is included in the total evapotranspiration estimated by the PM equation, and the role of sunshine in increasing plant transpiration is also considered by the PM equation. In addition, the greater the RHU of the atmosphere, the weaker the plant transpiration; on the contrary, the lower the RHU of the atmosphere, the faster the plant transpiration rate. However, the Hargreaves equation does not take into account the effects of these meteorological factors on evapotranspiration and hence underestimates the actual evapotranspiration compared to the PM equation on different time scales.  This study calculated the 10-years (i.e., 2010-2019) linear trend in monthly ET0 estimated by the two equations using the linear regression method 44 (see Fig. 3). Statistically significant (90% confidence level, hereinafter) increasing trends in monthly ET0 were found in March for both of the two estimated ET0 series (e.g., the trend is 4.12 mm/year for the PM equation and trend is 1.80 mm/year for the Hargreaves equation) and in November only for the PM equation (e.g., the trend is 1.51 mm/year). There is no statistically significant trend in the estimated monthly ET0 in other months for both of the equations during 2010-2019. Furthermore, there are statistically significant increasing trends in the yearly ET0 estimated by the PM equation (i.e., 17.28 mm/year) and the Hargreaves equations (i.e., 3.09 mm/year) during 2010-2019 (Fig. 2b). The reasons for these trends reflected in monthly and yearly ET0 series during 2010-2019 are analyzed in the discussion section.

Sensitivity analysis of ET0
The daily and monthly sensitivity coefficients (SCs) during 2010-2019 were calculated using the sensitivity analysis method (see Methods). Net radiation (Rn) is a necessary parameter in the PM equation. In this study, however, it was theoretical estimated based on the observed Tmax, Tmin and SSD (see Methods) rather than actual observations like Tmax, Tmin Tmean, RHU, WSP and SSD. Therefore, in this study, the sensitivity of ET0 to Rn was not analyzed. All of the six observed meteorological parameters (i.e., Tmax, Tmin Tmean, RHU, WSP and SSD) were considered in the PM equation, while only the temperature parameters (i.e., Tmax, Tmin Tmean) were considered in the Hargreaves equation. The variations of the daily and monthly SCs are presented in Fig. 4 and Fig. 5 and the corresponding statistics of the SCs are summarized in Table S1 and  Table S2     The monthly SC-PM-Tmax increases from January to May, reaches a peak in May, and then decreases from May to December (Fig. 5a). The monthly SC-PM-SSD increases from January to July, reaches a peak in July, and then decreases from July to December (Fig. 5f). The monthly SC-PM-WSP has a minimum value in July and two maximum values in May and November, respectively (Fig. 5e). In addition, the monthly SC-PM-Tmin indicates that the monthly ET0 is positively sensitive to Tmin in May, June, July, August, September and October while negatively sensitive to Tmin in January, February, March, November and December in the PM equation (Fig. 5b). The monthly SC-PM-Tmean is negative in all months except September; the minimum is found in January. The absolute values of the monthly SC-PM-Tmean are very small compared to other meteorological factors in all months (Fig. 5c). For the RHU, the monthly SC-PM-RHU is negative in all months, it increases from January to May, reaches a maximum in May, and then decreases from May to December (Fig. 5d). Moreover, the absolute values of the monthly SCs indicate that the monthly ET0 is more sensitive to RHU (i.e., mean = 0.51) than Tmax (i.e., mean = 0.  Fig. 5g indicates that in the Hargreaves equation the monthly ET0 is positively sensitive to Tmax in all months except in January. The monthly SC-H-Tmax increases from January to July, reaches a peak in July, and then decreases from July to December. In addition, the monthly ET0 is positively sensitive to Tmin in January, February, March, November and December and negatively sensitive to Tmin in May, June, July, August, September and October, with a dip in July in the Hargreaves equation (Fig. 5h). The monthly ET0 is positively sensitive to Tmean in March to October and negatively sensitive to Tmean in January, February, November and December in the Hargreaves equation (Fig. 5i). The monthly SC-H-Tmin and SC-H-Tmean show almost opposite phases by comparing Fig. 5h and Fig. 5i. The means of the absolute values of the monthly SCs indicate that the monthly ET0 is most sensitive to Tmax (i.e., mean = 0.56), followed by Tmin (i.e., mean = 0.50) and Tmean (i.e., mean = 0.41) in the Hargreaves equation.

Discussion
Since the PM equation takes into account more meteorological parameters affecting evapotranspiration (e.g., wind, air humidity and sunshine) and physical processes of crop evapotranspiration compared to other ET0 estimation equations 23 , previous studies confirmed that ET0 estimated by the PM equation is closest to actual evapotranspiration. In this study, ET0 estimated by the Hargreaves equation is lower than that estimated by the PM equation at different time scales. This underestimation associated with the Hargreaves equation is more apparent in March, April and May than those in other months ( Fig. 2a and Fig. 3). This is, possibly, because the ET0 is positively sensitive to the WSP and negatively to the RHU in the PM equation, and the WSP values in April and May are apparently larger than those in other months and the RHU values in March, April and May are apparently smaller than those in other months during 2010-2019 in the Zhangye farmland oasis (see Fig. S1). In addition, the Zhangye farmland oasis is a typical arid area in which the scarcity of meteorological data is still apparent compared to the Eastern China, thus, the Hargreaves equation has a great application potential in this area due to its inherent advantages of estimating ET0 in arid areas and the simplicity of data requirements 22,24,30 . However, results demonstrate that there are apparent deviations in the time domain between the two daily and monthly ET0 series respectively estimated by the two equations in the Zhangye farmland oasis though there are high correlations between the two ET0 series. Therefore, it is necessary to correct the coefficients of the Hargreaves equation and further to improve its accuracy for estimating ET0 in the Zhangye farmland oasis. In order to achieve this goal, this study established linear regression equations based on the estimated daily and monthly ET0 series. In these regression equations, the ET0 estimated by the PM equation is dependent variable and the ET0 estimated by the Hargreaves equation is independent variable. The established regression equation on daily ET0 is presented in Eq. (1). The established regression equations on monthly ET0 are presented in Table 1. When using the Hargreaves equation to estimate ET0 in other areas of the Zhangye farmland oasis without fully observed meteorological data, these regression equations can be used to correct the ET0 estimated by the Hargreaves equation so that the ET0 estimated by the Hargreaves equation is closer to those estimated by the PM equation.  Statistically significant increasing trends in monthly ET0 are found only in March for both of the two ET0 series (Fig. 3c). The reasons behind this result may be that the monthly Tmax and SSD show statistically significant increasing trends only in March during 2010-2019, and the monthly ET0 is positively sensitive to Tmax and SSD in March in the PM equation (Fig.  5a,f). The monthly RHU in March during 2010-2019 showed a statistically significant decreasing trend (Fig. S5) and the monthly ET0 is negatively sensitive to the RHU in the PM equation. In addition, although the monthly Tmean showed a statistically significant increasing trend in March during 2010-2019 and the monthly ET0 is negatively sensitive to the Tmean in the PM equation, the monthly SC-PM-Tmean is apparently smaller than the monthly SC-PM-Tmax and SC-PM-RUH, and thus result in an increasing trend in the monthly ET0 in March during 2010-2019 under the dominating influences of the Tmax and RHU. Moreover, the monthly ET0 in March is positively sensitive to the Tmax, Tmin and Tmean in the Hargreaves equation (Fig.  5g,h,i). The monthly Tmax and Tmean showed statistically significant increasing trends in March and the monthly Tmin showed no statistically significant trend in March during 2010-2019. Therefore, under the joint influences of Tmax and Tmean, the monthly ET0 in March estimated by the Hargreaves equation showed an increasing trend during 2010-2019. Furthermore, the monthly ET0 did not show statistically significant trends in other months, which may be due to the fact that there is no statistically significant trend in the six meteorological parameters in most of those months during 2010-2019 (Figs. S2-S7), especially for the two most sensitive parameters (i.e., Tmax and RHU) (Fig. S2 and Fig. S5). Finally, the two yearly ET0 series also showed statistically significant increasing trends. This may result from the joint influences of statistically significant increasing yearly Tmax, WSP and SSD and decreasing yearly RHU during 2010-2019 in the Zhangye farmland oasis (Fig. S8).
The common meteorological parameters contained in the PM equation and the Hargreaves equation are Tmax, Tmin and Tmean. However, although the SC-PM-Tmax is very similar to the SC-H-Tmax in term of monthly variation, the sensitivities of ET0 to Tmean and Tmin in the PM equation are apparently different from those in the Hargreaves equation. For example, from the Fig. 5, the monthly SC-PM-Tmean is negative in all months while the monthly SC-H-Tmean is positive in most months except January, February, November and December. The monthly SC-PM-Tmin is negative in January, February, March, November and December and is positive from April to October. While the variation of the monthly SC-H-Tmin is exactly opposite to that of the monthly SC-PM-Tmin, namely the monthly SC-H-Tmin is positive in January, February, March, November and December and is negative from April to October. This different sensitivities of ET0 to the same temperature variable in various equations implies that temperature has different or even opposite effects on estimating ET0 depending on the equation, because each ET0 method is based on a different physical mechanism. Therefore, when explaining the influence of temperature on estimating ET0, we may not rely on the results of a single equation.
In addition, although the sensitivity analysis method (i.e., the SC method) used in this study is a widely-used method to investigate the response of ET0 to meteorological parameters, it has some drawbacks. Specifically, the prerequisite of the SC method is that there is only one variable can be changed each time, while the other variables remain unchanged. However, in fact, there is usually a covariant relationship between meteorological variables. For example, a change of the mean temperature is often accompanied by simultaneous changes of the maximum temperature and the minimum temperature. In addition, there is a strong connection between relative humidity and air temperature. Increasing (decreasing) air temperature will lead to a decrease (increase) in relative humidity. Therefore, theoretically, the use of a global sensitivity analysis method 47 to study the response of the ET0 to the global changes of multiple meteorological factors in the Zhangye farmland oasis is more reasonable than the use of the SC method. However, this is an ongoing work which will be presented in future.
In summary, agricultural irrigation in the Zhangye farmland oasis requires a large amount of water. ET0 is a key variable in agricultural water management, farmland ecosystems and the hydrological cycle. Thus, it is of great significance to study ET0 variations and meteorological forcings behind these variations in the Zhangye farmland oasis. In this study, the variations of daily, monthly and yearly ET0 estimated by the PM and Hargreaves equations during 2010-2019 in the Zhangye farmland oasis were analyzed and compared with each other. In addition, the sensitivities of ET0 to six major meteorological parameters effecting crop evapotranspiration were also investigated on daily and monthly time scales using the SC method. Particularly, regression equations used to correct the underestimation associated with the Hargreaves equation for estimating ET0 in the Zhangye farmland oasis are also constructed. All findings provide a scientific basis for understanding ET0 variations and the sensitivity of the ET0 to main meteorological parameters in the Zhangye farmland oasis during the last decade, and can help improve agricultural water resources management and oasis environmental protection in the Zhangye farmland oasis.

Meteorological data
Daily data for the six meteorological parameters were downloaded from the China Meteorological Data Network (http://data.cma.cn/data/). These meteorological data were observed by the Zhangye national meteorological reference station (station number: 52652; longitude: 100.17 E; latitude: 39.05 N; altitude: 1461.10 m). Daily Tmean was calculated according to the mean of the temperature observed at 02:00, 08:00, 14:00 and 20:00 (Beijing time) rather than the mean of daily Tmax and Tmin. The daily, monthly and yearly ET0 from January/1/2010 to December/31/2019 were calculated based on these meteorological data.

The FAO-56 Penman-Monteith equation
The FAO-56 Penman-Monteith equation (hereinafter PM equation) is based on the energy balance and water vapor diffusion theory and overcomes several shortcomings associated with previous Penman equations 23 . Because the PM equation is currently the most accurate method for estimating ET0 in both humid and arid climatic conditions 24,26,29 , it was recommended by the FAO as the only method to estimate the ET0 when all the required meteorological data are available. In addition, the PM equation also has been considered as a standard method to evaluate other ET0 estimation equations 17 Where: ET0: reference evapotranspiration (unit: mm·d -1 ); es: saturation vapour pressure (unit: kPa); ea: actual vapour pressure (unit: kPa); ∆: slope vapour pressure curve (unit: kPa·℃ -1 ); γ: psychrometric constant (unit: kPa·℃ -1 ); u2: wind speed at 2 m height (unit: m·s -1 ); Tmean: mean daily air temperature at 2-m height (unit: ℃); G: soil heat flux (unit: MJ·m-2·day -1 ); Rn: net radiation at the crop surface (unit: MJ·m-2·day -1 ); All parameters in the above equation can be theoretically estimated according to the FAO methods provided that the Tmax, Tmin, latitude and altitude data are known 23 . The value of soil heat flux G is very small compared to Rn; in the case of time steps greater than or equal to 1 day, the value of G is approximately 0 and may often be ignored 3,23 . In this study, the daily radiation data of the Zhangye meteorological station from January/1/2010 to December/31/2019 did not form a continuous time series and also did not meet the criterion for time series interpolation because of a lot of missing values in the raw observation series. Therefore, this study used the theoretically estimated Rn, which was calculated based on the observed Tmax, Tmin and SSD data, the meteorological location and astronomical parameters 23 . The process of estimating Rn was as follows. First, this study calculated the total solar radiation reaching the ground (Rs, also known as solar short-wave radiation) based on the observed SSD and the Angstrom-Prescott equation 23  Ra can be calculated based on the meteorological location and astronomical parameters 23 . On the basis of the estimated Rs and Ra, the net radiation Rn can be gained by calculating the difference between the incoming net shortwave and the net outgoing longwave 23 . Rn is normally positive during the daytime and negative during the nighttime.

The Hargreaves equation
The Hargreaves equation was established for estimating evapotranspiration in the western arid regions of the USA 22 . The input variables required by this equation are very simple. All we need to know are Tmax, Tmin, Tmean and Ra. Tmean can be replaced by the mean of Tmax and Tmin in case of missing Tmean data 29 . All meteorological stations measure Tmax, Tmin, and Tmean, so they are easy to obtain. Ra can be computed based on the meteorological station location and astronomical parameters 23 Where: Tmax: daily maximum temperature at 2-m height (unit:℃); Tmin: daily minimum temperature at 2-m height (unit: ℃); Tmean: mean daily temperature at 2-m height (unit: ℃); Ra: extraterrestrial (solar) radiation (unit: MJ·m-2·day -1 );

Sensitivity analysis method
Sensitivity analysis is used to identify the extent to which the variability in the input values for a given variable (e.g., air temperature, wind speed, air humidity and sunshine duration in this study) will impact the result (e.g., ET0 in this study) in a mathematical model, and further calculate the degree of influence of these input variables on the result 45,46 . Currently, several sensitivity analysis methods have been proposed, of which the sensitivity coefficient (SC) method is a frequently-used method in earth sciences 47 , especially in studying the sensitivity of the ET0 to meteorological factors 16,36,43,[48][49][50] . The SC method was first proposed by McCuen 45 , and the SC is calculated by using the partial derivative of ET0 to various meteorological factors, namely the ratio between the relative variation of ET0 and the relative variation of a meteorological factor. The SC method varies one parameter (e.g., adding perturbation to this parameter) while keeping other parameters constant, and observes the variations of the model outputs before and after adding this perturbation to the parameter. However, currently, there is no universal rule for determining how much the perturbation will be added to the parameter. Therefore, in the sensitivity analysis of this study, specifically, the perturbation is β times the value of a meteorological variable, where β is a constant randomly sampled from a uniform distribution (0~20%) 16 . The form of the SC method is as follows: Where: : SC of the i th variable, unitless; Vi: ith variable; 0 : change amount of ET0; ∂V i : change amount in the i th variable; β: perturbation coefficient. A positive SC indicates that ET0 increases with the increase of meteorological parameters, and vice versa. Moreover, the larger the absolute value of SC associated with one meteorological parameter, the higher the sensitivity of ET0 to this meteorological parameter, and vice versa.