Meteorological Factors Affecting Pan Evaporation in the Haihe River Basin , China

Pan evaporation (Epan) is an important indicator of regional evaporation intensity and degree of drought. However, although more evaporation is expected under rising temperatures, the reverse trend has been observed in many parts of the world, known as the “pan evaporation paradox”. In this paper, the Haihe River Basin (HRB) is divided into six sub-regions using the Canopy and k-means (The process for partitioning an N-dimensional population into k sets on the basis of a sample is called “k-means”) to cluster 44 meteorological stations in the area. The interannual and seasonal trends and the significance of eight meteorological indicators, including average temperature, maximum temperature, minimum temperature, precipitation, relative humidity, sunshine duration, wind speed, and Epan, were analyzed for 1961 to 2010 using the trend-free pre-whitening Mann-Kendall (TFPW-MK) test. Then, the correlation between meteorological elements and Epan was analyzed using the Spearman correlation coefficient. Results show that the average temperature, maximum temperature, and minimum temperature of the HRB increased, while precipitation, relative humidity, sunshine duration, wind speed and Epan exhibited a downward trend. The minimum temperature rose 2 and 1.5 times faster than the maximum temperature and average temperature, respectively. A significant reduction in sunshine duration was found to be the primary factor in the Epan decrease, while declining wind speed was the secondary factor.


Introduction
Global warming has become an indisputable fact [1].Temperature records indicate that the earth has warmed by approximately 0.6 • C during the 20th century [2].This increase in global temperature has significantly impacted the natural environment, ecosystem, and social economy [3], and has led to a series of changes in hydrological factors, such as precipitation, evaporation, water infiltration, soil moisture, river runoff, and groundwater flow, all of which affect the global hydrological cycle.This, in turn, causes temporal and spatial redistribution of water resources, and thereby threatens water security, food security, social security, and national security [4,5].As a key component in the hydrological cycle, evapotranspiration is associated with water balance and water exchange, as well as surface energy balance; hence, of all components of the water cycle, evapotranspiration is the factor most directly affected by climate change [6].Therefore, analyzing the climate sensitivity of evapotranspiration has important theoretical and practical implications for understanding the impact of climate change on the hydrological cycle [7].
Evapotranspiration is the process of water transport from the earth's surface to the atmosphere [8].As a core process of the climate system, evapotranspiration closely links the hydrological cycle, energy budget, and carbon cycle [9].Pan evaporation (E pan ) [10] is the most universal and simplest way to measure evapotranspiration, which is often used to indicate the humidity level of a given regional climate [11].Although E pan cannot directly represent the evaporation of the water surface, it has a close correlation with water surface evaporation.Therefore, it has remained an important reference indicator in the assessment of water resources, water resources planning, and the design of irrigation systems, to name a few examples [12].As the global temperature rises, E pan should theoretically gradually increase.However, in reality, only certain regions in the world have an E pan value that is consistent with theoretical expectations, and the majority of the world's regions have been found to have declining E pan values.This phenomenon is called the "pan evaporation paradox" [3].Specifically, countries such as Spain [13], Iran [14], Israel [15], and Brazil [16] have been found to have increasing E pan values, and countries such as the former Soviet Union, the United States [17,18], New Zealand [19], China [20][21][22][23], Thailand [24], India [12], Nigeria [25], and Australia [26,27] have been found to have declining E pan values.Correctly interpreting the overall declining trend of E pan in the context of rising global temperatures and uncovering the main meteorological factors that affect the reduction of E pan is of great importance to accurately predict future hydrological cycles.
Many scholars have studied the temporal and spatial changes of E pan at global and regional scales, as well as the causes of such changes.According to their findings, the causes of E pan reduction can be categorized as follows.(1) An increase in humidity in the surrounding environment of the evaporation pan: Brutsaert and Parlange ascertained that the decrease in E pan value was due to an increase in the volume of evaporation from the land surface, considering the difference between evaporation from the land surface and the evaporation volume observed through the evaporation pan [28].Zuo et al. employed observational data from 62 conventional meteorological stations with solar-radiation observation equipment in China to analyze in detail the relationship between E pan and corresponding environmental factors, as well as the environmental factors' responses to global climate change.The researchers discovered that E pan was most correlated to atmospheric relative humidity [20].(2) Changes in precipitation: Tebakari et al. [24] analyzed the temporal and spatial variation of E pan in Thailand from 1982 to 2000 and concluded that both E pan and precipitation showed a declining trend.This conclusion was inconsistent with findings from the United States, where E pan was found to be decreasing while precipitation was increasing [29].Jaswal et al. utilized evaporation and rainfall data from 1971 to 2000 from 58 stations that were evenly distributed in India to analyze the overall correlation between evaporation and rainfall in a year, as well as their correlation in winter, summer, monsoon season, and post-monsoon season.The results showed that, in southern India, the evaporation trend had a complementary relationship with rainfall during the same period [30].
(3) A decrease in the diurnal temperature range: Peterson et al. compared data from both the United States and the former Soviet Union from 1950 to 1990 and found a steady decline in E pan values in all investigated regions (except Central Asia), as well as a decline in diurnal temperature range.E pan and diurnal temperature range were thus clearly correlated.Therefore, the researchers concluded that the reduction in the diurnal temperature range, caused by an increase in cloud cover, consequently caused the reduction in E pan [17].(4) A reduction in solar radiation: Roderick and Farquhar found that E pan values observed in many parts of the world over the past 50 years showed a clear downward trend and asserted that such a decline was caused by the reduction in overall solar radiation resulting from an increase in cloud cover and aerosol concentrations [3].(5) A reduction in wind speed: Burn and Hesch conducted a trend analysis on the evaporation data of 48 sites in the Canadian Prairies over three analysis periods and concluded that wind speed has a substantial influence on the decreasing trend of evaporation, while vapor-pressure deficit has a significant influence on the increasing trend of evaporation [31].Hoffman et al. studied the changes in E pan , rainfall, wind speed, temperature, and vapor-pressure deficit from 1974 to 2005 taken from 20 climate stations in the Cape Floristic Region (CFR), South Africa, and suggested that the reduction in E pan was likely due to a reduction in wind speed [32].Yang and Yang analyzed the daily E pan , temperature, wind speed, solar radiation, and relative humidity of 54 meteorological stations in China for 1961 to 2001 and concluded that the reduction in E pan in the majority of regions in China is due to a decrease in wind speed [33].(6) The comprehensive impact of meteorological elements: Roderick and Farquhar analyzed data from Australia for 1970 to 2002 and found that E pan values showed a downward trend.The results showed that such a change might be related to a decrease in solar radiation, wind speed, and diurnal temperature range [26].Sheng examined E pan data and other meteorological factors from 468 meteorological stations in China, measured simultaneously from 1957 to 2001, and found that the main influential factors of E pan were solar radiation, diurnal temperature range, and wind speed, while the influence of humidity was the weakest factor [21]. Liu et al. investigated data for 1955 to 2001 taken from 671 sites in China.The results revealed an overall decline in E pan .In addition, diurnal temperature range and wind speed were found to have the greatest correlation with such a decline [22].Based on the aforementioned studies, the causes of the reduction of E pan appear to be very complicated.Owing to the location, climate, atmospheric differences, and even the differences in the length of the data series, the conclusions of these studies are inconsistent.Therefore, identifying the impact of various meteorological variables on E pan trends is critical to quantifying the impact of global warming.
The HRB is located in a region with a warm semi-arid climate and a continental monsoon climate.This area is sensitive to climate change and is a region with a fragile ecological environment.Owing to the area's dense population and rapid economic development, as well as its status as one of China's major wheat producers, the contradiction between water supply and water demand is prominent in the area.Water shortages have become a major factor restricting sustainable economic and social development in the HRB [34].E pan in the HRB generally exhibits a decreasing trend [35][36][37], which is largely consistent with E pan trends in other regions of China [38][39][40][41][42][43][44][45][46][47][48].However, scholars have differing views on the causes of the E pan trend in the HRB.Zheng et al. analyzed the effects of temperature, wind speed, solar radiation, and atmospheric pressure on E pan in the HRB for 1957 to 2001 and concluded that wind speed is the main factor leading to the decrease of E pan in the region [49].Hao et al. selected eight meteorological elements from 34 climate stations for 1958 to 2011 in order to analyze the spatial and temporal variations in the HRB.The results showed that the potential evapotranspiration in the region was negatively correlated with relative humidity and was positively correlated with diurnal temperature range [50].Guo and Ren examined data observed from the evaporation pans of 117 meteorological stations for 1956 and 2000 and analyzed the changes in evaporation in the Huang-Huai-Hai River Basin.The findings showed that the direct climatic cause of the decrease in evaporation may be a reduction in sunshine duration and solar radiation.In addition, a reduction in wind speed and diurnal temperature range may also play an important role [51].In summary, though most of the papers consider that the decrease of wind speed is a main factor causing the E pan declining, but different literatures have different conclusions on the influence of other meteorological factors on the E pan decreasing.So, this study aimed to analyze the trend of changes in E pan in the HRB using data collected by 44 meteorological stations for 1961 to 2010, and to explore the temporal and spatial variation laws of E pan , as well as the main driving forces of declining E pan trend in the region.

Materials and Methods
The HRB is located between 112-120 • E longitude and 35-43 • N latitude, with the Bohai Sea to the east, the Yellow River to the south, the Yunzhong and Taiyue Mountains to the west, and the Inner Mongolia Plateau to the north.The total area of the HRB is 320,600 km 2 , accounting for 3.3% of the total area of the country.The HRB spans eight provinces: Beijing, Tianjin, Hebei, Shanxi, Shandong, Henan, Inner Mongolia, and Liaoning.It is a political and cultural center and an economically developed region of strategic significance in China.The HRB has two major rivers: the Hai River and the Luan River.The Hai River, which is the main water system for the area, consists of the Ji Canal River, Chaobai River, North Canal, Yongding River, Daqing River, Ziya River, and Zhangwei River, as well as plains rivers, such as the Tuhai River and Majia River, each of which enters the sea individually.The Luan River includes itself and the rivers along the eastern coast of Hebei Province.The annual average temperature range of the basin is between 1.5 • C and 14 • C, the annual average relative humidity is between 50% and 70%, and the average annual precipitation is 539 mm (semi-arid climate).The annual average land-surface evaporation is 470 mm, and the water surface evaporation is 1100 mm.The geographical location and topographic distribution of the study area are shown in Figure 1.
annual average relative humidity is between 50% and 70%, and the average annual precipitation is 539 mm (semi-arid climate).The annual average land-surface evaporation is 470 mm, and the water surface evaporation is 1100 mm.The geographical location and topographic distribution of the study area are shown in Figure 1.Meteorological data from the HRB and 55 meteorological stations in the surrounding area, provided by the National Meteorological Center of the China Meteorological Administration, were used in this study, including the daily average temperature, highest, and lowest temperatures, average relative humidity, sunshine duration, wind speed, precipitation, and daily evaporation from an evaporation pan 20 cm in diameter.In terms of missing data, the following rules were respected: when the daily data for five or more days were missing for a specific meteorological element in a month, the data of the entire month were considered missing; when the data for one or more months were missing, the data of the entire year were considered invalid.The time series of the data was from 1961 to 2010, and the length of the time series was 50 years.After excluding the station data that did not satisfy the time series requirements, data for 44 stations were retained.
In order to better analyze the seasonal changes of the elements, the data were divided into spring (March to May), summer (June to August), fall (September to November), and winter (December to February).The annual temperature (average, maximum, and minimum), annual average relative humidity, annual average sunshine duration, and annual average wind speed of each station were calculated based on the mean of the daily data.The annual Epan and precipitation were calculated by summing the daily data.The same methods were applied to obtain the seasonal data for each element.

Canopy and k-means Clustering
In order to explore patterns in the spatial and temporal variations of the HRB climate, as well as geomorphological differences in the region and the climatic characteristics, this study selected eight representative indicators as references to categorize the HRB into several sub-regions.The indicators included geodetic coordinates (X and Y values), elevation, average temperature, precipitation, relative humidity, sunshine duration, and wind speed.The Canopy and k-means clustering technique was adopted.This method performs clustering in two stages.In Stage 1, the Canopy clustering Meteorological data from the HRB and 55 meteorological stations in the surrounding area, provided by the National Meteorological Center of the China Meteorological Administration, were used in this study, including the daily average temperature, highest, and lowest temperatures, average relative humidity, sunshine duration, wind speed, precipitation, and daily evaporation from an evaporation pan 20 cm in diameter.In terms of missing data, the following rules were respected: when the daily data for five or more days were missing for a specific meteorological element in a month, the data of the entire month were considered missing; when the data for one or more months were missing, the data of the entire year were considered invalid.The time series of the data was from 1961 to 2010, and the length of the time series was 50 years.After excluding the station data that did not satisfy the time series requirements, data for 44 stations were retained.
In order to better analyze the seasonal changes of the elements, the data were divided into spring (March to May), summer (June to August), fall (September to November), and winter (December to February).The annual temperature (average, maximum, and minimum), annual average relative humidity, annual average sunshine duration, and annual average wind speed of each station were calculated based on the mean of the daily data.The annual E pan and precipitation were calculated by summing the daily data.The same methods were applied to obtain the seasonal data for each element.

Canopy and k-means Clustering
In order to explore patterns in the spatial and temporal variations of the HRB climate, as well as geomorphological differences in the region and the climatic characteristics, this study selected eight representative indicators as references to categorize the HRB into several sub-regions.The indicators included geodetic coordinates (X and Y values), elevation, average temperature, precipitation, relative humidity, sunshine duration, and wind speed.The Canopy and k-means clustering technique was adopted.This method performs clustering in two stages.In Stage 1, the Canopy clustering algorithm is used to calculate the similarity of the objects and to categorize similar objects in the same subset (canopy).In Stage 2, the k-means clustering algorithm [52] is used to cluster the points in each canopy.Once Stage 1 is complete, the algorithm only needs to accurately cluster the points in each canopy, which greatly reduces the time spent on the accurate calculation of all data points that was performed in a traditional clustering algorithm.In addition, the number of canopies obtained in Stage 1 can be used as the K value in K-means clustering, which may minimize the irrational selection of the K value to a certain degree.The Canopy and k-means clustering technique not only greatly reduces the calculation of distance between points, the result is also more accurate when compared to general clustering methods [53,54].
MATLAB software (MATLAB 9.0, R2016a, MathWorks, Natick, MA, USA) was used to train the algorithm.The clustering results of the stations and their spatial distribution are presented in Figure 2. According to the clustering results, the HRB could be divided into six sub-regions (Figure 3).Detailed information for each sub-region is shown in Table 1.algorithm is used to calculate the similarity of the objects and to categorize similar objects in the same subset (canopy).In Stage 2, the k-means clustering algorithm [52] is used to cluster the points in each canopy.Once Stage 1 is complete, the algorithm only needs to accurately cluster the points in each canopy, which greatly reduces the time spent on the accurate calculation of all data points that was performed in a traditional clustering algorithm.In addition, the number of canopies obtained in Stage 1 can be used as the K value in K-means clustering, which may minimize the irrational selection of the K value to a certain degree.The Canopy and k-means clustering technique not only greatly reduces the calculation of distance between points, the result is also more accurate when compared to general clustering methods [53,54].MATLAB software (MATLAB 9.0, R2016a, MathWorks, Natick, MA, USA) was used to train the algorithm.The clustering results of the stations and their spatial distribution are presented in Figure 2. According to the clustering results, the HRB could be divided into six sub-regions (Figure 3).Detailed information for each sub-region is shown in Table 1.

Trend-Free Pre-Whitening Mann-Kendall Test (TFPW-MK)
The Mann-Kendall (M-K) test [55,56] is a non-parametric statistical method.Compared to parametric statistical methods, the M-K test does not require samples to follow a certain distribution, the results are not subject to interference from a few outliers, and the method is simple and efficient in calculations.Therefore, it is commonly used to detect trends in a series of values.For that reason, the M-K test is suitable for examining the trend of the hydrological variables in this study [57][58][59].Assuming  ,  , ⋯ ,  is a time series, n is the length of the time series; then, the M-K method defines the statistical variable S as follows: where xj and xk are the measured values of years j and k, respectively; and ,  ≤  and  ≠ .
When the number of samples is greater than 10, Z is calculated as follows: where, Z is a normally distributed statistic, and Var(S) is the variance.If the Z value is positive, the data shows an increasing trend; if the Z value is negative, the data shows a decreasing trend [60].

Trend-Free Pre-Whitening Mann-Kendall Test (TFPW-MK)
The Mann-Kendall (M-K) test [55,56] is a non-parametric statistical method.Compared to parametric statistical methods, the M-K test does not require samples to follow a certain distribution, the results are not subject to interference from a few outliers, and the method is simple and efficient in calculations.Therefore, it is commonly used to detect trends in a series of values.For that reason, the M-K test is suitable for examining the trend of the hydrological variables in this study [57][58][59].Assuming X 1 , X 2 , • • • , X n is a time series, n is the length of the time series; then, the M-K method defines the statistical variable S as follows: where x j and x k are the measured values of years j and k, respectively; and k, j ≤ n and k = j.
When the number of samples is greater than 10, Z is calculated as follows: where, Z is a normally distributed statistic, and Var(S) is the variance.If the Z value is positive, the data shows an increasing trend; if the Z value is negative, the data shows a decreasing trend [60].Given the level of significance α, if |Z| ≥ Z 1−α/2 , the null hypothesis is rejected, and the trend of the time series data (increasing or decreasing) is statistically significant at α.The existence of serial correlation increases the probability that the M-K test will detect a significant trend [61,62].The meteorological and hydrological data are mostly skewed and do not follow the same distribution, and there may be autocorrelation.Thus, in this paper, the TFPW method proposed by Yue et al. [62] is used to limit the influence of serial correlation; then, the significance of the time series is assessed by the M-K test.
The slope of a trend is estimated using the TSA [63][64][65][66], and it is estimated as follows: where β is the estimate of the slope of the trend, and X i is the ith observation.The slope determined by the TSA is a robust estimate of the magnitude of a trend.Since the publication of Hirsch et al. [67], the TSA has been popularly employed to identify the slope of trends in hydrological time series [68][69][70].
Step 2. If β = 0, there is no need to continue trend analysis; if β = 0, it is assumed to be linear, and the sample data are detrended as: Step 3. The lag-1 serial correlation coefficient r 1 of Y t is calculated using Equation ( 7), and then the autocorrelation is removed by Equation ( 9).
Step 4. The identified trend T t and the residual Y t are blended by Step 5. Verify the significance of trend of the blended series using the MK test.

Spearman Correlation Coefficient
The Spearman correlation coefficient [71] is a nonparametric test method independent of distribution and can be used as an indicator to measure the relationship between two variables.If there are no repeated values in the data, the Spearman correlation coefficient is +1 or 1 when two variables are monotonously correlated.For a sample of size n, the n raw scores X i , Y i are converted to ranks x i , y i , and the Spearman correlation coefficient ρ can be calculated as [72,73].
where d i = x i − y i is the difference between the two ranks of each observation.The correlation degree between X i , Y i can be used according to the grading standards of ρ shown in Table 2 [74].

Trend and Significance Analysis
This study adopted the M-K test to analyze E pan , temperature (average, maximum, and minimum), precipitation, relative humidity, sunshine duration, and wind speed of the HRB.The TFPW method was used to eliminate the trends and autocorrelation of meteorological sequence data before the M-K Test was applied.If the value of Z was positive, the data exhibited an upward trend, while if the value of Z was negative, the data exhibited a downward trend.The threshold of the significance level was defined as α = 0.05.If the change trend of a given meteorological variable was found to be significant at this level, then |Z| > Z α 2 = 1.96 [75].The results of significance test of the interannual variations are shown in Figure 4, and that of the seasonal variations are shown in Table 3.
For TFPW-MK test detrending, the TSA method is used to calculate the magnitude of the trend of the meteorological variables.The rates of the meteorological elements of each sub-region for 1961 to 2010 are shown in Table 4.
Water 2019, 10, x FOR PEER REVIEW 8 of 18 where is the difference between the two ranks of each observation.The correlation degree between i X , i Y can be used according to the grading standards of ρ shown in Table 2 [74].

Trend and significance analysis
This study adopted the M-K test to analyze Epan, temperature (average, maximum, and minimum), precipitation, relative humidity, sunshine duration, and wind speed of the HRB.The TFPW method was used to eliminate the trends and autocorrelation of meteorological sequence data before the M-K Test was applied.If the value of Z was positive, the data exhibited an upward trend, while if the value of Z was negative, the data exhibited a downward trend.The threshold of the significance level was defined as 0.05 α = . If the change trend of a given meteorological variable was found to be significant at this level, then  3.
For TFPW-MK test detrending, the TSA method is used to calculate the magnitude of the trend of the meteorological variables.The rates of the meteorological elements of each sub-region for 1961 to 2010 are shown in Table 4.As shown in Figure 4, the interannual average temperatures (Tmean) in sub-regions I to VI presented a significant upward trend.With the exception of sub-region V, the average maximum As shown in Figure 4, the interannual average temperatures (T mean ) in sub-regions I to VI presented a significant upward trend.With the exception of sub-region V, the average maximum temperatures (T max ) of all sub-regions also increased significantly.The significance level of the trend of the average minimum temperatures (T min ) was consistent with that of the T mean .The T mean in spring and winter in sub-regions I to VI increased significantly; in summer, it increased significantly only in sub-regions II and III.In fall, it increased significantly in the majority of sub-regions (except for sub-region I).Only the significance levels for T max in winter of sub-regions I to VI were consistent with those of T max .The T min of sub-regions I to VI in all four seasons significantly increased.However, it can be seen from Table 2 that the T min rose more rapidly, followed by the T mean , with the T max rising the most slowly.
The interannual variations of average precipitation (P mean ) showed a downward trend of the six sub-regions: the decline rates were 22.50, 5.02, 4.07, 23.49, 5.35 and 19.71 mm/10a, respectively.However, the trend is not significant.The P mean in spring showed an upward trend; only the trends of sub-regions II and VI were significant.In summer, it was found to be declining, except for sub-region V; the declining trends of sub-regions I, II, and VI were significant.The P mean in fall showed an upward trend, except for sub-regions IV and V; The P mean in winter showed an upward trend except for sub-region VI.The trends were not statistically significant in fall and in winter.As the decline of precipitation in summer offsets the increase of precipitation in spring, the interannual P mean of the region showed a general downward trend.
The interannual variations of average relative humidity (RH mean ) of all sub-regions showed a downward trend, except for sub-region I; however only the trends of sub-regions IV and VI were statistically significant.The change rates of sub-regions I to VI were 0.20%/10a, −0.35%/10a, −0.30%/10a, −0.71%/10a, −0.56%/10a, and −0.78%/10a, respectively.As the increase in the RH mean in fall and winter was larger than the sum of the reduction of the RH mean in spring and summer, the RH mean in sub-region I was found to be increasing.The RH mean in spring, summer, and fall of sub-regions I to VI was found to be declining, except for sub-region I in fall and sub-region V in summer; however, the trends of the majority of sub-regions were not significant.The changes in the RH mean of winter were not consistent across sub-regions I to VI, and none of the analyzed trends were significant.
The interannual variations of the average sunshine duration (SD mean ) in sub-regions I to VI showed a significant downward trend; the decline rates were 0.27, 0.13, 0.17, 0.29, 0.31, and 0.32 h/10a, respectively.In addition, the significance values for the changes in the SD mean per season in sub-regions I to VI were consistent with those of the interannual variation.
The interannual average wind speed (U mean ) of sub-regions I to VI had significantly decreased; the rates were 0.16, 0.15, 0.28, 0.05, 0.17, and 0.16 m/(s•10a), respectively.The seasonal variation of the U mean was showed a downward trend in all sub-regions, except for sub-region VI in summer.The U mean in summer and fall in sub-region IV was not significant, while the declining trends of the U mean per season of the other sub-regions were found to be significant.
The interannual E pan of sub-regions I to VI was found to be decreasing over the research period at speeds of 55.2, 11.7, 24.4,30.7, 82.6, and 65.7 mm/10a, respectively.However, the trends of E pan for all sub-regions were significant except for sub-regions II and III.The E pan of sub-region II presented an upward trend in summer, fall, and winter, while the E pan of other sub-regions showed a downward trend in all seasons.Apart from sub-region IV, the E pan of other regions was significantly decreased in spring.The significance test results in summer aligned with those of the interannual variations.The changes of E pan in fall and winter in the majority of sub-regions were not statistically significant.

Sensitivity Analysis
Changes in meteorological variables led to temporal and spatial fluctuations in E pan , and their roles appeared to be different in different sub-regions.The influence of each meteorological element on changes in E pan depended on two factors: the sensitivity of the meteorological elements toward E pan and the change trend and corresponding significance level of the meteorological element.Therefore, it was necessary to analyze the meteorological elements that drive changes in E pan for each sub-region, to qualitatively evaluate the contribution of each meteorological element on the changes of E pan .The Spearman correlation coefficient was applied to qualitatively analyze the effects of T max , T mean , T min , P mean , RH mean , SD mean , and U mean on E pan .The results are shown in Table 5.
The sensitivity factors that caused changes in E pan in the HRB were different.It can be seen from Table 5 that there were significant correlations between E pan and T min , P mean , RH mean , SD mean , and U mean in sub-region I.Although the correlations between E pan and P mean and RH mean were significant, the trends of P mean and RH mean were not significant.Therefore, the meteorological elements that affected the interannual changes in E pan of sub-region I were T min , SD mean , and U mean .E pan was found to have a negative correlation with T min , which does not align with what is commonly expected.E pan and SD mean and U mean were positively correlated.The significant decline in SD mean and U mean directly led to the significant decline of E pan .Based on the above analysis, the significant decrease in E pan in sub-region I was mainly due to the significant reduction in the SD mean and U mean .The Spearman correlation coefficient between E pan and SD mean was 0.75, indicating that the correlation is strong (see Table 2).The Spearman correlation coefficient between E pan and U mean was 0.46, indicating a moderate correlation.Thus, the primary factor affecting E pan decline in sub-region I was SD mean , followed by U mean .According to the same analysis method, the factors responsible for E pan decline in each sub-region were also analyzed.The primary factor affecting the decline of E pan in sub-regions I, III, IV, V, and VI was SD mean , followed by U mean , while the reasons for the changes in E pan in sub-regions II was the significant reduction in the SD mean .
In addition, the contributing factors to E pan in sub-regions I to VI for each season were also analyzed.In spring, the decrease in E pan in sub-regions I to VI was primarily caused by a significant decrease in SD mean , followed by a decline in U mean .In summer, significant decreases in SD mean and U mean were the primary and secondary driving factors of E pan in sub-regions I, V, and VI, respectively; the decrease in E pan for sub-regions III and IV was mainly due to the significant decrease in SD mean .Moreover, the change trend of E pan in sub-region II was found to be the opposite to that for other sub-regions, which could be due to a significant increase in temperature, as well as a significant drop in P mean and RH mean .In fall, the decline in E pan in sub-regions I, III, V, and VI was attributed to a significant reduction in SD mean and U mean , while the decline in E pan in sub-region IV was caused by a significant reduction in SD mean .In addition, the change trend of E pan in sub-region II was the opposite of that for other sub-regions, which was mainly caused by a significant increase in T max and T mean .In winter, apart from sub-region II, where a significant increase in temperature resulted in an upward trend in E pan , the reduction of E pan in all sub-regions was caused by a significant reduction in SD mean and U mean .
In summary, the primary factor responsible for the decline in E pan in the HRB was a reduction in sunshine duration, followed by a reduction in wind speed.The factors responsible for the reduction in E pan in each sub-region were consistent with the overall reduction in E pan of the HRB.However, the correlation between U mean and E pan was not significant in sub-region II, where the only influential factor was the decrease of sunshine duration.

Discussion
As noted above, when global temperature increases, the overall temperature of the HRB increases.However, there are differences in the spatial and temporal distribution.From 1961 to 2010, the lowest temperature increased approximately two times faster than the highest temperature, and approximately 1.5 times faster than the average temperature.This result is consistent with the results from an analysis of variations in annual temperatures by Zheng et al. [49] in the HRB for 1957 to 2001, but it is different from the results obtained by Salinger and Griffiths [76], which indicated that the lowest temperature rose approximately three times faster than the highest temperature globally from 1951 to 1998.Meanwhile, the average temperature, highest temperature, and lowest temperature in spring and winter increased much more quickly than the corresponding values in summer and autumn, for both the HRB and every sub-region.The temperatures in winter increased most significantly, three times faster than those in summer, and two times faster than those in autumn.Additionally, temperatures decreased to some extent in some areas of the HRB.For example, unlike the highest temperature in other sub-regions, which increased, the highest temperature in area V decreased.
The "evaporation paradox" also exists in the HRB.With respect to both the whole area and the sub-regions in the HRB, except for the slight increase in E pan in sub-region II in autumn and winter, E pan decreased, and this decline mostly occurred in spring and summer.This result agrees with the conclusions of Zheng et al. [49] and Liu et al. [77].However, it contradicts with the findings in Liu et.al. [78] that from 1992 to 2007, E pan significantly increased in North China, including the HRB.The correlation between the general decrease in E pan and the temperature variation was weak, suggesting that the increase in temperature was not directly related to the decrease in E pan .The main factors responsible for the decline in E pan in the HRB included decreases in sunshine duration, which was the main factor, and in wind speed.This result differed from the conclusions of Zheng et al. [49], which stated that the decrease in wind speed is the main factor responsible for the decrease in E pan .In addition, this result does not agree with the conclusions from Liu et al. [79] that in the semi-humid/semi-arid region of China (including the HRB), decreases in diurnal temperature range, sunshine duration and wind speed were found to be the main factors contributing to the pan evaporation declines.Liu et al. [78] concluded that wind speed and solar radiation are the main factors that led to the decline in pan evaporation in North China, which differs from our findings.However, it is generally accepted that wind speed is one of the main driving factors of the decrease in E pan in the HRB.This statement is consistent with the conclusion that wind speed is one of the main factors driving the decrease in E pan in areas such as the Canadian Prairies [31], the Cape Floristic Region in South Africa [32], and Australia [26].
The factors attributed to the decrease in sunshine duration and wind speed also vary among studies [36,[80][81][82].For sunshine duration, multiple studies have concluded that this decrease may be related to the increase in aerosols and other air pollutants [3,83].In other studies, it is argued that the decrease may be related to the increase in cloud cover.In addition, several studies have reported a correlation between a decrease in sunshine duration and urbanization [84].Recently, Wei Pan [85] reported that the number of haze days significantly increased in North China (including in the HRB).Therefore, a decrease in sunshine duration may be related to the increase in haze days in the HRB.Regarding the decrease in wind speed, conclusions of various areas also differ; however, the main consensus is that the decrease in wind speed may be related to variations in global circulation [86,87], as well as the increase in surface roughness caused by afforestation and urbanization near the observation sites [88].Based on the present studies, it is difficult to determine the reasons for the decrease in wind speed, and further studies are required.

Conclusions
In this study, the Canopy and k-means clustering method was employed to categorize the HRB into six sub-regions.Then, 44 out of the 55 meteorological stations in the surrounding area that had relatively complete data were selected, and the trends and significance of the interannual and seasonal variations of the pan evaporation, temperature, precipitation, relative humidity, sunshine duration, and wind speed for 1961 to 2010 were analyzed using TFPW-MK.Based on this analysis, the sensitivities of the average, maximum, and minimum temperatures, and precipitation, relative humidity, sunshine duration, and average wind speed to E pan were qualitatively analyzed using the Spearman correlation coefficient.In the whole basin, the primary cause of declining E pan was a significant reduction in sunshine duration, followed by a significant reduction in wind speed.In sub-regions, E pan showed a downward trend; however, the influential factors on E pan reduction per sub-region were slightly different from those of the entire region.Except for sub-region II, which was only affected by sunshine duration, reductions in E pan in other sub-regions were due to the joint influence of decreasing sunshine duration and wind speed.
In this paper, only a qualitative analysis was performed on the reduction of sunshine duration and wind speed in the HRB, and an explanation for this reduction is still lacking, and thus further research is needed.

Figure 1 .
Figure 1.Location and topography of the HRB.

Figure 1 .
Figure 1.Location and topography of the HRB.

Figure 2 .
Figure 2. The Canopy and k-means clustering results and spatial distribution of the stations in the HRB.

Figure 2 .
Figure 2. The Canopy and k-means clustering results and spatial distribution of the stations in the HRB.

Figure 3 .
Figure 3.The sub-region categorization of the HRB.

Figure 3 .
Figure 3.The sub-region categorization of the HRB.
. The results of significance test of the interannual variations are shown in Figure 4, and that of the seasonal variations are shown in Table

Figure 4 .
Figure 4. Results of significance test on the interannual variations of the meteorological elements.

Figure 4 .
Figure 4. Results of significance test on the interannual variations of the meteorological elements.

Author Contributions:
This paper is a joint effort by several authors.S.W. conceived and designed the paper's structure; B.L. performed the Canopy-Kmeans clustering; Z.Y. compiled the code for the Mann-Kendall Test; Z.Y. and D.M. analyzed the data; H.L. and S.L. drew the figures; and Z.Y. wrote the paper.Funding: This work was supported by the National Natural Science Foundation of China (Grant No. 71573274; grant No. 51879066), the National Water Pollution Control and Management Science and Technology Major Project of China (2014ZX07203-008), the Science and Technology Project of Hebei Province (Grant No. 15227005D), and the Science and Technology Research and Development Program of Handan (1723209055-4).

Table 1 .
Basic information for the sub-regions of the HRB.

Table 1 .
Basic information for the sub-regions of the HRB.

Table 3 .
Z-Value of the Mann-Kendall test on the meteorological variables (α = 0.05).
* Trends statistically significant at the 95% confidence level.

Table 4 .
Climate tendency rates of the meteorological elements per sub-region from 1961 to 2010.

Table 5 .
The Spearman correlation coefficients between Meteorological Elements and E pan per Sub-Region.
* Trends statistically significant at the 95% confidence level.