Differences in Reference Evapotranspiration Variation and Climate-Driven Patterns in Different Altitudes of the Qinghai–Tibet Plateau (1961–2017)

: Reference evapotranspiration ( ET 0 ) in the hydrological cycle is one of the processes that is significantly affected by climate change. The Qinghai–Tibet Plateau (QTP) is universally recognized as a region that is sensitive to climate change. In this study, an area elevation curve is used to divide the study area into three elevation zones: low (below 2800 m), medium (2800–3800 m) and high (3800–5000 m). The cumulative anomaly curve, Mann–Kendall test, moving t-test and Yamamoto test results show that a descending mutation occurred in the 1980s, and an ascending mutation occurred in 2005. Moreover, a delay effect on the descending mutation in addition to an enhancement effect on the ascending mutation of the annual ET 0 were coincident with the increasing altitude below 5000 m. The annual ET 0 series for the QTP and different elevation zones showed an increasing trend from 1961 to 2017 and increased more significantly with the increase in elevation. Path analysis showed that the climate-driven patterns in different elevation zones are quite different. How-ever, after the ascending mutations occurred in 2005, the maximum air temperature (Tmax) became the common dominant driving factor for the whole region and the three elevation zones.


Introduction
Rapid global climate change has had a vast impact on the ecological environment and biological living systems [1,2]. Changes in water resources and water demand are among the important signals of climate change [3][4][5][6]. As the most active factor of the water cycle, evapotranspiration is not only the unique junction of the water transport and energy cycle in the earth-atmospheric system but also a hydrological cycle process, which is significantly affected by climate change [7][8][9][10]. As a result of the unfeasibility of the direct prediction of evapotranspiration, people usually estimate the reference evapotranspiration (ET0) first when formulating crop irrigation plans [11]. ET0 is a metric that refers to the maximum evaporation that can be achieved by a fixed underlying surface when the water supply is sufficient, as defined by the Food and Agriculture Organization of the United Nations (FAO) [12]. ET0 is one of the important indicators to reflect the regional evapotranspiration capacity and crop water requirement [13]. The change of ET0 is only related to climate factors and location. Studying the changing trend of ET0 plays an essential role in the calculation of the water balance in a watershed when developing crop irrigation plans and drought warning [14][15][16][17].
Climate change has significantly changed the spatial-temporal patterns of reference evapotranspiration in many regions [18][19][20][21][22][23]. Since the 1960s, the increasing global temperature has resulted in an obvious change in the trend of reference evapotranspiration in different regions of the world [24][25][26][27]. In recent years, the annual series of ET0 showed a significant increasing trend in the Amazon basin, Kazakhstan, Romania and the Ecuador coast [24][25][26][27]; however, a decreased ET0 appears in many other regions around the world, such as the Qilian Mountains, the Huaihe River Basin, the upper reaches of the Yellow River, and a major part of Northern China [28][29][30][31][32][33][34][35]. The abovementioned phenomenon is called "the evaporation paradox". Xing et al. [36] found that this was due to the limitation of relative humidity. Furthermore, relative humidity is the main driving factor of ET0 change in Tunisia, the Loess Plateau, Sichuan Plateau, Guangxi Basin, Yunnan Guizhou Plateau and many other regions of China [37][38][39][40]. It has also been found that wind speed is the dominant variable in the Huaihe River Basin and Southwest China [34,41]. McVicar et al. [29] concluded that the decrease in ET0 was mainly related to aerodynamic factors (the global wide performance of wind speed decline) and relative humidity (the significant humidity increase in the warm season in Australia, Central and Southern Canada) at the global scale. In addition, radiation composition, temperature and pollutant concentration were also discovered to be key factors motivating the change of ET0 [42][43][44]. Because different regions respond to climate change quite differently, it is necessary to conduct regional research. Moreover, most of the relevant research at present focuses on the change trends and its impact factors, only a few people have studied the impact relationships between these factors.
Liu et al. [45] found that regions with a higher altitude are more sensitive to climate change. As the most sensitive and vulnerable region in the context of global environmental change, the warming rate of the Qinghai-Tibet Plateau (QTP) is significantly higher than other parts of the world during the same period [9,[45][46][47]. A large number of studies have shown that the dramatic warming of the QTP has led to the rapid degradation of glaciers and subsequently strengthened the hydrological cycle process and the hydraulic connectivity in this region [48][49][50][51]. Furthermore, Huang et al. [52] suggested that compared with human factors, climate factors have a greater impact on vegetation dynamics in the QTP. Climate warming leads to the degradation of alpine grassland and changes the soil nutrients as well as the components of the alpine grassland ecosystem, which has a strong impact on the productivity of the plant community and soil microbial community [53][54][55][56][57][58]. In addition to being a significant characteristic of climate change, ET0 is also an important index of crop water demand and a key referential factor of the hydrological cycle. Consequently, the analysis of the ET0 changing trend and climate-driven mode of the QTP is not only of great value as a reference for guiding regional agricultural activity and investigating the structure of the ecosystem, but also of great significance for monitoring the dynamic changes of the world's water resource reserves. However, most of the existing research is mainly based on the spatial scale of the whole region of the QTP, ignoring the huge altitude difference in the region.
The objectives of this study are as follows: (1) to identify the temporal variation trend and mutation year of ET0 at different elevations in the QTP over the past 57 years, (2) to analyze spatial characteristics of ET0 and its trend, and (3) to explore the driving pattern and influence process between climate factors in different elevation zones on ET0 change before and after mutation. The results can help people understand the response mechanism of ET0 to climate change, clarify the complex relationship between multiple climate factors and provide a reference for agricultural and livestock production and water resource management.

Study Area and Data
The Qinghai-Tibet Plateau (26°00′ N~39°47′ N, 73°19′ E~104°47′ E), which has a total area of approximately 2.5 million km 2 , involving 8 countries (Figure 1), with an average altitude of more than 4000 m, is the plateau with the highest altitude in the world and is known as the "roof of the world" [59]. Furthermore, with approximately 70% of the global alpine permafrost, the QTP is the source of some major Asian rivers, including the Indus, Ganges, Salween, Brahmaputra, Mekong, Irrawaddy, Yangtze and Yellow rivers, known as "the world water tower" [60][61][62]. Because of its unique geographical conditions, fragile ecological environment and sensitive climate change, the response of the QTP to climate change is indicative of other regions.
The datasets were obtained from the China Meteorological Data Service Center (CMDC) (available online at: http://data.cma.cn/, accessed on 13 April 2020), which have undergone quality control and outlier processing before release (available online at: http://data.cma.cn/, accessed on 13 April 2020), . About the few missing data, we interpolate it by the nearest neighbor method. We selected 46 meteorological stations with longterm, relatively continuous and representative data in the QTP ( Figure 1). The daily meteorological data from 1961 to 2017 contain eight climatic elements, including daily actual atmospheric pressure (AP), mean temperature (Tmean), maximum temperature (Tmax), minimum temperature (Tmin), precipitation (P), relative humidity (RH), sunshine duration (SD) and wind speed (U2). The altitude of the QTP ranges from 732 m to 8848.43 m (Figure 1). In this study, the area-elevation curve of the QTP was drawn (Figure 2), and 46 meteorological stations were assigned into three elevation zones (low, medium, and high) according to the curve to analyze the variation trend of ET0 in different elevation zones in the QTP. Immerzeel et al. [60] used this method to capture the spatial climate heterogeneity in the QTP.

FAO Penman-Monteith Formula
The Penman-Monteith (P-M) formula is recommended by FAO [12] as the sole standard method to calculate the ET0, and it is currently one of the most widely used, together with the pan evaporation method and some formulas, in the world [63,64]. In this paper, the P-M formula is used to calculate daily ET0 values: where ∆ is the slope of the relationship curve between saturated water pressure and temperature (kPa °C −1 ); is the net radiation on the crop surface (MJ m −2 ·day −1 ); G is the soil heat flux (MJ m −2 ·day −1 ); is the hygrometer constant (kPa °C −1 ); is the average air temperature (°C); is the wind speed at a height of 2 m above the ground (ms −1 ); is the air saturation vapor pressure (kPa); is the actual vapor pressure (kPa). The net radiation is calculated using the following formula: (60) [ sin sin + cos cos sin ] where is the shortwave radiation (MJ m −2 ·day −1 ); is the canopy reflectance of vegetation with a value of 0.23; is the sunshine duration (h); is the maximum possible sunshine duration (h); is the solar radiation at the top of the atmosphere (MJ m −2 ·day −1 ); is the solar constant with a value of 0.082 (MJ m −2 min −1 ); is the inverse relative distance of Earth-Sun; is the angle of the sunset sun (rad); is the latitude (rad); is the solar declination (rad); is the net long wave radiation (MJ m −2 ·day −1 ); is the Stefan-Boltzmann constant, the value is 4.903 × 10 −9 (MJ K −4 m −2 day −1 ); , is the maximum daily air temperature (K); , is the minimum daily air temperature (K); is the value of short wave radiation on the surface of vegetation on sunny days (MJ m −2 ·day −1 ). Ye, J.S. et al. [65] found that = 0.24, = 0.6 and = (0.64 + 5.48 * 10 ) are suitable for the QTP, where Z is the altitude (m).

Data Preprocessing
Before trend analysis for ET0, the Ljung-Box test (LB test) was used to test the autocorrelation of the data, and the Cochrane-Orcutt method was used to transform the variables of the sites with autocorrelation.
The LB test is widely used in white-noise monitoring of time series, especially in the field of meteorology [66,67]. First, assuming that the sequences are completely uncorrelated, this method constructs a statistic Q: where n is the number of the sample, is the autocorrelation coefficient of the k-order lag of the sample, and the statistic Q obeys the chi-square distribution. At the α significance level, the reject domain is Q > , . When the p-value of Q (h) is less than α, the original hypothesis is rejected and the sequence receives autocorrelation.
At a 0.05 significance level, this study calculated the p-value at the first order lay of the original sequence by the LB test (Table S1 in supplementary information). If the sequences with autocorrelation feature, the Cochrane-Orcutt method [68] was used for variable transformation: Execute the LB test again on the new sequence after conversion. If autocorrelation still exists, execute the variable conversion again by Cochrane-Orcutt until autocorrelation is removed. Table S1 (in supplementary information) shows the p-value of the original sequence and the converted sequence. It can be seen that many sites have autocorrelation, but the autocorrelation can be removed after one conversion.

Cumulative Anomaly Curve
A cumulative anomaly curve can intuitively indicate the long-term evolution trend and continuous change of a sequence, and it is a common method to judge the change trend of hydrological and meteorological sequences [69,70]. The calculation procedure of the cumulative anomaly is as follows: where is the meteorological dataset, n is the length of the dataset. The anomaly is the distance to the average value, which reflects the degree of data dispersion. The rising of the cumulative anomaly curve showed an upward trend, while the opposite showed a downward trend. Since the mutation point must appear near the peak or inflection point of the cumulative anomaly curve, the mutation year can be roughly tested.

Mann-Kendall Test
Recommended by the WMO, the Mann-Kendall (MK) test is widely used as a nonparametric test of temperature, precipitation and other factors [71]. The advantage is that the samples do not need to follow a certain distribution and are not disturbed by small fluctuations. The MK test can as well reveal the trend and mutation of a dataset [72].
Define two statistics named S and Z, the calculation process is as follows: where and are two sequential values of the meteorological dataset (1 < i < j < n), n is the length of the dataset; ( ) is a symbolic function; Z is an indicator for the severity of the change trend. The positive value indicates that the trend is rising, while a negative value indicates that the trend is falling. If | | ≥ ( )/ , it means that at an significance level, the data series change significantly.
Define as the sum of the cumulative numbers when > (1 < i < j < n). The mathematical expectation and variance of are: where 2 ≤ k ≤ n, = 0; at an α significance level, if ≥ ⁄ , the data series change significantly.
is the negative sequence of . When the intersection point of the curve and curve is within the confidence interval, the point is a mutation point.

Moving T-test
The moving t-test examines whether the difference between the average values of the two datasets are significant to detect mutations, which is widely used in the study of time series [73]. The calculation method of the statistic t is as follows: where and are two parts of the datasets, and are their lengths, and are their average values, and and are their standard deviations.
The statistic t is submitted to the t-distribution, for which freedom v = + − 2. At significance level, if t ≥ | |, the point is a mutation point.

Yamamoto Test
The Yamamoto test defines the signal-to-noise ratio (SNR) of meteorological factors, which is simple, effective and intelligible [74]. Its calculation formula is as follows: Similar to the moving t-test, where and are two parts of the datasets, and are their lengths, and are their average values, and and are their standard deviations.
If SNR > 1.0, there is a mutation in the point. Furthermore, if SNR > 2.0, it can be considered a strong mutation year.
Since the current methods for identifying mutation points have their own advantages and disadvantages, this study believes that when the results of more than two methods are consistent, specific mutation points are determined.

Path Analysis
Path analysis is an extension of the multiple linear regression analysis, which can identify and quantify the direct and indirect effects of independent variables on dependent variables [75]. It is a powerful tool used to analyze the relationship between multiple variables and is widely used in many fields. The analysis process is as follows: First establish a multiple regression equation: Then, standardize the variables of Equation (19): where , , , … and are different meteorological factors, is ET0, is the estimated value obtained by the least squares regression, is the standard deviation, is the estimated value of the independent variable , is the standard deviation of , and is a partial regression coefficient. The direct path coefficient is calculated according to the following formula: Establish a path analysis model: where is the correlation coefficient between and , is the correlation coefficient between and indicates the direct influence degree of an independent variable on a dependent variable.
is the indirect path coefficient, which indicates the indirect influence degree of on through .
indicates the influence of factors not considered by the dependent variables. Figure 3 shows the relationship between variables in path analysis. In this study, path analysis was used to identify and quantify the direct and indirect effects of the nine climatic elements, including the daily actual atmospheric pressure (AP), mean temperature (Tmean), maximum temperature (Tmax), minimum temperature (Tmin), precipitation (P), relative humidity (RH), sunshine duration (SD), wind speed (U2) and net radiation ( ) on ET0.
To realize these methods, we used R programming and the SPSS software. Figure 4 shows the data processing and research procedure of this study.  ... ...

Mutation Analysis of Annual ET0
The cumulative anomaly curve, Mann-Kendall (MK) test (α = 0.05), moving t-test (α = 0.05) and Yamamoto test were used to distinguish the mutation year of the ET0 series annually and for four seasons from 1961 to 2017. Figure 5 shows that there was a turning point from increase to decrease in 1988; nevertheless, a turning point from decrease to increase in the QTP in 1968 and 2005 can be found according to the cumulative anomaly curve. Furthermore, in 1988, the moving t-test curves of the three steps all exceeded the upper critical line, and their SNR reached a peak level. Moreover, the moving t-test and Yamamoto test curves of 5-year and 10-year sliding steps exceeded the critical line in 2005. In addition, the intersection points of the MK curves in the confidence interval were in 1964 and 2012, since the sequences before and after them were relatively short, neither of them can be determined as a true mutation point. In summary, it can be considered that a descending mutation occurred in 1988, while an ascending mutation occurred in 2005 in the QTP. Figure 5 illustrates that the results of the moving t-test and Yamamoto test in the medium and high-elevation zones were more inclined to 1988 and 2005. Although the positive peak value of the cumulative anomaly curve in the medium-elevation zone was 1981, 1988 was also its inflection point and the decline thereafter became more pronounced. The year 2005 was the negative peak value of the cumulative anomaly curve. Therefore, it can be determined that the descending mutation occurred in 1988, and the ascending mutation occurred in 2005. However, in the low-elevation zone, 1981 was the positive peak of the cumulative anomaly curve and the only intersection of the MK curve in the confidence interval. Both the moving t-test and the Yamamoto test exceeded the critical line at the point. Therefore, it can be concluded that the descending mutation occurred in 1981. In conclusion, the ET0 series in the QTP had two mutations during 1961-2017, with a descending mutation occurring in the 1980s and an ascending mutation occurring in 2005. Zhao et al. [76] and Li et al. [77] also found that mutations of climatic factors occurred in China in the 1980s and the early 21st century. Furthermore, with the increase in altitude below 5000 m, the period of the descending mutation of ET0 was prolonged, and the amplitude of ascending was enhanced.

Mutation Analysis of Four Seasons of ET0
Figure S1 (in supplementary information) shows that the results of the QTP are consistent with the low and the medium-elevation zones in the spring, the descending mutations occurred approximately in 1981, and the ascending mutations occurred in 2005. However, in the high-elevation zone, the peak value of the cumulative anomaly curve was not prominent, the MK curve had more intersections in the confidence interval, and the sliding t-value and SNR value of the three sliding steps did not exceed the critical line, meaning no mutation was detected in the high-elevation zone. Table S2 and Figure S2 (in supplementary information) show that in the summer, the descending mutations in the QTP and the medium-elevation zone occurred in 1988-1991 because 1991 was the intersection point of the MK curve in the confidence interval, and 1988-1991 was the peak segment of their cumulative anomaly curve and passed the moving t-test for the 15-year sliding step. However, in the low-elevation zone, the intersection point of the MK curve in the confidence interval was 1985 and 1987. The other three methods were more inclined to 1981, which can roughly suggest that the descending mutation occurred in 1981-1987. In addition, in the high-elevation zone, the ascending mutations in 2005 were more significant than those in the lower ones. Figure S3 (in supplementary information) shows that in autumn, the QTP and the medium-elevation zone had a descending mutation in 1987 and that for the low-elevation zone was in 1994; furthermore, the medium-elevation zone in 2005 showed a significant ascending mutation. Figure S4 (in supplementary information) shows that the descending mutation of the winter ET0 in the QTP, low and medium-elevation zones failed to pass the test since no more than two methods had the same results, and there was an ascending mutation in the high-elevation zone in 2004.
In conclusion, with the increasing altitude below 5000 m, the delay phenomenon of the descending mutation of ET0 was verified in the summer, and the enhancement effect on the ascending mutation was verified in the summer and winter series.

Temporal Trend of Annual ET0
The trend ratios and Z-values from the MK test of annual ET0 in different periods are listed in Table 1 and Table S3. From 1961 to 2017, ET0 on the QTP showed an increasing trend, which is consistent with the findings of Wang et al. [78], but different from the decreasing trend obtained by Thomas [79] and Kuang and Jiao [80] because the period of the ET0 series was selected differently. Table 1 also shows that as the altitude increases, the increasing trend of the ET0 sequence becomes more significant. Before the first mutation, ET0 in different elevation zones of the QTP showed an increasing trend, and the increasing trend in the low-elevation zone was the most significant by 28.56 mm per decade before 1981. Between the two mutations, there are negative Z-value and change ratios of ET0 in the QTP, low and high -elevation zones, indicating that their first reduction mutation was a transit mutation, and the middle-elevation zone was a rate or mean value mutation. After the second mutation occurred in 2005, the ET0 of each elevation zone of the QTP showed an increasing trend, and with the increase in altitude, this trend became more significant, which was consistent with 1961-2017. Thomas [79] also observed the positive relation between evapotranspiration change and station altitude in the mountains of Southwest China. Trend ratios is the slope of linear regression between meteorological factors and year; the slope unit of ET0 is mm/decade. The three periods before and after the mutation for the QTP, the middle and high-  Figure 6 shows the distribution of the Z-value of ET0 in the four seasons from the MK test from 1961 to 2017. The graph shows that the change trend of ET0 varied greatly in different seasons. In winter, the change trend of ET0 in the whole region was the most consistent, with 89.13% of the stations showing an increasing trend, and 82.93% of the stations showing a significant increasing trend (p < 0.05). Kang et al. [81] and Zhang et al. [82] also concluded that the warming trend in the winter has become more significant. In the spring, the small Z-value bubbles of ET0 on the stations indicate that the change trend was the gentlest, and the spring ET0 shows a decreasing trend at the low-elevation stations of 63.64%, while most of the middle and high-elevation zones show an increasing trend, with the proportion of stations being 68.18% and 84.62%. The low elevation zone of the QTP is mainly located in the Qaidam Basin, where the ET0 decreased significantly in summer and autumn, which is obviously different from the significant increase in other elevation zones.  25, the ET0 change is extremely significant. The red bubble represents the increasing trend, the blue represents the decreasing trend, and the larger the bubble is, the more significant the change trend is.

Climate-Driven Pattern of ET0 in Different Altitudes
According to the mutation year determined in Section 3.1, the annual ET0 dataset is divided into three series, which are before the first mutation, between two mutations and after the second mutation. Path analysis was used to quantify the relationship between ET0 and climate factors. Figure 7 shows that before the first mutation, ET0 in the QTP was mainly related to the factors, including Tmean, U2, SD, AP and Tmax, of which U2 and Tmax had a strong promotion effect on the change of ET0, and AP directly inhibited ET0 and indirectly inhibited ET0 by suppressing wind speed. Moreover, after the descending mutation occurred in 1988, the main climatic factors affecting ET0 were Tmax, U2 and Rn. The direct promotion effect of Tmax and U2 was greatly enhanced. Furthermore, after the ascending mutation in 2005, the main meteorological factor affecting the change of ET0 was only Tmax, and the combined effect of other factors accounts for only 31.69%.
In the low-altitude zone, the main influencing factors were Tmean, U2, SD and RH before the first mutation, among which the direct promotion of U2 was the largest, and RH directly inhibited ET0. Liu et al. [41] and Fan and Thomas [39] came to a similar conclusion in China. After the descending mutation occurred in 1981, the main impact factors were RH, Rn, U2 and Tmax. After 2005, the main impact factors were Tmax, with a contribution rate of 63.16%. In the medium-elevation zone, Tmean, U2 and SD directly promoted the extent of the change of ET0, and Tmin and AP played an inhibitory role before the first mutation. From 1987 to 2004, the direct promotion of Tmax and U2 was greatly enhanced. After 2005, Tmax and Rn were the main climatic factors, and both contributed. In the high elevation zone, the main climatic factors and driving patterns were consistent with the low-elevation zones before the first mutation. From 1988 to 2004, the direct promotion of Tmean and U2 was enhanced, and Tmin had a significant inhibitory effect. After 2005, the main impact factors became Tmax, RH, U2 and SD.
In general, before the first mutation, Tmean and U2 played a major role in the change of ET0 in QTP and the three elevation zones. Between the two mutations, the promotion effect of Tmean and U2 in the QTP, low and medium-elevation zones was significantly enhanced, and Rn became the main influencing factor in the QTP, low and high-elevation zones. After the ascending mutation occurred in 2005, the direct promotion of Tmax in QTP and the three elevation zones contributed the most to the change of ET0. Zhang et al. [83] also found that Tmax had the greatest impact on China's ET0.  25, the ET0 change is extremely significant. The red bubble represents the increasing trend, the blue represents the decreasing trend, and the larger the bubble is, the more significant the change trend is.

Conclusions
Based on the daily meteorological data of 46 meteorological stations in the QTP from 1961 to 2017, the daily ET0 was calculated by the Penman-Monteith formula. Before trend analysis, we used the Ljung-Box test to verify the autocorrelation of data and used the Cochrane-Orcutt method to transform the variables of the sites with autocorrelation. Furthermore, the cumulative anomaly curve, Mann-Kendall test, moving t-test and Yamamoto test were implemented to identify the mutation year of the ET0 series in different elevation zones. The combination of the four methods can avoid the error caused by the accidental fluctuation of the sequence and make the recognition results more accurate. In this study, the mutation of the ET0 series, temporal trend of different elevation zones, and the driving mode of meteorological factors before and after the mutation were explored. The main conclusions are as follows: (1) The annual and four season ET0 series in the QTP experienced a descending mutation in the 1980s and an ascending mutation in 2005. The descending mutation of the low-elevation zone (below 2800 m) was in 1981. As the altitude (below 5000 m) increased, the descending mutation of the annual ET0 tended to be delayed, while the ascending mutation was more significant. The descending mutation of ET0 in the spring appeared in 1981. The mutation features of the QTP and three elevation zones in the summer were similar to those of the annual scale, but not obvious in winter.
(2) The annual ET0 series of the Qinghai Tibet Plateau and different elevation zones showed an increasing trend from 1961 to 2017, and it increased more significantly with the increase in elevation below 5000 m. The increasing trend was the most consistent and intense in the winter, with 89.13% of stations showing an increasing trend. Regarding the impact of mutations on trend changes, ET0 was more prone to mutation of the rate or mean value in the middle-elevation zone, and transit mutation, in which the change trend before and after the mutation is opposite, was more likely to occur in other elevation zones.
(3) Before and after the two mutations, the climate-driven patterns of different elevation zones were quite different, but on the whole, before the descending mutation, the Tmean and U2 played a major role in promoting the change of ET0 in the QTP and the three elevation zones. Between the two mutations, the promotion effect of Tmean and U2 in the QTP, low and medium-elevation zones was significantly enhanced, and Rn became the main influencing factor in the QTP, low and high-elevation zones. After the ascending mutation occurred in 2005, Tmax became the most important driving factor for the QTP and the three elevation zones.
This study explored the different trends of evapotranspiration at different altitudes, which adds new understanding to the regional hydrology research. The results of the study find that there is an increasing trend of reference evapotranspiration in some areas, where the department should conduct drought warning and water resources management to avoid a water shortage. In the future, regional agricultural irrigation plans can be formulated in combination with snow melting or river discharge research.

Supplementary Materials:
The following are available online at www.mdpi.com/article/10.3390/w13131749/s1, Table S1: The p-value of the raw sequence and the converted sequence by LB test, Figure S1: The results of mutation analysis of spring ET0 on different elevation zones, Figure S2: The results of mutation analysis of summer ET0 on different elevation zones, Figure S3: The results of mutation analysis of autumn ET0 on different elevation zones, Figure S4: The results of mutation analysis of winter ET0 on different elevation zones, Table S2. The mutation years for ET0 of four seasons as a result of section 3.1.2, Table S3. Z-values from the MK test of ET0 in four seasons in different elevation zones in different periods.
Author Contributions: Conceptualization, methodology and software, Y.L.; validation, X.Y. and Q.W.; formal analysis and writing-original draft preparation, Y.L. and Q.W.; investigation, L.L.; resources and data curation, W.J. and Q.J.; supervision, project administration and funding acquisition, J.Y. All authors participated in the review and editing. All authors have read and agreed to the published version of the manuscript. Institutional Review Board Statement: Not applicable.

Informed Consent Statement: Not applicable.
Data Availability Statement: Raw datasets are available online at http://data.cma.cn/. Readers can contact the corresponding author for processed data of this study.