Characteristics and Drivers of Reference Evapotranspiration in Hilly Regions in Southern China

: This paper has adopted related meteorological data collected by 69 meteorological stations between 1951 and 2013 to analyze changes and drivers of reference evapotranspiration (ET 0 ) in the hilly regions located in southern China. Results show that: (1) ET 0 in southern China’s hilly regions reaches its maximum in summer and its minimum in winter, and that the annual ET 0 shows an increasing trend. ET 0 happened abrupt change due to the impact of abrupt meteorological variables changes, and the signiﬁcant year of mutation were 1953, 1964 and 2008. Most abrupt changes of ET 0 in meteorological stations occurred in the 1950s and 1960s. (2) The low value of ET 0 was mainly captured in high-altitude areas. Spatially, the ET 0 in the east was higher than that in the west. With the exception of a handful of stations, the trend coe ﬃ cients of ET 0 were all positive, exhibiting a gradual rise. Changes in ET 0 in the east were much more sensitive than that in the west. Since ET 0 was a ﬀ ected by the cyclical changes in relative humidity, short-period oscillations were observed in all these changes. (3) In general, the ET0 was negatively correlated with relative humidity, and positively correlated with temperature and sunshine percentage. ET 0 is most sensitive to changes in average temperature, with a sensitivity coe ﬃ cient of 1.136. ET 0 showed positive sensitivity to average temperature and sunshine hours, which were notable in the northeastern, and uniform in the spatial. ET 0 showed negatively sensitivity to relative humidity, and the absolute value of sensitivity coe ﬃ cient in the northwestern is smaller. The highest contribution to ET 0 is the average temperature (6.873%), and the total contribution of the four meteorological variables to the change of ET 0 is 7.842%. The contribution of average temperature, relative humidity, and sunshine hours to ET 0 is higher in the northern and eastern, northern, northern and eastern areas, respectively. Climate indexes (Western Paciﬁc Index (WP), Southern Oscillation Index (SOI), Tropical Northern Atlantic Index (TNA), and El Niño-Southern Oscillation (ENSO)) were correlated with the ET 0 . In addition, the ET 0 and altitude, as well as the latitude and longitude were also correlated with each other.


Introduction
As the population grows, the threat of climate change to global food security has become a major challenge in the 21st century [1]. Evapotranspiration (ET) is an important means of surface water consumption, and it is also an important variable of the hydrological cycle in the atmosphere terrestrial system [2]. The changes of evapotranspiration influence various regional water equilibrium components, including the production and life in the region, and agriculture production especially. Reference evapotranspiration (ET 0 ) is a reflection of the energy available to evaporate water under adequate water supply conditions. Meteorological factors, such as wind speed, temperature, solar radiation and precipitation have various changes under the conditions of global warming and climate change in China, and these factors are closely related to changes in reference evapotranspiration [3]. With wind available to transport the water vapor from the ground up into the lower atmosphere, evapotranspiration accounts for the movement of water within a plant, and the subsequent loss of water as vapor through stomata in its leaves. Reference evapotranspiration is said to equal actual evapotranspiration when there is ample water, and it is an important part of the water cycle [4]. Although the actual evapotranspiration is an objective variable for measuring water changes, the lack of long-term and large-scale actual evapotranspiration data makes the study of reference evapotranspiration important for understanding the actual evapotranspiration [5]. As an important indicator to measure the regional water status, reference evapotranspiration not only has an impact on the sustainable development of the social economy, but also has a profound impact on regional agricultural production, especially in terms of food security [6]. Therefore, a scientific basis can be provided for the Equationtion of China's agricultural change by analyzing the dynamic and spatial characteristics of reference evapotranspiration in the assessment area, and exploring its response to global climate change.
Quite a few studies have been conducted on reference evapotranspiration in the world, (especially in selected ecologically fragile areas). The research on reference evapotranspiration is mainly in terms of two aspects. One is the temporal and spatial variation of reference evapotranspiration, and the other is the influence of meteorological elements on the reference evapotranspiration. Some scholars focus on the analysis of spatiotemporal variation of reference evapotranspiration. Wang et al. [7] researched the variety characters reference evapotranspiration in the Zoige wetland, and Xu et al. [8] analyze the spatial distribution and temporal trends in the Changjiang catchment of China. There are still many scholars who are concerned with the analysis of the influencing factors of reference evapotranspiration. While Yang et al. [9] adopted the Thornthwaite and Penman-Monteith methods to conduct sensitivity analysis on global reference evapotranspiration, Roderick et al. [10] concluded that evaporation in the Northern Hemisphere underwent an overall reduction over the past 50 years and discussed the possible causes. De la Casa et al. [11] studied the change of reference evapotranspiration in central Argentina. Results shows that beneficial change in agricultural suitability reported for Central Argentina region was almost exclusively produced by increasing plus precipitation. Croitoru et al. [12] used the monthly meteorological data of 57 sites in Romania to analyze the variation characteristics of reference evapotranspiration under the change of meteorological variables in Romania from 1961 to 2007. In addition, some scholars have also studied the reference evapotranspiration in the Northeastern [13,14], Northwestern [15] and Southwestern areas [16] as well as, the Changjiang Basin [8,17], of China. Huang Huiping [18][19][20] and other scholars have even pointed out that China's annual average reference evapotranspiration underwent an overall reduction over the past 50 years. In short, due to the harsh climatic conditions and fragile ecological environment in northern China [21], especially the northwestern region, numerous studies on surface water have been successfully conducted in the northwestern region. However, the research is often overlooked in terms on agricultural climate resources and water balance particularly in the Jiangnan hilly areas in the humid climate zone due to the so-called "excellent climatic conditions and good ecological environment." The precipitation is abundant in the Jiangnan hilly area, as the most important grain producing area in China. However, it has seen an upward trend of the reference evapotranspiration and extreme hydrological events, such as droughts and floods in higher frequency and on a larger scale in recent decades [22,23]. The purpose of the paper is to analyze the spatial and temporal trends of reference evapotranspiration and the driving factors of reference evapotranspiration in the Jiangnan hilly area for more than 60 years. In this regard, hence the study on the reference evapotranspiration in the hilly areas of the south of the Yangtze River can provide a scientific basis for the Equationtion of regional agricultural planning and agricultural development strategies and agricultural response to climate change.

Overview of Research Areas
Hilly regions in southern China usually refer to the areas in south of the Yangtze River, north of Nanling, east of Xuefeng Mountain, and west of Mount Wuyi (Figure 1), which is approximately about 370,000 km 2 in area. As these areas fall within the typical subtropical monsoon climate zone, they are characterized by hot and rainy summers, as well as mild and dry winters. The annual precipitation is about 1300-1800 mm while annual average temperature is about 16-20 • C. Research areas feature abundant water sources, numerous rivers and lakes, mainly including the Dongting Lake Basin and Poyang Lake Basin. Dongting Lake Plain and Poyang Lake Plain, with a long history in agriculture production, are also significant grain producing areas in southern China.
Water 2019, 11, 1914 3 of 22 and the driving factors of reference evapotranspiration in the Jiangnan hilly area for more than 60 years. In this regard, hence the study on the reference evapotranspiration in the hilly areas of the south of the Yangtze River can provide a scientific basis for the Equationtion of regional agricultural planning and agricultural development strategies and agricultural response to climate change.

Overview of Research Areas
Hilly regions in southern China usually refer to the areas in south of the Yangtze River, north of Nanling, east of Xuefeng Mountain, and west of Mount Wuyi (Figure 1), which is approximately about 370,000 km 2 in area. As these areas fall within the typical subtropical monsoon climate zone, they are characterized by hot and rainy summers, as well as mild and dry winters. The annual precipitation is about 1300-1800 mm while annual average temperature is about 16-20 °C . Research areas feature abundant water sources, numerous rivers and lakes, mainly including the Dongting Lake Basin and Poyang Lake Basin. Dongting Lake Plain and Poyang Lake Plain, with a long history in agriculture production, are also significant grain producing areas in southern China.

Data Source
This paper has selected various meteorological data collected between 1951 and 2013 by 69 meteorological stations located in the hilly regions in southern China. The data includes different

Data Source
This paper has selected various meteorological data collected between 1951 and 2013 by 69 meteorological stations located in the hilly regions in southern China. The data includes different meteorological variables such as temperature, wind speed, relative humidity and sunshine hours, etc. Meteorological data had been provided by the National Meteorological Science Data Sharing Network.

Reference Evapotranspiration
Several methods/Equations, such as Penman-Monteith [24], Hargreaves and Samani [25], Thornthwaite [26] and Makkink [27], can be applied to calculate reference evapotranspiration. At present, the Penman-Monteith method, recommended by FAO (Food and Agriculture Organization of the United Nations) in 1998, is the most widely adopted one. This paper has also applied the Penman-Monteith method to calculate reference evapotranspiration, hence taking various meteorological elements, including temperature, wind speed, relative humidity and sunshine hours, etc., into account.
In Equation (1), ET 0 represents reference evapotranspiration (mm/day), R n the net surface radiation amount (MJm −2 day −1 ), G the soil heat flux (MJm −2 day −1 ), T the average temperature ( • C), u 2 the wind speed (ms −1 ) at the height of 2 m, e s the saturated vapor pressure (KPa), e a the actual vapor pressure (KPa), e s − e a the saturated vapor pressure difference (KPa), ∆ the slope of the vapor pressure curve (KPa) • C −1 ), and γ the constant of the hygrometer.
There are many methods for calculating net radiation. The calculation method of net radiation used in this paper is the method of net short-wave radiation minus net long-wave radiation: where R ns is short-wave radiation and R nl represent long-wave radiation. Also: where a is surface reflectivity, and the value is 0.23; R a means solar radiation at the top of the atmosphere (MJm −2 day −1 ); n is sunshine hours (h); N is maximum sunshine hours (h); a s represents the amount of solar radiation reaching the surface of the Earth in cloudy days (n = 0); b s imply the amount of solar radiation reaching the surface of the earth in a sunny day (n = N); the recommended values of a s is 0.25 and b s is 0.50 in this paper.
where G SC is solar constant (0.0820 MJm −2 min −1 ), d r represents reciprocal of the relative distance between the sun and the earth, w s is the solar time angle (rad), ϕ is the geographic latitude and δ is the solar magnetic declination (rad). R nl is long-wave radiation, the long-wave radiation is affected by various variables such as water vapor and temperature. According to the Stefan-Boltzmann law, the long-wave radiation calculation Equation used in this paper is: where σ is the Stefan-Boltzmann constant (4.90 × 10 −9 MJK −4 m −2 day −1 ), T 4 max,k is the maximum absolute temperature; T 4 min,k is the lowest absolute temperature (K = • C + 273.6), e a is the actual vapor pressure (KPa).
where G is soil heat flux (MJ m −2 day −1 ), C S is the soil heat capacity (2.1 MJ m −3 • C −1 ), ∆Z is effective soil depth (0.2 m), T I means the average air temperature of the ith day, and T i−1 is the average temperature of the (i−1)th day.

Accumulative Anomaly
Cumulative anomaly is a method used to visually judge the trend of climate change by curve, and its accumulative anomaly is an index to distinguish the changing tendency of discrete data. It is the accumulation of the difference between a value and the average of a series of values. The accumulative anomaly can be calculated by the Equation (7); where LP i is the anomaly cumulative value of the i-th year; R i is the ET 0 of the i-th year; R is the mean value of the series R i , and N is the number of discrete points. Concerning the positive and negative characteristics of the anomaly, when the LPi continues to increase, it indicates that the anomaly continues to be positive during the period, that is, the ET 0 is higher than the average; When the LPi decreases to zero, it indicates that the anomaly of the period is constant, that is, the ET 0 is kept at an average level. When the LPi continues to decrease, it indicates that the anomaly continues to be negative during the period, that is, the ET 0 is lower than the average level. Therefore, the interannual variation phase of ET 0 can be determined more intuitively and accurately [28].

Climatic Trend Coefficients
The study of climatic change trends is based on climatic trend coefficients, hence changes in these coefficients will indicate any increment/drop in various meteorological elements [14]. The value interval is [−1,1]. The closer r xt is to −1, the more significant the drop is; the closer r xt is to 1, the more significant the increment is.
In Equation (8), x i represents the element's value in year i, and x is the sample average; when t = (n + 1)/2, and r xt is positive (negative), it indicates that the element tends to undergo a linear increment (drop) within n years.

Sensitivity Coefficient and Contribution
The calculation method of the sensitivity coefficient in this paper is the method proposed by McCucn [29]. It is used to calculate the sensitivity coefficient of reference evapotranspiration to key climatic variables, by using the partial derivative of ET 0 for each meteorological factor.
where S Vi is the sensitivity coefficient; ∆ET 0 is the variations of ET 0 . V i , ∆V i is meteorological variables and their variations respectively. If S Vi > 0, it means that the reference evapotranspiration increased as the variable increases. If S Vi < 0, it indicates that the reference evapotranspiration decreased as the variable increases. The larger S Vi , the greater the effect a given variable has on reference evapotranspiration [17]. Based on literature review, and as used by other workers [30], meteorological parameters (temperature, wind speed, sunshine hours and relative humidity) have been selected for a range of specific variability (±20%) for the analysis.
Multiplying the sensitivity coefficient of a single meteorological factor by the multi-year relative change of the variable can obtain the change of reference evapotranspiration caused by this variable, that is, the contribution of the element to the change of reference evapotranspiration [31]. The specific expression is as follows: where, Con Vi is the contribution of meteorological factor Vi to ET 0 change, S Vi is the sensitivity coefficient of V i , RC Vi (%) is the relative change of V i for many years, a v is the multi-year average of V i , and Trend is the year-to-year rate of change of V i , Trend is calculated by the trend analysis method [32]. Reference evapotranspiration is affected by multiple meteorological variables. In this paper, the contribution of the four variables of average temperature, relative humidity, sunshine hours and wind speed is added to obtain the total contribution to the change of reference evapotranspiration [33]. The Equation is: Con ET 0 = Con T + Con RH + Con n + Con U (12) where, Con ET0 represents the change of ET 0 caused by the interaction of four meteorological variables, also called the estimated change of ET 0 . Con T , Con RH , Con n , and Con U represent the contribution of average temperature, relative humidity, sunshine hours, and wind speed to changes in ET 0 respectively.

Wavelet Analysis
Wavelet analysis (also known as multi-resolution analysis) is a commonly used method for analyzing the scale and trend of time series, studying the evolution of different scales (cycles) over time, with multi-resolution analysis and adaptive characteristics to the signal because the periodic transformation characteristics of meteorological elements are complex. It contains periodic changes of various time scales in the same period, and they will exhibit multiple time scale features. Therefore, the use of wavelet analysis to study the evolution of meteorological elements at different scales (cycles) over time has become an important method for studying the long-term changes of meteorological elements [34].   Figure 4 shows the seasonal changes and annual average change in reference evapotranspiration in the hilly regions in southern China over the last 63 years. A large gap between the reference evapotranspiration in summer and winter (the peak in summer and the lowest point in winter) can be observed due to the large differences in sunlight and temperature in summer and winter. From the linear change trend, reference evapotranspiration increased with time throughout all seasons. However, the most significant change occurred in spring with the highest correlation coefficient (Figure 4a).

Seasonal Changes and Annual Average Change in Reference Evapotranspiration
Between 1951 and 2013, reference evapotranspiration showed a fluctuating and rising trend, passing the significance test of 0.01 ( Figure 4b). The average reference evapotranspiration for several years was 1284 mm, with its highest value in 2013 (indicating a surface evapotranspiration peak) and lowest value in 1952 (indicating the lowest point in surface evapotranspiration). Although average reference evapotranspiration had seen increments and reductions over 5 years, it revealed an overall fluctuating and rising trend, and it increased significantly in 2001, mainly due to the increase in solar radiation after 2001, resulting in significant fluctuations in temperature and evapotranspiration [35].  Figure 4 shows the seasonal changes and annual average change in reference evapotranspiration in the hilly regions in southern China over the last 63 years. A large gap between the reference evapotranspiration in summer and winter (the peak in summer and the lowest point in winter) can be observed due to the large differences in sunlight and temperature in summer and winter. From the linear change trend, reference evapotranspiration increased with time throughout all seasons. However, the most significant change occurred in spring with the highest correlation coefficient (Figure 4a).

Seasonal Changes and Annual Average Change in Reference Evapotranspiration
Between 1951 and 2013, reference evapotranspiration showed a fluctuating and rising trend, passing the significance test of 0.01 ( Figure 4b). The average reference evapotranspiration for several years was 1284 mm, with its highest value in 2013 (indicating a surface evapotranspiration peak) and lowest value in 1952 (indicating the lowest point in surface evapotranspiration). Although average reference evapotranspiration had seen increments and reductions over 5 years, it revealed an overall fluctuating and rising trend, and it increased significantly in 2001, mainly due to the increase in solar radiation after 2001, resulting in significant fluctuations in temperature and evapotranspiration [35]. Furthermore, changes in cumulative anomalies also reflect the changing trend of the research object. The reference evapotranspiration cumulative anomaly refers to the accumulation of the difference between reference evapotranspiration values and the multi-year average reference evapotranspiration during the calculation period. Between 1951 and 2013, the reference evapotranspiration cumulative anomaly curve went through three main stages: gradual fluctuating rise (before 1966), gradual fluctuating decline (1966 to 2002), rapid rise (after 2002). According to these three stages, ET 0 was segmentally fitted, and partition fitting also indicated the changing trend of reference evapotranspiration featuring a rise-slow decline-sharp rise pattern.
Water 2019, 11, 1914 9 of 22 Furthermore, changes in cumulative anomalies also reflect the changing trend of the research object. The reference evapotranspiration cumulative anomaly refers to the accumulation of the difference between reference evapotranspiration values and the multi-year average reference evapotranspiration during the calculation period. Between 1951 and 2013, the reference evapotranspiration cumulative anomaly curve went through three main stages: gradual fluctuating rise (before 1966), gradual fluctuating decline (1966 to 2002), rapid rise (after 2002). According to these three stages, ET0 was segmentally fitted, and partition fitting also indicated the changing trend of reference evapotranspiration featuring a rise-slow decline-sharp rise pattern.

Spatial Changes in Reference Evapotranspiration and Trend Coefficient
Hilly regions located to the south of the Yangtze River are large and wide, featuring complex and diverse terrains as well as densely-distributed waters. As different landscape types result in differences in the spatial distribution of climatic variables, these differences will always influence the spatial distribution of reference evapotranspiration. As shown in Figure 5, the reference evapotranspiration in northwest and southeast in spring was relatively higher, and the reference evapotranspiration in Southwest and central east was lower, with the trend coefficient between 0.31 and 0.6. Stations with a positive trend coefficient accounted for 89.9%. In most areas, reference evapotranspiration showed a rising trend. Only a handful of areas in the east had revealed a downward trend. In summer, reference evapotranspiration decreased from the southwestern to the northeastern area, and the northwest is also lower; most trend coefficients fell between −0.29 and 0.0, and the stations with a positive trend coefficient accounted for 47.8%. The scope of the rising reference evapotranspiration trend was comparable to that of its downward trend, with the rising trend mainly occurring in the northeastern and southwestern areas. In autumn, reference evapotranspiration was generally higher in the east than that in the west, with the trend coefficient between 0.01 and 0.4. Stations with positive trend coefficient accounted for 69.6%, indicating an overall rising trend of the regional reference evapotranspiration in autumn. However, reference evapotranspiration in the southeast and northeast was relatively higher in winter, which was similar to the spatial distribution pattern of reference evapotranspiration in spring. Meanwhile, most trend coefficients fell between 0.01 and 0.3, and stations with a positive trend coefficient accounted for 85.5%. Although reference evapotranspiration in most areas revealed a rising trend, stations showing a downward trend were mainly distributed in the northwestern areas.
Areas with low annual average reference evapotranspiration coincided with that of the seasonal reference evapotranspiration, but large differences remained in high-value areas. Annual average reference evapotranspiration was between 1134 mm and 1387 mm, and reference evapotranspiration in the eastern region was higher than that in the western region, which coincided with the spatial distribution of reference evapotranspiration in autumn. Areas with lower annual average reference evapotranspiration were mainly distributed in the western mountains with higher altitudes, mainly

Spatial Changes in Reference Evapotranspiration and Trend Coefficient
Hilly regions located to the south of the Yangtze River are large and wide, featuring complex and diverse terrains as well as densely-distributed waters. As different landscape types result in differences in the spatial distribution of climatic variables, these differences will always influence the spatial distribution of reference evapotranspiration. As shown in Figure 5, the reference evapotranspiration in northwest and southeast in spring was relatively higher, and the reference evapotranspiration in Southwest and central east was lower, with the trend coefficient between 0.31 and 0.6. Stations with a positive trend coefficient accounted for 89.9%. In most areas, reference evapotranspiration showed a rising trend. Only a handful of areas in the east had revealed a downward trend. In summer, reference evapotranspiration decreased from the southwestern to the northeastern area, and the northwest is also lower; most trend coefficients fell between −0.29 and 0.0, and the stations with a positive trend coefficient accounted for 47.8%. The scope of the rising reference evapotranspiration trend was comparable to that of its downward trend, with the rising trend mainly occurring in the northeastern and southwestern areas. In autumn, reference evapotranspiration was generally higher in the east than that in the west, with the trend coefficient between 0.01 and 0.4. Stations with positive trend coefficient accounted for 69.6%, indicating an overall rising trend of the regional reference evapotranspiration in autumn. However, reference evapotranspiration in the southeast and northeast was relatively higher in winter, which was similar to the spatial distribution pattern of reference evapotranspiration in spring. Meanwhile, most trend coefficients fell between 0.01 and 0.3, and stations with a positive trend coefficient accounted for 85.5%. Although reference evapotranspiration in most areas revealed a rising trend, stations showing a downward trend were mainly distributed in the northwestern areas.
Areas with low annual average reference evapotranspiration coincided with that of the seasonal reference evapotranspiration, but large differences remained in high-value areas. Annual average reference evapotranspiration was between 1134 mm and 1387 mm, and reference evapotranspiration in the eastern region was higher than that in the western region, which coincided with the spatial distribution of reference evapotranspiration in autumn. Areas with lower annual average reference evapotranspiration were mainly distributed in the western mountains with higher altitudes, mainly due to the lower temperatures, higher relative humidity, smaller surface water evaporation and favorable moisture conditions [16], the annual average reference evapotranspiration in the hilly regions in southern China was mainly between 0.01 and 0.3, and stations with positive trend coefficients accounted for 84.1%, showing an overall rising trend. As areas with larger trend coefficients were mainly located in the western, northern [36] and northeastern regions, a significant rising trend could be observed in these areas. On the whole, reference evapotranspiration to the east of the hilly regions in southern China was higher than that in the western regions, and the reference evapotranspiration in the western mountains was relatively lower. In most areas, reference evapotranspiration showed a rising trend, and a gradual downward trend was observed in a handful of high-altitude areas including Mount Qixian and Mount Wuyi.
Water 2019, 11,1914 10 of 22 due to the lower temperatures, higher relative humidity, smaller surface water evaporation and favorable moisture conditions [16], the annual average reference evapotranspiration in the hilly regions in southern China was mainly between 0.01 and 0.3, and stations with positive trend coefficients accounted for 84.1%, showing an overall rising trend. As areas with larger trend coefficients were mainly located in the western, northern [36] and northeastern regions, a significant rising trend could be observed in these areas. On the whole, reference evapotranspiration to the east of the hilly regions in southern China was higher than that in the western regions, and the reference evapotranspiration in the western mountains was relatively lower. In most areas, reference evapotranspiration showed a rising trend, and a gradual downward trend was observed in a handful of high-altitude areas including Mount Qixian and Mount Wuyi.

Analysis of Abrupt Changes in Reference Evapotranspiration
To conduct further discussion on changes in reference evapotranspiration, the MK (Mann-Kendall) Test Method was applied to analyze abrupt changes in reference evapotranspiration in the hilly regions in southern China over the last 63 years, since the MK test analysis method is relatively accurate for the mean mutation test, it was selected as the mutation analysis method, and the temporal sub-sequence was 5 years. For an introduction to this method, see the paper by Wei et al. [37]. The average reference evapotranspiration calculated from the multi-year average reference evapotranspiration of all sites in the region is used to represent the reference evapotranspiration status of the entire region, and the MK test of the reference evapotranspiration is performed. Figure 6 shows the MK Test of ET 0 . UF referred to the sequential statistical curve (full line), and UB is the inverted sequential statistical curve (dotted line). If the value of UF or UB is greater than 0, the sequence indicates an upward trend, and less than 0 indicates a downward trend. The range in which they exceed the critical line is determined as the time zone in which the mutation occurs. If there is an intersection between the two curves of UF and UB, and the intersection is between the critical lines, then the moment corresponding to the intersection is the time at which the mutation begins [37] To conduct further discussion on changes in reference evapotranspiration, the MK (Mann-Kendall) Test Method was applied to analyze abrupt changes in reference evapotranspiration in the hilly regions in southern China over the last 63 years, since the MK test analysis method is relatively accurate for the mean mutation test, it was selected as the mutation analysis method, and the temporal sub-sequence was 5 years. For an introduction to this method, see the paper by Wei et al. [37]. The average reference evapotranspiration calculated from the multi-year average reference evapotranspiration of all sites in the region is used to represent the reference evapotranspiration status of the entire region, and the MK test of the reference evapotranspiration is performed. Figure  6 shows the MK Test of ET0. UF referred to the sequential statistical curve (full line), and UB is the inverted sequential statistical curve (dotted line). If the value of UF or UB is greater than 0, the sequence indicates an upward trend, and less than 0 indicates a downward trend. The range in which they exceed the critical line is determined as the time zone in which the mutation occurs. If there is an intersection between the two curves of UF and UB, and the intersection is between the critical lines, then the moment corresponding to the intersection is the time at which the mutation begins [37]   In order to further explore the temporal and spatial distribution of mutations in various meteorological stations, Figure 7 counts the number of mutations of reference evapotranspiration at each meteorological site and classifies the number of mutations by age. In total, 65 meteorological stations had shown abrupt changes in reference evapotranspiration, accounting for 94.2%. There was a significant fluctuation in reference evapotranspiration in most areas over the 63 years. The number of stations involved in such abrupt changes peaked in the 1950s and 1960s, indicating the significant sensitivity of reference evapotranspiration in most stations in the 1950s and 1960s. This result is consistent with the result of Figure 6 that the mutation of the reference evapotranspiration in the entire region occurred in 1953 and 1964. Spatially, stations with the most abrupt changes were mainly distributed in the eastern region, especially in the southeastern and northern areas in the vicinity of Anhui Province. Only a handful of areas located in the northeastern and northwestern regions were free from abrupt changes. In contrast to Figure 5e, areas with the most abrupt changes in reference evapotranspiration coincided with that with large annual average reference evapotranspiration, showing that areas with larger reference evapotranspiration were also sensitive to changes in reference evapotranspiration.

Periodic Analysis of Reference Evapotranspiration
Similar to various climatic and hydrologic elements, reference evapotranspiration has corresponding periodic features as well. Thus, this paper has applied the continuous wavelet to analyze and discuss periodic changes in reference evapotranspiration in the hilly regions in southern  (Figure 8e). In addition, the annual average reference evapotranspiration had experienced a relatively weak vibration in the end of the 1990s (since the 1970s), indicating that changes in reference evapotranspiration were stable during this period. This finding was similar to that of the annual average reference evapotranspiration shown in Figure 4b. In short, seasonal and annual average reference evapotranspiration had experienced 1.5-5.0 years of short-period changes. Based on the wavelet analysis on various meteorological variables, the periodic changes in reference evapotranspiration were in sync with that of relative humidity, proving that periodic changes in reference evapotranspiration were mainly caused by the short-period vibration of relative humidity.

Meteorological Variables Pearson Correlation Analysis
Judging from the definitions of reference evapotranspiration, various meteorological variables, which influence these values include precipitation, percentage of sunshine, average wind speed, relative humidity and temperature [38].Usually, the greater the percentage of sunshine, the longer the duration of actual sunshine; the higher the temperature, the faster the evapotranspiration of surface moisture. Similarly, an increase in wind speed would also accelerate evapotranspiration and may create significantly more moisture as well. Meanwhile, an increase in relative humidity would lead to a corresponding increase in surface water and moisture in the air, but reduce reference evapotranspiration [39]. As such, the Pearson correlation analysis was used to analyze the correlation between reference evapotranspiration and various meteorological elements [40] and r value was between [−1,1]. When r > 0, it indicates positive correlation, while r < 0 indicates negative correlation. The ranges of 0.3 > |r| ≥ 0, 0.5 > |r| ≥ 0.3, 0.8 > |r| ≥ 0.5 and |r| ≥ 0.8 indicate weak correlation, medium correlation, strong correlation and extremely strong correlation respectively.
Based on Table 1, it is evident that precipitation has a negative correlation with reference evapotranspiration in the hilly regions in southern China on an annual and seasonal basis (p < 0.01). In addition, relative humidity was negatively correlated with reference evapotranspiration, and passed the 0.01 significance test. Although a weak positive correlation was observed on an annual basis and in autumn, it did not pass the significance test. In general, the wind speed was positively correlated with reference evapotranspiration, but, this correlation was not significant. Temperature was positively correlated with reference evapotranspiration. With the exception of an annual basis and in winter, sunshine percentage was positively correlated with reference evapotranspiration, and passed the 0.01 significance test.

Meteorological Variables Pearson Correlation Analysis
Judging from the definitions of reference evapotranspiration, various meteorological variables, which influence these values include precipitation, percentage of sunshine, average wind speed, relative humidity and temperature [38]. Usually, the greater the percentage of sunshine, the longer the duration of actual sunshine; the higher the temperature, the faster the evapotranspiration of surface moisture. Similarly, an increase in wind speed would also accelerate evapotranspiration and may create significantly more moisture as well. Meanwhile, an increase in relative humidity would lead to a corresponding increase in surface water and moisture in the air, but reduce reference evapotranspiration [39]. As such, the Pearson correlation analysis was used to analyze the correlation between reference evapotranspiration and various meteorological elements [40] and r value was between [−1,1]. When r > 0, it indicates positive correlation, while r < 0 indicates negative correlation. The ranges of 0.3 > |r| ≥ 0, 0.5 > |r| ≥ 0.3, 0.8 > |r| ≥ 0.5 and |r| ≥ 0.8 indicate weak correlation, medium correlation, strong correlation and extremely strong correlation respectively.
Based on Table 1, it is evident that precipitation has a negative correlation with reference evapotranspiration in the hilly regions in southern China on an annual and seasonal basis (p < 0.01). In addition, relative humidity was negatively correlated with reference evapotranspiration, and passed the 0.01 significance test. Although a weak positive correlation was observed on an annual basis and in autumn, it did not pass the significance test. In general, the wind speed was positively correlated with reference evapotranspiration, but, this correlation was not significant. Temperature was positively correlated with reference evapotranspiration. With the exception of an annual basis and in winter, sunshine percentage was positively correlated with reference evapotranspiration, and passed the 0.01 significance test.  Table 2 shows the relationship between the sensitivity coefficient of each meteorological variables, the amount of change over the years, and the contribution. Reference evapotranspiration has a higher sensitivity coefficient to the relative humidity and average temperature. Except to relative humidity, the sensitivity coefficient of reference evapotranspiration to temperature, wind speed and sunshine hours are positive. Many scholars have also proved that the sensitivity coefficient of reference evapotranspiration to relative humidity is usually negative in China [33,[41][42][43]. Although the annual variation of average temperature and relative humidity is much smaller than the sunshine hours and wind speed, the reference evapotranspiration has higher sensitivity coefficient to the average temperature (positive effect) and relative humidity (negative effect), while the sensitivity coefficient to the other two meteorological variables (positive effect) are too small. As a result, the highest contribution is the average temperature, followed by the relative humidity. Many scholars have also shown that relative humidity has a great contribution to reference evapotranspiration in most areas of China, but the contribution of temperature to reference evapotranspiration is small in most parts of China as different locations, the contribution of wind speed to reference evapotranspiration is greater than temperature in Northeast China and Northern China [41,42]. The average temperature and relative humidity are positive contributions to the reference evapotranspiration. It is similar to the Anhui province of China, which is closer to the study area [33]. The sunshine hours and wind speeds are negative contributions to the reference evapotranspiration, and the negative contributions of sunshine hours to the reference evapotranspiration is much greater than the wind speed. The total contribution of the four main variables affecting the reference evapotranspiration to the change of reference evapotranspiration is 7.84%, while the average annual change of reference evapotranspiration in the past 63 years is 11.13%, indicating that the change of other variables may contribute positively (3.288%) to the reference evapotranspiration. Note: T means average temperature; RH means relative humidity; n means sunshine hours; U means wind speed. Figure 9a shows the spatial ratio of the sensitivity coefficient of reference evapotranspiration to meteorological variables (average temperature, relative humidity, sunshine hours and wind speed) in the hilly area of Jiangnan. In space, the meteorological variables with strong sensitivity are temperature, relative humidity, sunshine hours and wind speed, but the sensitivity coefficient of ET 0 to wind speed is extremely low. Therefore, the spatial distribution of sensitivity coefficient of ET 0 to average temperature, relative humidity and sunshine hours is mainly considered. The sensitivity coefficient of ET 0 to the average temperature in the study area is positive (Figure 9b), shows that the reference evapotranspiration increases with the increase of the average temperature. Among them, the sensitivity coefficient of ET 0 to the average temperature in the north is generally larger than that in the south, and the sensitivity coefficient of ET 0 to the average temperature in the east is significantly larger than that in the west, indicating reference evapotranspiration is more sensitive to average temperature in the north and east. The sensitivity coefficient of ET 0 to relative humidity in the study area is negative (Figure 9c), indicating that the reference evapotranspiration decreases with the increase of relative humidity. The sensitivity coefficient' absolute value of ET 0 to the relative humidity in the northwestern region is generally smaller than other regions, indicating that the reference evapotranspiration in the northwestern has lower sensitivity with relative humidity (RH). The sensitivity coefficient of ET 0 to sunshine hours in the study area is positive (Figure 9d), the spatial distribution of the sensitivity coefficient of ET 0 to sunshine hours is different from the spatial distribution of average temperature and relative humidity, and the spatial distribution of sensitivity coefficient of ET 0 to sunshine hours is uniform, indicating that the sensitivity of reference evapotranspiration to sunshine hours is less affected by spatial location.

Sensitivity Coefficient and Contribution
Water 2019, 11,1914 15 of 22 temperature, relative humidity, sunshine hours and wind speed, but the sensitivity coefficient of ET0 to wind speed is extremely low. Therefore, the spatial distribution of sensitivity coefficient of ET0 to average temperature, relative humidity and sunshine hours is mainly considered. The sensitivity coefficient of ET0 to the average temperature in the study area is positive (Figure 9b), shows that the reference evapotranspiration increases with the increase of the average temperature. Among them, the sensitivity coefficient of ET0 to the average temperature in the north is generally larger than that in the south, and the sensitivity coefficient of ET0 to the average temperature in the east is significantly larger than that in the west, indicating reference evapotranspiration is more sensitive to average temperature in the north and east. The sensitivity coefficient of ET0 to relative humidity in the study area is negative (Figure 9c), indicating that the reference evapotranspiration decreases with the increase of relative humidity. The sensitivity coefficient' absolute value of ET0 to the relative humidity in the northwestern region is generally smaller than other regions, indicating that the reference evapotranspiration in the northwestern has lower sensitivity with relative humidity (RH). The sensitivity coefficient of ET0 to sunshine hours in the study area is positive (Figure 9d), the spatial distribution of the sensitivity coefficient of ET0 to sunshine hours is different from the spatial distribution of average temperature and relative humidity, and the spatial distribution of sensitivity coefficient of ET0 to sunshine hours is uniform, indicating that the sensitivity of reference evapotranspiration to sunshine hours is less affected by spatial location.  Figure 10a shows the spatial ratio of the contribution of meteorological variables (average temperature, relative humidity, sunshine hours and wind speed) in the hilly area of Jiangnan. Affected by the size of the sensitivity coefficient, the meteorological variables with a large contribution to the reference evapotranspiration are average temperature, relative humidity, sunshine hours, and wind speed. Since the wind speed has a very small sensitivity coefficient to reference evapotranspiration, the contribution of the wind speed to the reference evapotranspiration  Figure 10a shows the spatial ratio of the contribution of meteorological variables (average temperature, relative humidity, sunshine hours and wind speed) in the hilly area of Jiangnan. Affected by the size of the sensitivity coefficient, the meteorological variables with a large contribution to the reference evapotranspiration are average temperature, relative humidity, sunshine hours, and wind speed. Since the wind speed has a very small sensitivity coefficient to reference evapotranspiration, the contribution of the wind speed to the reference evapotranspiration is also small. Therefore, the spatial distribution of the average air temperature, relative humidity and sunshine hours to the reference evapotranspiration contribution is mainly analyzed. The contribution of the average temperature to the reference evapotranspiration is positively contributed by most sites, while the sensitivity coefficient of the average temperature is positive in Figure 10b. According to the Equation, the relative change of temperature in many sites is positive. The spatial distribution law of the contribution of average temperature is similar to the sensitivity coefficient of average temperature. The contribution of the north and east is greater than that of the south and the west, and the contribution of the southwest is the smallest. Relative humidity has a negative contribution to the reference evapotranspiration in most sites. Combined Figure 10c with Equation (10), the sensitivity coefficient of relative humidity is negative, indicating that the relative humidity of the study area has had a positive change over the years. The spatial distribution of the contribution of relative humidity was different from its sensitivity coefficient, and the contribution of relative humidity in the southeast is lower. The contribution of sunshine hours to reference evapotranspiration in space is mostly negative, while the sensitivity coefficient of sunshine hours in Figure 10d is positive, shows that the average change of sunshine hours is negative. The spatial distribution of the contribution of sunshine hours is far less uniform than the spatial distribution of its sensitivity coefficients, and the contribution of sunshine hours in the southwest is the smallest, mean that the average change of sunshine hours in the southwestern region also smallest.
Water 2019, 11, 1914 16 of 22 is also small. Therefore, the spatial distribution of the average air temperature, relative humidity and sunshine hours to the reference evapotranspiration contribution is mainly analyzed. The contribution of the average temperature to the reference evapotranspiration is positively contributed by most sites, while the sensitivity coefficient of the average temperature is positive in Figure 10b. According to the Equation, the relative change of temperature in many sites is positive. The spatial distribution law of the contribution of average temperature is similar to the sensitivity coefficient of average temperature. The contribution of the north and east is greater than that of the south and the west, and the contribution of the southwest is the smallest. Relative humidity has a negative contribution to the reference evapotranspiration in most sites. Combined Figure 10c with Equation (10), the sensitivity coefficient of relative humidity is negative, indicating that the relative humidity of the study area has had a positive change over the years. The spatial distribution of the contribution of relative humidity was different from its sensitivity coefficient, and the contribution of relative humidity in the southeast is lower. The contribution of sunshine hours to reference evapotranspiration in space is mostly negative, while the sensitivity coefficient of sunshine hours in Figure 10d is positive, shows that the average change of sunshine hours is negative. The spatial distribution of the contribution of sunshine hours is far less uniform than the spatial distribution of its sensitivity coefficients, and the contribution of sunshine hours in the southwest is the smallest, mean that the average change of sunshine hours in the southwestern region also smallest.

Climatic Indicators
Reference evapotranspiration was affected by various variables. Driven by various meteorological variables, they also have a teleconnection with climatic indexes. From Table 3, a significant correlation existed between WP, SOI, TNA, and ENSO, and the reference evapotranspiration on an annual and seasonal basis.

Climatic Indicators
Reference evapotranspiration was affected by various variables. Driven by various meteorological variables, they also have a teleconnection with climatic indexes. From Table 3, a significant correlation existed between WP, SOI, TNA, and ENSO, and the reference evapotranspiration on an annual and seasonal basis. Reference evapotranspiration was negatively correlated with WP on an annual and seasonal basis due to WP being closely related to the winter monsoon in the Northern Hemisphere. The stronger the WP in winter and spring, the lower the temperature, the larger the precipitation, the smaller the reference evapotranspiration [32,44]. Reference evapotranspiration and SOI were positively correlated on an annual and seasonal basis, and the significance test was passed on annual basis in winter, indicating that the stronger the SOI, the larger the reference evapotranspiration. Reference evapotranspiration was positively correlated with TNA on an annual and seasonal basis, and the significance test was passed on an annual basis, in spring and autumn, due to the positive correlation of TNA in summer with the South China Sea monsoon [44]. Reference evapotranspiration was negatively correlated with ENSO on an annual and seasonal basis, and the significance test was only passed in winter due to the significant increase in precipitation during most seasons in most parts of Southern China during El-Nino [45]; precipitation was negatively correlated with reference evapotranspiration. The correlation between reference evapotranspiration and ENSO was contrary to that between reference evapotranspiration and SOI indexes as SOI indexes were negative during El-Nino, and El-Nino3. 4 Zone was negatively correlated with SOI indexes [45].

Other Drivers
Besides meteorological elements and climatic factors, the geographical distribution of surface water was also affected by the geographic location. By screening various geographic factors related to surface water and applying the multivariate regression analysis, the prediction models of annual and seasonal reference evapotranspiration in the hilly regions in southern China were created [46]. All five models had passed the significance test with confidence coefficient at α = 0.01.
Analyses of the models showed that annual and seasonal reference evapotranspiration was negatively correlated with the altitude, indicating that a rise in altitude would cause a decline in reference evapotranspiration. As for latitude, a negative correlation could be observed in all seasons other than spring. This means that an increase in latitude would cause a decline in reference evapotranspiration as higher latitudes would lead to lower temperatures with lower evapotranspiration rates. With regards to longitude, a weak positive correlation was found, and an increase in longitude would cause an increase in reference evapotranspiration (Table 4). These laws are consistent with the contents shown in Figure 5. The reference evapotranspiration is relatively lower in areas with high altitudes; the reference evapotranspiration in the south is higher than in the north in the spring, summer, winter, and year; the seasonal and annual average reference evapotranspiration are higher in the east than in the west. Figure 11 shows the Variation Fitting for Annual and Seasonal Simulated Reference Evapotranspiration and Actual Reference Evapotranspiration with Latitude, Longitude and Altitude (69 stations were numbered from the lowest latitude, longitude and altitude to the highest as the horizontal axis). Most estimated values of annual and seasonal reference evapotranspiration fell within the relative error range (±20% of the actual measured values), proving that it is better to use regression distribution models to simulate changes in reference evapotranspiration with latitude, longitude and altitude. However, the annual fitting equation could achieve the best effect, while the effect of such a fitting equation in winter was slightly poor. Compared to latitude and longitude, altitude would almost exert a synchronized influence on reference evapotranspiration both on a seasonal and an annual basis. On the contrary, reference evapotranspiration would decline as altitude rose. In particular, a sharp decline was observed from No.58 station (altitude above 350 m). influence on reference evapotranspiration both on a seasonal and an annual basis. On the contrary, reference evapotranspiration would decline as altitude rose. In particular, a sharp decline was observed from No.58 station (altitude above 350 m).

Summary and Conclusions
Between 1951 and 2013, changes in the annual average reference evapotranspiration in hilly regions in southern China had featured a gradual fluctuating rise-gradual fluctuating decline-rapid rise pattern. However, an overall gradual rising trend could be observed; reference evapotranspiration reached its peak in summer and its lowest point in winter, and the rising trend was extremely significant in spring. Reference evapotranspiration to the east of the hilly regions in southern China was higher than that in the west, a rising trend could be observed in most areas and the reference evapotranspiration in a handful of high-altitude mountains was relatively lower with a gradual downward trend.
In 1953, 1966 and 2008, abrupt changes occurred in reference evapotranspiration. In fact, 94.2% of the meteorological stations were involved in such abrupt changes in reference evapotranspiration in the 1950s and 1960s; As the spatial distribution of stations involved in the abrupt changes in reference evapotranspiration was different, stations that were involved in the most abrupt changes in reference evapotranspiration were mainly distributed in the eastern region. Reference evapotranspiration was affected by periodic changes in relative humidity, and showed short-period vibration in most cases.
Various meteorological variables had played a driving role in reference evapotranspiration. Amongst these variables, relative humidity was negatively correlated with reference evapotranspiration on the whole, while temperature and sunshine percentage were positively correlated with reference evapotranspiration. The meteorological variables that have a great influence on the sensitivity of reference evapotranspiration are average temperature, relative humidity, sunshine hours, and wind speed. The sensitivity coefficient of reference evapotranspiration to temperature, sunshine hours and wind speed are positive, while the sensitivity coefficient to relative humidity is negative. It is roughly consistent with the correlation of reference evapotranspiration and various meteorological variables. The sensitivity coefficient of reference evapotranspiration to average temperature in the northern and eastern regions are larger than that in the south and west, and the sensitivity coefficient to relative humidity in the northwest is the smallest. However, the spatial distribution of the sensitivity coefficient of reference evapotranspiration to sunshine hours is relatively uniform, its less affected by spatial location. The spatial distribution of the contribution is similar to sensitivity coefficient, however, the spatial distribution of the contribution of sunshine hours is far less uniform than its sensitivity coefficient due to the difference of relative changes over the years. WP and ENSO are positively correlated with reference evapotranspiration, while SOI and TNA are negatively correlated with reference evapotranspiration. In addition, a negative correlation, weak positive correlation and negative correlation could be observed between the reference evapotranspiration and the latitude/longitude/altitude respectively.