Quantifying Climatic Impact on Reference Evapotranspiration Trends in the Huai River Basin of Eastern China

Reference evapotranspiration (ETref) is an important study object for hydrological cycle processes in the context of drought-flood risks of the Huai River Basin (HRB). In this study, the FAO-56 Penman–Monteith (PM) model was employed to calculate seasonal and annual ETref based on 137 meteorological station data points in HRB from 1961 to 2014. The Mann–Kendall (MK) trend analysis was adopted together with Theil–Sen’s estimator to detect tendencies of ETref and climate factors. Furthermore, a developed differential equation method based on the FAO-56 PM model was applied to quantify the sensitivities of ETref to meteorological factors and their contributions to ETref trends. The results showed that the ETref demonstrated a strong spatially heterogeneity in the whole HRB at each time scale. ETref showed a significant decreasing trend in the upper-middle HRB and Yi-Shu-Si River Basin, especially at the annual time scale, in growing season and summer, while a generally increasing trend in ETref was detected in the lower HRB, and the significance only showed in spring. These phenomena could be reasonably explained by a significantly increasing mean temperature (TA), a significantly decreasing wind speed (WS), solar radiation (SR), and a slightly decreasing relative humidity (RH). The most sensitive factor to ETref was RH in most sub-regions and most time scales, except in the growing season and summer. Based on the developed differential equation method, the dominant factor of the decreasing ETref was WS in the annual time scale, spring, autumn, and winter in most sub-regions, except the lower HRB, which then shifted to SR in the growing season and summer. However, in the lower HRB, the significantly decreasing RH was the most dominant factor, especially in the annual time scale, growing season, and spring, which might be responsible for the slightly increasing ETref there.


Introduction
Evapotranspiration (ET), as a vital component in water cycle, is a critical variable for hydrometeorological studies [1].At the same time, the reference evapotranspiration (ET ref ) refers to the evapotranspiration (ET) from the crop reference surface, which has a 70 s•m −1 unified surface resistance, a 0.12 m hypothetical reference height, and a 0.23 albedo [2,3].It determines the surface runoff, soil moisture, groundwater recharge, as well as the crop growth [4], and plays a crucial role in climatology, hydrological processes, irrigation management, planning, and scheduling [5,6].Moreover, the global average temperature across the land and ocean surface increased by 0.85 • C from 1880 to 2012 [7].In the context of global warming, the hydrological cycle processes were altered dramatically in many parts of the world.Furthermore, crop productivity was also influenced by the climate change to some extent [8].
Nowadays, several methods have been developed to compute the ET ref presently, mainly based on the water budget, mass-transfer, radiation, and temperature [9].Among these methods, the FAO-56 PM model, as proposed by the Food and Agriculture Organization (FAO), is universally accepted as the most precise one to estimate the reference evapotranspiration [2].However, as global warming has become a historical truth and scientific fact in recent decades, the ET ref , as well as pan evaporation (E pan ), exhibited a significant downward trend globally [10][11][12].Meanwhile, the decreasing trend of ET ref or E pan generally appears to be an "evaporation paradox" [13], one of the most controversial scientific problems.Likewise, decrease of ET ref trends has been discovered from many areas in China [14][15][16].Apart from these reports, many researchers have shown the increase of ET ref in various places around the world, including Asia [17,18] and Europe [19,20].These results are echoed with similar observations in China; for instance, the upper, middle, and entire Yellow River Basin [21], Jianghuai area, and Sunan area in Jiangsu province [16], and most parts of Wei River Basin [22].Additionally, during the past several decades, a significant number of investigators have identified the influence of meteorological parameters to the changing ET ref [23][24][25].ET ref variation, and its spatiotemporal distribution, are not only influenced by the increasing air temperature, but also by some primary climate factors, namely, relative humidity, wind speed, and solar radiation.Therefore, it is an important study object to quantify the influences of various climate factors on the ET ref change, especially in the background of global climatic change [26].
Different attribution analysis methods are widely employed for calculating contributions of climate factors to ET ref variations.The sensitivity coefficients method, multiple regressions, and the detrending method are the three most used attribution analysis methods.Mccuen [27] first introduced the sensitivity coefficient method to recognize the relative change of ET ref according to meteorological factors.Later, Yin et al. [28] and Huo et al. [14] successfully applied this approach to identify the anticipated change of ET ref caused by variations of climate factors.Similarly, other scholars in China revealed relative contributions of climate variables to the ET ref by employing sensitivity coefficient analysis [29][30][31].Contrarily, Ye et al. [32], Shan et al. [33], and Wang et al. [34] adopted the multiple regression analysis method to establish the relations between the ET ref (dependent variable) and climate parameters (independent variables), then determined relative contributions of climatic variables to the ET ref trends.Recently, the stepwise regression analysis and partial correlation analysis were employed for investigating influences of climate factors to ET ref change in the West Liao River Basin [35].Moreover, the detrending method, a statistical and mathematical modeling approach, first proposed by Xu et al. [1], was used to evaluate changes of ET ref to the related climate factors.This method has been applied to assess the influences of climate factors to the ET ref change in various regions of China such as the Heihe River Basin [36], Yellow River Basin [37], northwest area of China [14], Jiangsu province of Eastern China [16], and also other countries like Iran [38].The above-mentioned attribution analysis methods only consider the relative contributions of meteorological factors to ET ref changes, however, actual contributions to the ET ref changes are still unclear and needed to be explored.Roderick et al. [39] have pioneered the use of a physical and mathematical model-namely, the PenPan model-which is based on mass and energy balances, to evaluate the contributions of climatic parameters to E pan .In the PenPan model, the E pan trend is divided into the radiative and aerodynamic items.Furthermore, the E pan trend caused by the aerodynamic item is also partitioned into three components: wind speed, vapor pressure deficit, and mean air temperature.Various researchers have adopted this differential equation method to quantitatively analysis the change in E pan [40][41][42], however, few studies have used it for ET ref analysis.Additionally, the PenPan model was employed by Liu and Zhang [43] to determine the changes in ET ref in Northwest China.Some previous studies [40,43] also adopted the differential equation method based on the FAO-56 PM model to determine contributions of four meteorological variables (mean temperature, actual vapor pressure, solar radiation, and wind speed) to the ET ref trend.However, the actual vapor pressure is not an independent variable, but a dependent variable on mean air temperature and relative humidity.Thus, considering the independence of each selected variable, we decided to replace the actual vapor pressure with the relative humidity in this study, and hypothesize that the combination of the four new variables is reasonable, then we verify this hypothesis.Considering the accuracy and reliability of the differential equation method, we developed this method based on the FAO-56 PM model to verify our hypothesis, and also to quantify the attribution of each climatic factor to ET ref trends in the Huai River Basin (HRB).
HRB, adjoining the Yangtze River Delta, as a fast economic growth region in Eastern China, is susceptible to confronting extreme weather and hydro-meteorological events, particularly flood calamities [44].These climatic challenges will exacerbate a potential impact on the hydrometeorology of the HRB.Moreover, the HRB is in a decisive status in the agricultural production of China and the density of population here is nearly five-fold higher than the national average level.However, there are only a few studies about long-term quantitative attribution analysis of meteorological factors to the ET ref in the HRB of Eastern China.Thus, research on the climate factors attribution analysis to regional ET ref changes is needed to understand the spatiotemporal trends of regional climatic change and the influence of climatic change on agricultural irrigation scheduling.
Based on the hypotheses above, the aims of this study are (1) to estimate the ET ref changes in seasonal and annual time scales in HRB from 1961 to 2014 by means of the FAO-56 PM model; (2) to identify the spatiotemporal tendencies of seasonal and annual ET ref series together with its relative climate factors through the Mann-Kendall (MK) test, together with Sen's slope estimator, over the past 54 years; (3) to explore the sensitivities of ET ref to each climatic factor at seasonal and annual time scales using sensitivity coefficients analysis; and (4) to quantify the contributions from climate factors to ET ref change trends by using a developed differential equation model and its spatial distribution patterns in seasonal and annual time scales.This study would offer a theoretical guidance for the management of agricultural irrigation resource usage, crop production, as well as the countermeasure to climatic change perspective in this vital region.

Study Area
The HRB is located between the Yangtze River and Yellow River Basin in Eastern China (111 • 55 E-121 • 25 E and 30 • 55 N-36 • 36 N) and has a watershed area of 270,000 km 2 .The length of the main stream is about 1000 km and the total river fall is about 200 m (Figure 1).Meanwhile, the HRB is situated in the climatic transition zone, which is recognized as the demarcation line (zero isotherm) between warm temperate and northern subtropical zones in China, accompanied with four distinctive seasons and moderate climate.The annual average precipitation is about 911 mm, with an uneven distribution pattern.More than 50% of the rainfall is concentrated during the monsoon season from June to September, accompanied by complex weather systems.Furthermore, the HRB is the most densely populated area and also a major grain-producing base in Eastern China [45], with the main crops-including wheat, rice, maize, potato, soybean, cotton, and canola-accounting for 17.3% of the total national grain yield.

Data Source
The datasets employed in the present study are mainly daily weather data obtained from 137 stations during 1961-2014 in the HRB, including the daily mean temperature (TA, • C), maximum temperature (T max , • C) and minimum temperature (T min , • C), relative humidity (RH, %), wind speed at 10 m (u 10 , m•s −1 ), and sunshine duration (SD, h).We obtained the weather data from the National Meteorological Information Center (NMIC) of the China Meteorological Administration (CMA) (http://data.cma.cn/site/index.html).The censored data are complemented through establishing linear interpolation between the neighboring weather stations.In this study, we defined five seasons as the growing season (April-October), spring (March-May), summer (June-August), autumn (September-November), and winter (December-February of the following year).The DEM (digital elevation model) datasets of 90 m spatial resolution (Figure 1) were acquired from the Internet (http://srtm.csi.cgiar.org/).
Water 2018, 10, x FOR PEER REVIEW 4 of 23 (http://data.cma.cn/site/index.html).The censored data are complemented through establishing linear interpolation between the neighboring weather stations.In this study, we defined five seasons as the growing season (April-October), spring (March-May), summer (June-August), autumn (September-November), and winter (December-February of the following year).The DEM (digital elevation model) datasets of 90 m spatial resolution (Figure 1) were acquired from the Internet (http://srtm.csi.cgiar.org/).

Reference Evapotranspiration Calculation
The FAO-56 Penman-Monteith (PM) model proposed by the Food and Agriculture Organization (FAO) has been regarded as an international recognized method to calculate ETref [2].The FAO-56 PM model has a full theoretical basis which shows high accuracy [16].The detailed formula is given as where ET ref refers to the daily reference crop evapotranspiration (mm•d −1 ), ∆ refers to the slope of vapor pressure curve (kPa

Reference Evapotranspiration Calculation
The FAO-56 Penman-Monteith (PM) model proposed by the Food and Agriculture Organization (FAO) has been regarded as an international recognized method to calculate ET ref [2].The FAO-56 PM model has a full theoretical basis which shows high accuracy [16].The detailed formula is given as where ET re f refers to the daily reference crop evapotranspiration (mm•day −1 ), ∆ refers to the slope of vapor pressure curve (kPa• • C −1 ), R n refers to the net radiation at the crop height (MJ•m −2 •day −1 ), G refers to the soil heat flux density (MJ•m −2 •day −1 ), γ refers to the psychrometric constant (kPa• • C −1 ), T refers to the daily mean air temperature at 2 m ( • C), u 2 refers to the wind speed at two-meter height (m•s −1 ), e s and e a refer to the saturation vapor pressure (kPa) and the actual vapor pressure (kPa), respectively.For the convenience of ET re f calculation, a conversion formula for wind speed was proposed by Allen et al. [2] as where u 2 refers to wind speed at two-meter height (m•s −1 ), u z refers to the wind speed at z meters above ground level.
where R s refers to daily total solar radiation (MJ•m −2 •day −1 ), R a refers to extraterrestrial radiation (MJ•m −2 •day −1 ), n refers to daily sunshine hours (h), N refers to maximum possible sunshine hours (h), a s and b s are regression coefficients, according to Chen et al. [47], the values were set to 0.19 and 0.53, respectively [16].
Meanwhile, some sunshine duration data were missing in several inland stations of the HRB.In this regard, the radiation equation proposed by Hargreaves and Samani [48] was employed for calculating solar radiation (R s ) through T max and T min as where k RS refers to the empirical coefficient and usually adopted the value of 0.16 for inland regions [48][49][50].In order to verify the availability of the Hargreaves-Samani method, complete sunshine duration (SD) datasets of 39 meteorological stations located in the upper and middle HRB were employed in this study.The estimated SR and ET ref calculated with Hargreaves-Samani (HS) method were compared with that calculated with SD during 1961-2014.The specific accuracy test can be found in Appendix A.
Finally, the net radiation (R n ) is estimated by in which R ns refers to the incoming net shortwave radiation (MJ•m −2 •day −1 ), α refers to the albedo of reference crop (value of 0.23), R nl refers to the net outgoing long-wave radiation (MJ•m −2 •day −1 ), σ refers to the Stefan-Boltzmann constant (4.903 × 10 −9 MJ•K −4 •m −2 •day −1 ), T max,K represents the maximum absolute temperature during 24 h (K = • C + 273.16),T min,K represents the minimum absolute temperature during 24 h (K = • C + 273.16), and R so refers to clear-sky radiation.Furthermore, more specific approach procedures are provided in Allen et al.'s literature [2].

Trend Analysis
The worldwide recommended nonparametric test, Mann-Kendall (MK) test, was built to calculate the significance of a trend in hydro-meteorological time series [51,52].The MK test is employed in this research to analyze temporal change trends of ET ref and related meteorological variables [53][54][55].It has an advantage which lies in its ability to test the linearity of a trend [31].
In order to identify slopes of ET ref trends and climate factors, Theil-Sen's estimator [55,56], a widely applied method to calculate the slope magnitude of the trend line, was adopted in this study where β refers to the estimated trend slope of the data series, x j and x i refer to the sequential data corresponding to time j and i, respectively.Prior to applying the M-K test in this study, a pre-whitening method should be employed to eliminate the effect of a serial correlation on the trend value [57].Yue et al. [58] improved this method Water 2018, 10, 144 6 of 24 and proposed an alternative approach, namely TFPW (trend-free pre-whitening), which is widely applied in many related study areas [17,30,59].

Spatial Interpolation Method
For investigating the spatial patterns of ET ref and their trends, some spatial interpolation methods, such as spline, ordinary kriging (OK), and the inverse distance weighted (IDW), were implemented on account of datasets of meteorological stations being widely employed recently.Among these interpolation approaches, the IDW approach was adopted for the spatial pattern analysis in our study since it provides the lower mean error than the others [60,61].The main justification for using the IDW technique is that it computes the spatially-interpolated values very quickly and accurately.All spatial distribution maps were prepared using the ArcGIS (version 9.3) software (Environmental Systems Research Institute (ESRI), Redlands, CA, USA).

Sensitivity Coefficient Estimation
For identifying sensitivities of ET ref to meteorological factors, the differential method proposed by McCuen [27] is applied for this study where S(v i ) refers to the sensitivity coefficient of ET re f in relation to climatic variables (v i ).
The positive (negative) sensitivities illustrate the ET re f increases (decreases) are accompanied with the increases in v i .The absolute value of S(v i ) represents the effect of a given v i on ET re f .The specific derived calculation formula can be found from other literature [16].

Contributions of Climate Factors on ET ref
The differential equation method indicates that the ET ref is the function of meteorological variables.In previous studies [40,43], four meteorological variables (mean temperature (TA, • C), actual vapor pressure (e a , kPa), wind speed at two-meter height (WS, m•s −1 ), solar radiation (SR, MJ•m −2 •day −1 )) were selected as main influence factors to ET ref , however, e a was not an independent variable, but a dependent variable function of TA and relative humidity (RH, %), considering the independence of each meteorological variable, we decided to replace the e a by RH in this study.Thus, the improved the differential equation method was used to estimate the individual contribution of each climate factor (TA, RH, WS, and SR) to ET ref change trend based on the FAO-56 PM model.The detailed differential equation is or, written briefly, where C_ET re f refers to long term estimated trend of ET re f , C(TA), C(RH), C(WS), and C(SR) are contributions to ET re f due to changes in TA, RH, WS, and SR, respectively.ε refers to the error term between the C_ET re f and ET re f trends estimated by Theil-Sen's estimator (T_ET re f ).The contribution of each selected climate factor to the ET re f trend can be obtained by calculating the averaged partial derivatives (Equations ( 9) and ( 10)) multiplied by the annual/seasonal tendency rates and the length of time (namely 365 or 366 for annual, 214 for growing season, 92 for spring and summer, 91 for autumn, and 90 or 91 for winter).

Variations of Climate Factors and Their Trend Analysis
Variations of monthly climate factors are shown in Figure 2. From Figure 2, monthly TA and SR exhibited similar variation patterns.The peak value occurred in June and July.The largest regional differences of SR were shown in spring and summer seasons, while the largest regional difference of TA exhibited in winter.Meanwhile, RH gradually increased before August and then decreased.Monthly WS values represent a fluctuating variation trend, the peak and valley values occurred in April and September, respectively.
The effects evaluation of climate factors to ET ref was implemented by analyzing temporal trends of four meteorological parameters from 1961 to 2014 in the HRB.As from Figure 3, TA showed an overall upward trend in the past 54 years.It was noteworthy that a slow decreasing trend was observed before the 1990s and then a significantly increasing trend was observed.The RH fluctuated obviously, with a non-significant declining trend in the whole HRB.The WS and SR exhibited more significant declining trends in comparison with TA and RH.
The similar results can also be detected from Table 1.At the seasonal time scale, TA revealed a significant increasing trend in each season except the summer.Opposite to TA, RH exhibited a reverse change trend.Meanwhile, significant decreasing RH trends were found in the lower HRB in each season compared to the other regions.WS always exhibited a significant downward trend in all regions at each time scale, which could also be verified from Figure 3.The SR revealed a significant decreasing trend in the whole HRB in all time periods except spring.The similar results can also be detected from Table 1.At the seasonal time scale, TA revealed a significant increasing trend in each season except the summer.Opposite to TA, RH exhibited a reverse change trend.Meanwhile, significant decreasing RH trends were found in the lower HRB in each season compared to the other regions.WS always exhibited a significant downward trend in all regions at each time scale, which could also be verified from Figure 3.The SR revealed a significant decreasing trend in the whole HRB in all time periods except spring.

Variations of ET ref and Their Trend Analysis
In order to obtain insights into the spatial distributions of annual and seasonal ET ref in the HRB, the annual averaged ET ref in each station from 1961 to 2014 was spatially interpolated using the IDW method (Figure 4).The annual ET ref during the past 54 years in the HRB ranged from 881 mm to 1091 mm, with an average value of 969.56 mm.The higher values were observed in the north of the middle and Yi-Shu-Si River basins.Similar spatial distributions of ET ref were identified during the growing season and spring, ranging from 689 mm to 863 mm, and from 245 mm to 345 mm, respectively.Meanwhile, the lower values of ET ref mainly distributed in the upper and lower HRB.However, in the summer, the ET ref was much higher than other seasons (336~427 mm), and the higher ET ref values expanded to most parts of the upper and middle HRB also the Yi-Shu-Si River Basin.In autumn and winter, differences of ET ref were smaller in the whole HRB.The higher values were observed in the southeast in autumn, and in the northwest in winter.In general, the ET ref demonstrated a strong spatial heterogeneity in the whole HRB.respectively.Meanwhile, the lower values of ETref mainly distributed in the upper and lower HRB.However, in the summer, the ETref was much higher than other seasons (336~427 mm), and the higher ETref values expanded to most parts of the upper and middle HRB also the Yi-Shu-Si River Basin.In autumn and winter, differences of ETref were smaller in the whole HRB.The higher values were observed in the southeast in autumn, and in the northwest in winter.In general, the ETref demonstrated a strong spatial heterogeneity in the whole HRB.As shown in Table 2, at the annual time scale, a significant downward trend of ETref was found in the whole HRB (−1.894 mm•a −2 ), the upper and middle HRB and the Yi-Shu-Si River Basin.A non-significant upward trend was detected in the lower HRB (0.177 mm•a −2 ).The magnitude of the ETref trend followed the order: Middle > Upper > Yi-Shu-Si > Lower.During the growing season, the magnitude of the change rate was slightly lower than that at the annual time scale, while the spatial As shown in Table 2, at the annual time scale, a significant downward trend of ET ref was found in the whole HRB (−1.894 mm•a −2 ), the upper and middle HRB and the Yi-Shu-Si River Basin.A non-significant upward trend was detected in the lower HRB (0.177 mm•a −2 ).The magnitude of the ET ref trend followed the order: Middle > Upper > Yi-Shu-Si > Lower.During the growing season, the magnitude of the change rate was slightly lower than that at the annual time scale, while the spatial distribution was similar to that of the annual time scale.At the seasonal time scale, the magnitude of the change rates of ET ref ranked in sequence as: summer > autumn > spring > winter in most of the HRB, except the lower HRB, while in the lower HRB the magnitude of change rates of ET ref ranked in a sequence as: spring > summer > autumn > winter.In spring, only ET ref in the lower HRB exhibited a significant increasing trend.In summer, ET ref in all sub-regions showed significant decreasing trends.In autumn, ET ref in the middle HRB and Yi-Shu-Si River Basin exhibited significant decreasing trends.However, in winter, no sub-region showed significant change trends.In the whole region, significant decreasing trends were observed at most time scales, except spring and winter, during which the change rate of ET ref was the largest in the annual time scale.
The spatial distribution of seasonal ET ref trends is displayed in Figure 5. Similar spatial distribution (increased from northwest to southeast) was observed at the annual time scale and during the growing season, ranging from −4.87 mm•a −2 to 1.79 mm•a −2 and −4.38 mm•a −2 to 1.28 mm•a −2 , respectively.Among which, most decreasing trends were significant.In spring, the significant downward ET ref trends were detected in the north parts of middle HRB and Yi-Shu-Si River Basin stations.However, ET ref in the majority of stations of the lower HRB exhibited significant upward trends.Similar significant increasing trends could also be detected in some stations in the southern upper and middle HRB.In summer, more than 90% of stations demonstrated the significant decreasing trends of ET ref , except the lower HRB.A general downward trend in ET ref was observed in autumn and winter, with about 52% and 26% of stations exhibiting significant decreasing trends, respectively.decreasing trends.However, in winter, no sub-region showed significant change trends.In the whole region, significant decreasing trends were observed at most time scales, except spring and winter, during which the change rate of ETref was the largest in the annual time scale.
The spatial distribution of seasonal ETref trends is displayed in Figure 5. Similar spatial distribution (increased from northwest to southeast) was observed at the annual time scale and during the growing season, ranging from −4.87 mm•a −2 to 1.79 mm•a −2 and −4.38 mm•a −2 to 1.28 mm•a −2 , respectively.Among which, most decreasing trends were significant.In spring, the significant downward ETref trends were detected in the north parts of middle HRB and Yi-Shu-Si River Basin stations.However, ETref in the majority of stations of the lower HRB exhibited significant upward trends.Similar significant increasing trends could also be detected in some stations in the southern upper and middle HRB.In summer, more than 90% of stations demonstrated the significant decreasing trends of ETref, except the lower HRB.A general downward trend in ETref was observed in autumn and winter, with about 52% and 26% of stations exhibiting significant decreasing trends, respectively.

Sensitivity Coefficient of ET ref to Climate Factors
To explore the effects of climate factors on ET ref , sensitivity coefficients of air temperature, relative humidity, wind speed, and solar radiation in each time period of the HRB are shown in Table 3.Among all climate factors, only the RH showed negative sensitivity coefficients.Moreover, it declared that the ET ref showed more sensitivity to RH, with SR, TA, and WS followed closely at the annual time scale in the whole HRB and every sub-region.In growing season, ET ref in the upper and lower HRB showed the highest sensitivity to RH, while the whole HRB and other sub-regions showed the highest sensitivity to SR.At seasonal time scales, RH was the most sensitive factor in the whole HRB and all sub-regions except in summer.However, in summer, SR was the most sensitive factor in the whole HRB and most sub-regions, except the lower HRB.Meanwhile, RH was still the most sensitive factor in the lower HRB.In addition, the sensitivities of ET ref to TA and SR were higher in the summer and during the growing season, but lower in the winter.The sensitivities of ET ref to RH in autumn and winter seasons were much greater than that of other seasons.A detailed spatial distribution of sensitivity coefficients in each time period can be found in Figure 6.
Water 2018, 10, x FOR PEER REVIEW 12 of 23

Sensitivity Coefficient of ETref to Climate Factors
To explore the effects of climate factors on ETref, sensitivity coefficients of air temperature, relative humidity, wind speed, and solar radiation in each time period of the HRB are shown in Table 3.Among all climate factors, only the RH showed negative sensitivity coefficients.Moreover, it declared that the ETref showed more sensitivity to RH, with SR, TA, and WS followed closely at the annual time scale in the whole HRB and every sub-region.In growing season, ETref in the upper and lower HRB showed the highest sensitivity to RH, while the whole HRB and other sub-regions showed the highest sensitivity to SR.At seasonal time scales, RH was the most sensitive factor in the whole HRB and all sub-regions except in summer.However, in summer, SR was the most sensitive factor in the whole HRB and most sub-regions, except the lower HRB.Meanwhile, RH was still the most sensitive factor in the lower HRB.In addition, the sensitivities of ETref to TA and SR were higher in the summer and during the growing season, but lower in the winter.The sensitivities of ETref to RH in autumn and winter seasons were much greater than that of other seasons.A detailed spatial distribution of sensitivity coefficients in each time period can be found in Figure 6.Note: The bold fonts represent the most sensitive factors in each time period.

Quantitative Contributions of Climate Factors to ETref
The contributions from climate factors to the long-term ETref change trends were calculated by Equation (10).As observed from Figure 7, the change trends of the annual and seasonal ETref estimated using Equation (10) (C_ETref) were well fitted with ETref trends estimated by the Theil-Sen's estimator (T_ETref).The R-squared value was larger than 0.94 in each time period and the slopes of the linear fit approximated to 1.The specific ETref change trends calculated from two approaches in all sub-regions in all time scales can be found in Table 4.The well-fitted relationship between C_ETref and T_ETref suggested that it was reasonable and reliable to adopt the four climate factors in Equation (10) to estimate its contributions to the change trends in ETref.Based on the accuracy and applicability of the modified differential equation method, as shown from Table 4, TA and RH had positive contributions to ET ref in the whole HRB and all sub-regions,

Temporal Change Trends of ET ref and Climate Factors in HRB
Global warming has become an accepted fact in the present era [31].The warming rate in the global mean surface temperature from 1951 to 2012 was 0.012 • C•a −1 .However, the significant increasing rate in TA (0.021 • C•a −1 ) in the entire HRB was two-fold larger than the global average [7], which was echoed by Sun et al. [62].An interesting phenomenon was observed in this study that the TA in the upper and middle HRB demonstrated a slow and non-significant decreasing trend in summer.
Recently, the most significant decreasing rate of WS (−0.021 m•s −1 •a −1 ) was detected in each time scale in the HRB.Similar results were also found elsewhere [63].The decreasing WS rate agreed well with result of Guo et al. [64].However, the magnitude of decreasing rate of WS in the HRB was much greater than those in the other river basins [1,21] and different areas in China [65,66], and also in East Asia [67].The recent declining WS trend is mostly a result of large-scale atmospheric circulation patterns on account of climate warming [68], weakened Siberian High and East Asia Monsoon [69], descending temperature and pressure gradient [69][70][71], and intensified Asian zonal circulation patterns [72].In addition, increased land surface roughness caused by the increased vegetation cover can give a reasonable explanation for the decreasing WS trend.Shan et al. [33] also found this was the reason for declining of WS rate, especially in Northwest China, where Chinese forestry and ecological projects have been implemented in large-scale forestation and reforestation.Furthermore, accelerated urbanization caused by population expansion and economic growth may be reasons for decreased WS, particularly in the regions near weather stations.Additionally, the decreasing WS will reduce the migration of water molecules, reducing the impact of aerodynamic force on ET ref , therefore, lead to a downward trend in ET ref [73].However, the significant downward trend in the WS at each station of the HRB is still a debatable matter and needs further study.
Solar radiation (SR) as one crucial energy source, has found to be in a decreasing trend [1,21,43] which is similar to our results in the HRB.The radiation observation data showed that the SR decreased before 1980s [74,75], then gradually increased after 1990s, this phenomenon was named "from dimming to brightening" by Wild et al. [76].As seen from Figure 3d, this phenomenon is not obvious in the HRB, and SR exhibited a significant downward trend at the rate of −0.03 MJ•m −2 •day −1 •a −1 in the annual time scale in the whole HRB.Possible reasons for this decreasing tendency are (1) the increased regional cloud cover [77]; (2) the increased aerosol loading due to the pollutants from anthropogenic emission [78,79]; (3) increased atmospheric moisture and the aggravated air-pollution caused by the energy depletion; and (4) the significant decreasing WS rate which restrained the spread of aerosols and pollutants [80].
In this study, we have drawn a conclusion that the TA increased significantly with a rate of 0.021 • C•a −1 in last 54 years in the HRB, whereas, the ET ref exhibited a significant downward trend (−1.894 mm•a −2 ) in the whole region.The significant decreasing ET ref in the HRB is in disagreement with climate warming result [81,82].This differs from the earlier findings by Hulme et al. [83] and Goyal [82].In fact, the two adverse change trends between the TA and ET ref are titled the "evaporation paradox" proposed by Brutsaert and Parlange [13].This paradoxical situation occurred in the HRB, especially in the upper and lower HRB and the Yi-Shu-Si River Basin.This finding is echoed with similar works in major basins of the Yangtze River [1], lower Yellow River [21], and Haihe River [84] in China and in many regions worldwide [5,11,85].However, our results about the "evaporation paradox" slightly differ from the outcomes proposed by Thomas [15] which showed an increasing trend of ET ref in the northeast and southwest region of China.This is because of the absence of a long-term meteorological dataset, a shortage of weather stations, and different study periods and areas in his research; whereas our study data source is a long-time series dataset (1961-2014) from 137 weather stations.McVicar et al. [63] conducted a worldwide research to study the ET ref trend in the last three decades, then the downward tendency was detected in many parts worldwide.However, the lower HRB demonstrated an increasing ET ref , this is also proved by results from Chu et al. [16].

Impact of Climate Factors to ET ref Trends
The outcomes exhibited that the RH was the most sensitive factor of ET ref in HRB, followed by SR, TA, and WS, respectively, which coincided with the conclusions proposed by Gong et al. [86] and Li et al. [31] in the Yangtze River Basin and also the Loess Plateau of Northern Shaanxi, China.However, in summer and during the growing season, the ET ref showed the most sensitivity to SR, with TA, RH, and WS following, respectively.Similar findings could also be found in Spain [87] and the United States [88].Additionally, the RH was the most sensitive factor in ET ref in each time period of the lower HRB.All these discoveries are in accordance with the findings from Chu et al. [16] in Jiangsu Province, Eastern China.In fact, the meteorological factors influence ET ref differently in different study areas.Huo et al. [14] indicated that the WS had the strongest sensitivity in arid Northwest China.Ye et al. [32] found SD was the dominating factor, with RH, WS, and TA following in the Poyang Lake catchment.Gao et al. [30] studied the sensitivities of ET ref to the climate factors in the crop growth period in West Liao River Basin, then found SR was the most sensitive factor, with RH following.Moreover, Gao et al. [35] also have demonstrated SR was the most dominant factor for ET ref in the same study area at the annual time scale, but varied at the seasonal time scale.The maximum temperature (T max ) was the dominant climate factor for decreasing ET ref in the Jing River Basin [8].In other countries, Estévez et al. [87] proposed that ET ref was most sensitive to TA, followed by RH, SR, and WS in Spain.Irmak et al. [88] found ET ref was most sensitive to VPD (Vapor pressure deficit) in five locations of the United States, followed by WS in semiarid regions and SR in humid locations during the summer months.Thus, we concluded that diverse climatic conditions, different latitudes and elevations, and local environmental conditions are responsible for the various sensitivity analysis results of meteorological factors to ET ref changes in various areas of the world.
In our research, we also found that WS was the dominant contributing climatic factor in the HRB.This result was echoed with results of Wang et al. [37], Zhang et al. [89], and McVicar et al. [63].As shown above, RH was the most sensitive factor for ET ref in the whole HRB, while WS was the least sensitive factor to the ET ref in the same regions.This phenomenon may be due the change rate of WS was much greater than that of RH.Ultimately, WS exhibited a greater contribution to ET ref trend.Thus, the higher contribution of RH in the lower HRB can be explained reasonably by the significant decreasing trend of RH in each time period shown in Table 1.However, the SR contributed more than WS in the growing season and summer, followed by TA and RH, respectively.This is mainly the result of the high values of SR concentrated in the growing season and summer (Figure 2), and its significantly decreasing trend during these periods (Table 1).Though the monthly change trend of TA was similar to SR, TA exhibited a significant increasing trend, except in summer and during the growing seasons, and these might be the main reason that SR exhibited a higher contribution to ET ref than TA in the same period.However, the non-significance of the TA trend in summer and during the growing season was still unclear and needed to be further investigated.Similar to the sensitivities, the different meteorological factors contributions to ET ref trends were discovered from different areas and in different climatic conditions.Fan et al. [29] reported that the ET ref had different responses to climate factors across various climatic zones in China.Wang et al. [42] declared the decreasing net radiation was responsible for the ET ref decrease in the three-river source area in China.The significantly increasing TA was the main reason for decreasing ET ref in the entire Yellow River Basin [21].Thus, effects of climate factors on ET ref need to be further analyzed, especially in areas with varied climatic conditions, such as the HRB.Even though outcomes of this research indicated that the sensitivity and quantitative contribution analysis can delineate the importance of model parameters to output findings, which allows using it for making an accurate prediction, the methods of sensitivity and quantitative contribution analysis can be extended to other regions.

Uncertainties
We have taken into consideration four main climate factors from the FAO-56 PM model and analyze their effects on ET ref trends.Although we have developed the differential equation method by using RH instead of e a in this study to reduce the estimation error, some errors may still exist due these four meteorological parameters influencing each other to some extent.Meanwhile, the systematic and random errors, inherent in meteorological measurements, could also affect the values of climate factors for estimating ET ref [87].At the same time, the influence of anthropogenic activities has acted as a primary part in recent years (e.g., Peterson et al. [81], Wang et al. [37]) and more uncertain factors are rising.Furthermore, attention should be drawn to some of the uncertainties that exist in the partial derivatives of the ET ref , such as relative and fluctuating values during the long-term study period.The mean values of the partial derivatives of each climatic factor were used to calculate the contributions to ET ref in our study which might generate more uncertainties [43].In addition, on account of the lack of observed solar radiation data, the empirical parameters (a s , b s ) of the Angstrom equation listed in Equation (3) might have changed during the past 54 years' time series due to the changed aerosol optical thickness in this area [90].This is a crucial indication and may also constrain the measured effects of the changes in SR on ET ref , ultimately leading to the underestimation of results by adopting this approximate value for calculating solar radiation from sunshine duration [16].

Conclusions
We found the variation of ET ref could not generate obvious spatial diversity, with the higher values concentrated in the north of middle HRB and Yi-Shu-Si River Basin.ET ref showed a significant downward trend (p < 0.001) in the upper and middle HRB, and also the Yi-Shu-Si River Basin, particularly in summer and during the growing season, while ET ref exhibited a slight increasing trend in the lower HRB, especially significantly in spring.All these phenomena could be explained reasonably by the variations in climate factors.TA had a significant upward trend in most time scales except in summer; RH fluctuated obviously and showed a significant decreasing trend in the lower HRB in each time scale; WS exhibited a significant decreasing trend in all sub-regions in each time scale; and SR had a significant downward trend in most time scales, except in spring.According to the sensitivities of ET ref to climate factors in the HRB, the negative RH showed as the most sensitive factor in most sub-regions, with the positive SR, TA, and WS following, respectively.However, the ET ref was most sensitive to SR in the growing season and summer in most sub-regions except the lower HRB.Based on the developed differential equation method, WS was the most contributing factor, followed by SR, TA, and RH in most time scales, while the SR contributed most to the ET ref trends in growing season and summer in most sub-regions, except the lower HRB.Generally, the negative contribution values of the WS and SR would offset the positive ones of TA and RH, which lead to an overall significant decreasing trend in ET ref in relevant sub-regions and time scales.However, in the lower HRB, the RH contributed most to ET ref , with WS, TA, and SR following, which could be responsible for the slightly increasing trend in ET ref here.
This research deepens insights into the change patterns of regional climate and ET ref , which will offer a scientific theoretical foundation for future investigation on the actual evapotranspiration (ET) of this basin.This research also plays a pivotal role in analogical agricultural crop production and irrigation management.The estimation of ET by using the physical model, as well as remote sensing and land use data, in recent years is an important study, which will delineate the spatiotemporal patterns of ET and its controlling factors in the HRB in the context of rapid urbanization and climatic change.It deserves further investigation and can be a promising extension of this present research work.

Figure 1 .
Figure 1.The location of the Huai River Basin (HRB) in China and its spatial distributions of the elevation, meteorological stations and the river system, four sub-regions represent ①: Upper HRB, ②: Middle HRB; ③: Yi-Shu-Si River Basin; and ④: Lower HRB, respectively.

Figure 1 .
Figure 1.The location of the Huai River Basin (HRB) in China and its spatial distributions of the elevation, meteorological stations and the river system, four sub-regions represent 1 : Upper HRB, 2 : Middle HRB; 3 : Yi-Shu-Si River Basin; and 4 : Lower HRB, respectively.
Water 2018, 10, x FOR PEER REVIEW 7 of 23 obviously, with a non-significant declining trend in the whole HRB.The WS and SR exhibited more significant declining trends in comparison with TA and RH.

Figure 2 .
Figure 2. Box plots show the monthly variations of (a) TA; (b) RH; (c) WS; and (d) SR in the HRB.Figure 2. Box plots show the monthly variations of (a) TA; (b) RH; (c) WS; and (d) SR in the HRB.

Figure 2 .
Figure 2. Box plots show the monthly variations of (a) TA; (b) RH; (c) WS; and (d) SR in the HRB.Figure 2. Box plots show the monthly variations of (a) TA; (b) RH; (c) WS; and (d) SR in the HRB.

Figure 2 .
Figure 2. Box plots show the monthly variations of (a) TA; (b) RH; (c) WS; and (d) SR in the HRB.

Figure 5 .
Figure 5. Spatial distribution maps of (a) annual; (b) growing season; (c) spring; (d) summer; (e) autumn; and (f) winter ETref trends during 1961-2014 in the HRB.The upward and downward triangles represent increasing and decreasing ETref trends The solid black triangle represents a significance level of 0.05.

Figure 5 .
Figure 5. Spatial distribution maps of (a) annual; (b) growing season; (c) spring; (d) summer; (e) autumn; and (f) winter ET ref trends during 1961-2014 in the HRB.The upward and downward triangles represent increasing and decreasing ET ref trends The solid black triangle represents a significance level of 0.05.
The contributions from climate factors to the long-term ET ref change trends were calculated by Equation(10).As observed from Figure7, the change trends of the annual and seasonal ET ref estimated using Equation (10) (C_ET ref ) were well fitted with ET ref trends estimated by the Theil-Sen's estimator (T_ET ref ).The R-squared value was larger than 0.94 in each time period and the slopes of the linear fit approximated to 1.The specific ET ref change trends calculated from two approaches in all sub-regions in all time scales can be found in Table 4.The well-fitted relationship between C_ET ref and T_ET ref suggested that it was reasonable and reliable to adopt the four climate factors in Equation (10) to estimate its contributions to the change trends in ET ref .Water 2018, 10, x FOR PEER REVIEW 13 of 23

Figure 7 .
Figure 7. Relationship between the calculated C_ET ref and the estimated ET ref trend (T_ET ref ) by Theil-Sen's estimator.

Figure 8
Figure 8 exhibits the seasonal and annual contributions of four climate factors to long-term ET ref change trends in each meteorological station in the study period.WS was the dominate driving force to the changes in ET ref at the annual time scale and spring, autumn, and winter seasons, which occupied 93, 94, 109, and 122 stations among a total of 137 stations in the whole HRB, respectively.On the contrary, in the same time scales, TA, RH, and SR were the dominate driving force to the changes in ET ref , occupied 6, 6, 32 stations; 21, 21, 1 stations; 3, 8, 17 stations; and 10, 1, 4 stations among a total of 137 stations in the HRB, respectively.In the growing season and summer, SR contributed most to the change trends of ET ref , especially in summer, which played a dominant role and occupied about 98.5% of stations of the entire HRB.Meanwhile, the TA, RH, WS, and SR acted as dominant parts in the changes of ET ref at 1, 9, 45, and 82 stations during the growing season, respectively.

Figure 8
Figure 8 exhibits the seasonal and annual contributions of four climate factors to long-term ETref change trends in each meteorological station in the study period.WS was the dominate driving force to the changes in ETref at the annual time scale and spring, autumn, and winter seasons, which occupied 93, 94, 109, and 122 stations among a total of 137 stations in the whole HRB, respectively.On the contrary, in the same time scales, TA, RH, and SR were the dominate driving force to the changes in ETref, occupied 6, 6, 32 stations; 21, 21, 1 stations; 3, 8, 17 stations; and 10, 1, 4 stations among a total of 137 stations in the HRB, respectively.In the growing season and summer, SR contributed most to the change trends of ETref, especially in summer, which played a dominant role and occupied about 98.5% of stations of the entire HRB.Meanwhile, the TA, RH, WS, and SR acted as dominant parts in the changes of ETref at 1, 9, 45, and 82 stations during the growing season, respectively.

Figure 8 .
Figure 8. Spatial distribution patterns of contributions of climate factors to ETref in (a) annual; (b) growing season; (c) spring; (d) summer; (e) autumn; and (f) winter during 1961-2014 in the HRB.

Figure 8 .
Figure 8. Spatial distribution patterns of contributions of climate factors to ET ref in (a) annual; (b) growing season; (c) spring; (d) summer; (e) autumn; and (f) winter during 1961-2014 in the HRB.

Figure A1 .
Figure A1.(a) The fitting precision of the SR_HS (SR calculated with Hargreaves-Samani method), (b) compared with the SR_SD (SR calculated with SD), and (c) the fitting precision of ET ref _HS (ET ref calculated with Hargreaves-Samani method) compared with (d) the ET ref _SD (ET ref calculated with SD).
[46]to the lack of measured solar radiation (R s ) data (MJ•m −2 •day −1 ), R s is often calculated through sunshine hours, combed with the Angstrom formula[46]

Table 1 .
Seasonal and annual trend analyses of climatic factors in each region of the HRB during 1961-2014.

Table 2 .
Seasonal and annual trends of ET ref (mm•a −2 ) during 1961-2014 in each region of the HRB.
Note: *, **, and *** denote the significance levels of 0.05, 0.01, and 0.001, respectively.β refers to the estimated slope trend of ET ref , β > 0 and β < 0 signify an upward and a downward trend, respectively.Z is the Mann-Kendall test statistic.

Table 3 .
Annual and seasonal sensitivity coefficients of ET ref to the main climate factors in each region of the HRB during 1961-2014.
Note: The bold fonts represent the most sensitive factors in each time period.3.4.Quantitative Contributions of Climate Factors to ET ref

Table 3 .
Annual and seasonal sensitivity coefficients of ETref to the main climate factors in each region of the HRB during 1961-2014.
Note: The bold fonts represent the most contributing factor in each time period.