Characteristics of Potential Evapotranspiration Changes and Its Climatic Causes in Heilongjiang Province from 1960 to 2019

: Climate change refers to the statistically signiﬁcant changes in the mean and dispersion values of meteorological factors. Characterizing potential evapotranspiration ( ET 0 ) and its climatic causes will contribute to the estimation of the atmospheric water cycle under climate change. In this study, based on daily meteorological data from 26 meteorological stations in Heilongjiang Province from 1960 to 2019, ET 0 was calculated by the Penman–Monteith formula, linear regression method and the Mann–Kendall trend test were used to reveal the seasonal and inter-annual changing trend of ET 0 . The sensitivity-contribution rate method was used to clarify the climatic factors affecting ET 0 . The results showed that: (1) From 1960 to 2019, the maximum temperature ( T max ), minimum temperature ( T min ) and average temperature ( T mean ) showed an increasing trend, with climate tendency rate of 0.22 ◦ C per decade (10a), 0.49 ◦ C/(10a), 0.36 ◦ C/(10a), respectively. The relative humidity ( RH ), wind speed ( U ) and net radiation ( Rn ) showed a decreasing trend, with a climate tendency rate of − 0.42%/(10a), − 0.18 m/s/(10a), − 0.08 MJ/m 2 /(10a), respectively. (2) ET 0 showed a decreasing trend on seasonal and inter-annual scales. Inter-annually, the average climate tendency rate of ET 0 was − 8.69 mm/(10a). seasonally, the lowest climate tendency rate was − 6.33 mm/(10a) in spring. (3) ET 0 was negatively sensitive to T min , and RH , while positively sensitive to T max , T mean U and Rn , its sensitivity coefﬁcient of U was the highest, which was 1.22. (4) The contribution rate of U to ET 0 was the highest on an inter-annual scale as well as in spring and autumn, which were − 8.96%, − 9.79% and − 13.14%, respectively, and the highest contribution rate to ET 0 were Rn and T min in summer and winter, whose contribution rates were − 4.37% and − 11.46%, respectively. This study provides an understanding on the response of evapotranspiration to climatic change and further provides support on the optimal allocation of regional water resource and agricultural water management under climate change.


Introduction
Over the past 100 years , the global average surface temperature has increased by 0.74 • C and is expected to increase by at least 4 • C in 2100 if carbon dioxide emissions are not reduced [1].Climate change will inevitably lead to changes in the global hydrological cycle [2].Apart from precipitation and runoff, actual evapotranspiration (ET) is the most important climatic factor and a main component in the hydrological cycle [3].Evapotranspiration is affected by a variety of meteorological factors, such as temperature, relative humidity (RH), wind speed (U), and net radiation (Rn).Characteristics of ET changes and its climatic causes under climate change will reveal the water cycle process and its driving mechanism.
Due to the complexity of the ET process, it is difficult to directly measure the ET: however, with the deepening of vegetation ET research, obtaining regional continuous ET data by calculating the ET becomes particularly vital.Potential evapotranspiration (ET 0 ), as the basis for calculating the ET [4], represents the maximum ET that can be achieved on a fixed underlying surface with unlimited water supply under certain meteorological conditions.ET 0 is widely used in the analysis of climate dry and wet conditions, rational use and evaluation of water resources, crop water demand and agricultural production management, and ecological environment research [5].There are many ways to calculate ET 0 , such as temperature and the solar radiation-based Hargreaves method, the equilibrium evaporation-based Priestley-Taylor model and the FAO Blaney-Criddle formula [6].However, many scholars choose the Penman-Monteith method recommended by the FAO56 to calculate ET 0 .The advantage of the Penman-Monteith method is that standard meteorological data are easily obtained or obtained through routine observation, and all calculation procedures can be standardized by the calculation of available meteorological data and time scales [7].
ET 0 was heterogeneous in temporal and spatial distribution.Jung et al. [8] pointed out that the global average annual ET increased at a rate of 7.1 ± 1.0 mm/10a from 1982 to 1997, but since 1998, the global average ET decreased by −7.9 mm/10a.Roderick et al. [9] found that in the past 50 years, the ET 0 in the northern hemisphere was decreasing at a faster rate of 2 to 4 mm/a.ET 0 also showed an increasing trend in some areas, such as West Africa [10] and Turkey [11].Wu et al. [12] pointed out that the national annual average ET 0 decreased at a rate of 0.52 mm/a based on the daily meteorological data of 552 meteorological stations in China from 1961 to 2015, and the ET 0 of most stations in arid and humid regions showed a significant decreasing trend; however, neither the increasing nor decreasing trend of ET 0 is significant in semi-arid and semi-humid regions.The decrease or increase of ET 0 may affect the hydrological cycle differently.For example, when the ET 0 decreases, the transport of water vapor in the atmosphere is reduced, resulting in corresponding changes in precipitation patterns.For agriculture, the water use efficiency of crops may be improved by reducing the adverse effects of drought on crops.However, the increase of ET 0 will increase water consumption, resulting in increased land water loss and drought, and the atmospheric circulation may be affected, which may lead to strong rainfall, thus changing the distribution of water resources.
The spatio-temporal heterogeneity of ET 0 was attributed to the common effect of different climatic factors, including maximum temperature (T max ), minimum temperature (T min ), average temperature (T mean ), RH, U and sunshine hours, while the dominant factors affecting ET 0 varied in different regions.For example, Bandyopadhyay et al. [13] concluded that the main reason for the decrease of ET 0 in India was the significant increase in RH and decrease of U by the non-parametric method of Sen's slope.Roderick et al. [14] believed that the decrease in solar radiation caused by the increase in cloud cover and aerosol concentration in the southern hemisphere was the main reason for the decrease of ET 0 in New Zealand.Hossein et al. [15] believed that U has the greatest impact on ET 0 in Iranian region under arid climate condition by sensitivity analysis.Guo et al. [16] used the Sobol global sensitivity analysis method to analyze the sensitivity of ET 0 in Australia calculated by the Penman-Monteith and Priestley-Taylor models and obtained that air temperature was the most sensitive factor to ET 0 .China has seven geographic regions, with complex and diverse climatic characteristics.ET 0 has shown a fluctuating downward trend in recent decades, such as Henan Province in Central China [17], Anhui Province in East China [18], Sichuan Basin [19] and Huaihe basin [20] in Southwest China; the results of the sensitivity-contribution rate method showed that the contribution of U to ET 0 was higher than that of sunshine hours, indicating that U was the dominant meteorological factor in the above area.On the contrary, in North China, the results of the sensitivity coefficient method showed that sunshine hours were the dominant factor, followed by U [21].Zhou et al. [22] analyzed the partial correlation between ET 0 and meteorological factors by discussing climate attribution and concluded that the increase of the T max in the Three River Headwaters region located in Northwest China was the main reason for the increasing trend of ET 0 , however, the contribution of U was the smallest among all meteorological factors, which was different from the meteorological factors affecting the trend of ET 0 in other regions.
Heilongjiang Province, located in one of the four major black soil belts in the world, has an existing land area of 454,600 km 2 and an area of 23,900 km 2 of arable land.It is the largest province in terms of grain production in the country.In recent years, Heilongjiang Province has become one of the regions with the largest temperature rise in the country and one of the regions with more serious climate disasters, threatening the safe production of food crops.The impact of floods and droughts on the agricultural economy in Heilongjiang Province accounted for about 89.4% of the total impact of various natural disasters [23]; the reduction of grain production due to floods and droughts accounted for about 12% of the total grain output in the same period [24].Nie et al. [25] speculated that the climate in Heilongjiang Province would continue to warm in the next 30 and 50 years, and water resources would be scarcer.The probability of drought is increasing.ET 0 , as a key factor to characterize the regional dry and wet conditions, is of great significance for guiding water resource management and agricultural production management in Heilongjiang Province.Jiang et al. [26] expounded on the various characteristics of ET 0 in Heilongjiang Province from the scale of the crop growing season and pointed out that ET 0 was the most sensitive to RH, which provided an important reference index for formulating a reasonable irrigation system for the crop growing season in Heilongjiang Province.However, their study was limited to the growing season of crops, which could not fully reflect the changes of ET 0 in Heilongjiang Province between inter-annual and different seasons; moreover, contributing rates of different climate factors were not analyzed.Su et al. [27] pointed out that the contribution rate of water vapor pressure and U to the ET 0 change in Heilongjiang Province was the largest by multiple regression analysis, and it has a certain significance for understanding the influencing factors on the inter-annual change of ET 0 .However, they failed to study the impact of more meteorological factors on a shorter time scale.The weather of Heilongjiang Province changes significantly in four seasons, and ET 0 may be affected by different climatic factors in seasonal scales.
This study used a linear regression equation and the Mann-Kendall method to analyze the trend of ET 0 in Heilongjiang Province in the past 60 years on the inter-annual and seasonal scales.The sensitivity coefficient and contribution of different climate factors to ET 0 were comprehensively analyzed by the sensitivity-contribution rate method, aiming to provide an important basis for the rational and efficient use of water resources and agricultural water management in different regions of Heilongjiang Province under climate change.

Overview of the Study Area and Data Sources
Heilongjiang Province, located in the east of Eurasia, is the northernmost province with the highest latitude in China, starting at 121 • 11 in the west, 135 • 05 in the east, 43 • 26 in the south, and 53 • 33 in the north, spanning 14 longitudes from east to west and 10 latitudes from north to south, with an average altitude of 481 m.It has a temperate continental monsoon climate.It spans the four major river systems of the Heilongjiang River, Wusuli River, Songhua River and Suifenhe River.The average annual temperature of Heilongjiang Province is between −5 and 5 • C, gradually increasing from the north to the south.The annual precipitation is between 400 and 650 mm, with more in the central mountainous areas, followed by the east, and less in the west and north.The annual sunshine hours in the province are mostly between 2400 and 2800 h, where the west is more than the east.
The meteorological data in this study were from the daily meteorological data of 26 meteorological stations (Figure 1) in Heilongjiang Province from 1960 to 2019 recorded by the China Meteorological Data Network, including T max , T min , T mean , RH, U, and sunshine hours; Rn was calculated according to the method recommended by the Food and Agriculture Organization of the United Nations (FAO) [7].According to the meteorological division method, March to May is divided into spring, June to August is summer, September to November is autumn, and December to February is winter [28].
nental monsoon climate.It spans the four major river systems of the Heilongjiang River, Wusuli River, Songhua River and Suifenhe River.The average annual temperature of Heilongjiang Province is between −5 and ~5 °C, gradually increasing from the north to the south.The annual precipitation is between 400 and 650 mm, with more in the central mountainous areas, followed by the east, and less in the west and north.The annual sunshine hours in the province are mostly between 2400 and ~2800 h, where the west is more than the east.
The meteorological data in this study were from the daily meteorological data of 26 meteorological stations (Figure 1) in Heilongjiang Province from 1960 to 2019 recorded by the China Meteorological Data Network, including Tmax, Tmin, Tmean, RH, U, and sunshine hours; Rn was calculated according to the method recommended by the Food and Agriculture Organization of the United Nations (FAO) [7].According to the meteorological division method, March to May is divided into spring, June to August is summer, September to November is autumn, and December to February is winter [28].

Mann-Kendall Trend and Mutation Test
For the time series X, (containing n samples), an order column is constructed: where The order column Sk is the sum of the number of values when the value of the i moment is greater than the j moment.
Under the assumption that the time series is random, define the statistic: where UF1 = 0, E(Sk) and Var(Sk) are the mean and variance of Sk, respectively, and when X1, X2,…, Xn are independent of one other, they have the same continuous distribution, which can be deduced from the following equation:

Mann-Kendall Trend and Mutation Test
For the time series X, (containing n samples), an order column is constructed: where The order column S k is the sum of the number of values when the value of the i moment is greater than the j moment.
Under the assumption that the time series is random, define the statistic: where UF 1 = 0, E(S k ) and Var(S k ) are the mean and variance of S k , respectively, and when X 1 , X 2 , . . ., X n are independent of one other, they have the same continuous distribution, which can be deduced from the following equation: UF k is a standard normal distribution, it is a sequence of statistics calculated in the order of time series X (X 1 , X 2 , . . ., X n ); given the significance level α by checking the normal distribution table, if UF i > U α , it indicates that there is a significant trend change in the sequence.Then repeat the above process in the reverse order of time series X (X n , X n−1 , . . ., X 1 ) and make Generally, if the significance level α = 0.05, then the critical value Z = ±1.96.The curves of the two statistical sequences of UF, UB and the two straight lines of ±1.96 are plotted on one plot.If the values of UF and UB are greater than 0, it indicates an upward trend in the series, and if the values of UF and UB are less than 0, it indicates a downtrend.When UF and UB exceed the critical line, it indicates a significant upward or downward trend, and the range above the critical line is determined as the time zone in which the mutation occurred.If the two curves of UF and UB intersect and the intersection point is between the critical lines, then the moment corresponding to the intersection point is the time when the mutation begins [29].

Climate Tendency Rate
The climate tendency rate was calculated by the least squares method, and the unary linear regression equation of Y i and X i was established: where A and B are regression constants, and A × 10 is the climate tendency rate, representing the changing rate of each climatic factor every 10 years (10a).A positive value indicates an increasing trend of climate in the corresponding time series, while a negative value indicates a decreasing trend.

Calculation of ET 0
The daily ET 0 of 26 sites in Heilongjiang Province was calculated using the Penman-Monteith model recommended by the FAO.The formula is as follows: where ET 0 is the potential evapotranspiration, mm; ∆ is the slope of the temperature change with the saturated water vapor pressure, kPa/ • C; Rn is the surface net radiation, MJ/(m 2 •d); G is the soil heat flux, MJ/(m 2 •d); γ is the hygrometer constant, kPa/ • C; U 2 is the U at a height of 2 m above the ground, m/s; e s is the saturated vapor pressure, kPa; e a is the actual vapor pressure, kPa; T is the temperature, • C.

Sensitivity-Contribution Rate Method Based on Partial Derivatives
Sensitivity analysis makes each input variable change within the corresponding value range, and studies and predicts the influence degree of the changes of these input variables on the output value.The influence degree is called the sensitivity coefficient [30] and is used to judge the interference degree of the relative changes of climatic factors to the changes of ET 0 .This paper mainly analyzed the sensitivity and contribution rate of T max , T min , T mean , RH, U and Rn to ET 0 .The sensitivity coefficient Sv i of ET 0 to each climatic factor was calculated by the following formula: where v i is the climatic factor; Sv i is the sensitivity coefficient of the climatic factor.The positive and negative values of Sv i reflect the correlation between ET 0 and climatic factor.A negative sensitivity coefficient indicates that ET 0 decreases with the decrease of climatic factors, and vice versa.The absolute value reflects the impact of climatic factors on ET 0 .The greater the absolute value, the greater the impact, and vice versa.The contribution rate Gv i is equal to Sv i multiplied by the annual change rate of meteorological variables (Rv i ), which is: where Rv i is the annual relative change rate of v i ; n is the total number of years; Trend vi is the annual climate tendency rate of v i , which is the slope of the univariate linear regression equation between v i and n; av vi is the annual average value of v i ; Gv i is the contribution rate, the magnitude of the absolute value of Gv i reflects the contribution of the relative change of v i to the change of ET 0 .

Data Processing
The CROPWAT 8.0 (FAO, Rome, Italy) software was used to calculate the daily ET 0 by the Penman-Monteith formula, Matlab 2021a (MathWorks, Natick, MA, USA) was used to calculate the rate of change of ET 0 , T max , T min , T mean , RH, U and Rn, and perform a Mann-Kendall mutation test.The Mann-Kendall test was used to analyze the long-term change trend and mutation of ET 0 .The ArcMap 10.4 (ESRI, Redlands, CA, USA) toolkit was used for spatial analysis to perform spatial interpolation and mapping for each meteorological variable, the spatial interpolation method was Inverse Distance Weighting (IDW) and the spatial resolution was 500 dpi.

Spatial Distribution of the Mean Values of Meteorological Factors
Spatial distribution of average T max , T min , T mean , RH, U, Rn from 1960 to 2019 is shown in Figure 2. On the inter-annual scale, the T max , T min and T mean showed an increasing trend from north to south, their inter-annual ranges were 5.5~11.5 • C, 0.5~7.0• C and 1.0~5.5 • C, respectively.A higher RH was mainly distributed in the central region, while the RH in east and west region was relatively lower.A higher U was mainly distributed in the east and west regions.Spatially, Rn increased from north to south.
On the seasonal scale, the T max , T min and T mean also showed an increasing trend from north to south.The RH was higher in summer and lower in spring.However, the U was higher in spring and lower in summer, which were 3.71 m/s and 2.67 m/s, respectively.The Rn is larger in summer and smaller in winter; the ranges were 11.5~12.8MJ/m 2 and 0.30~1.90MJ/m 2 , respectively.

Temporal and Spatial Variation of the Climate Tendency Rate of Meteorological Factors
At time series on inter-annual and seasonal scales, the T max , T min , T mean showed a significant upward trend (p < 0.05) for T max in winter (Table 1).Whereas, the RH, U and Rn showed a significant downward trend (p < 0.05), except for RH in spring and summer, U in autumn and Rn in autumn and winter.The inter-annual climate tendency rate for the T max , T min , T mean , RH, U and Rn was 0.22 • C/10a, 0.49 • C/10a, 0.36 • C/10a, −0.42 (%/10a), −0.18 (m/s/10a), and −0.04 (MJ/m 2 /10a), respectively.In terms of spatial distribution on inter-annual scales, the T max , T min , and T mean showed an increasing trend, higher climate tendency rate of them were mainly distributed in the northern region (Figure 3).The RH showed a decreasing climate tendency rate except for some central and southern regions.The U was shown a decreasing trend except Jixi, and the climate tendency rate was lower in western and higher in central region.The Rn showed a decreasing climate tendency rate except for Hulin and Tonghe.On seasonal scales, The T max , T min , and T mean showed an increasing trend except for T max in Suifenhe and T mean in Hulin in summer.The climate tendency rate of T max , T min , T mean were greatest in winter compared with the other 3 seasons (Figure 3).The RH showed a decreasing climate tendency rate except in the central regions in spring, summer, and autumn.The climate tendency rate of U was lower in the western and eastern regions, and the climate tendency rate of Rn showed a decreasing trend except for some central and eastern regions.

Temporal and Spatial Variation of ET 0
ET 0 was the highest in summer and lowest in winter, with daily averages of 4.11 mm and 0.35 mm, respectively (Figure 4).ET 0 showed decreasing trends in seasons and inter-annual; the higher ET 0 values were mainly distributed in the southwest region.The average climate tendency rate of ET 0 in spring, summer, autumn, winter and interannual was −6.33 mm/(10a), −2.72 mm/(10a), −2.58 mm/(10a), −0.65 mm/(10a), and −8.69 mm/(10a), respectively (Figure 4).In the western region, the seasonal lower ET 0 climate tendency rate, especially in spring and summer, led to a more rapid decrease trend of inter-annual ET 0 .
The inter-annual ET 0 in Heilongjiang Province showed a significant increasing trend from 1977 to 1983 (Z > 1.96), while it showed a significant decreasing trend in 2017-2019 (Z < −1.96).The mutation point of ET 0 appeared in 2011, and the changing trend of ET 0 changed from increase to decrease (Figure 5e).On a seasonal scale, the mutation years range from 2006 to 2015 (Figure 5a-d).

Sensitivity of ET 0 to Meteorological Factors
The sensitivity coefficients of ET 0 to U, Rn, RH, T max , T min and T mean are shown in Table 2. On the inter-annual scale, in the case that other climatic factors remain unchanged, when the U, Rn, T max , T mean , RH and T min increase by 10%, ET 0 will increase by 12.2%, 4.0%, 4.2%, 1.4% or decrease by 11.5% and 1.4%, respectively.The sensitivity order of ET 0 change to each climatic factor was U > RH > T max > Rn > T min = T mean .On the seasonal scale (Table 2), the changes of ET 0 in spring, summer and autumn were all negatively sensitive to RH, positively sensitive to U, Rn, T max , and T mean , while the changes of ET 0 to RH, T max , T min and T mean were negatively sensitive in winter.ET 0 was most sensitive to U in spring and autumn, and most sensitive to RH and T min in summer and winter, with sensitivity coefficients of 1.35, 1.40, −0.91 and −1.76, respectively.
Figure 6 shows the spatial distribution of sensitivity coefficients of ET 0 .On the interannual scale, The sensitivity coefficients of ET 0 to T max , T min and T mean gradually increased from north to south.The sensitivity coefficients of ET 0 to RH decreased from west to east.A higher sensitivity coefficient to U was mainly distributed in the eastern and western regions, while the sensitivity coefficient to Rn was mainly distributed in the central region.
On the seasonal scales, the sensitivity of ET 0 to T max , T min , T mean was higher in the western region in spring and summer, and higher in the southeast and south in autumn and winter.ET 0 showed lower sensitivity to RH in the central region throughout the 4 seasons, the sensitivity coefficient of ET 0 to U was higher in spring and autumn, while the sensitivity coefficient of ET 0 to Rn was higher in summer.spectively (Figure 4).In the western region, the seasonal lower ET0 climate tendency rate, especially in spring and summer, led to a more rapid decrease trend of inter-annual ET0.

Spring
The inter-annual ET0 in Heilongjiang Province showed a significant increasing trend from 1977 to 1983 (Z > 1.96), while it showed a significant decreasing trend in 2017-2019 (Z < −1.96).The mutation point of ET0 appeared in 2011, and the changing trend of ET0 changed from increase to decrease (Figure 5e).On a seasonal scale, the mutation years range from 2006 to 2015 (Figure 5a-d).

Dominant Climatic Factors for ET 0 Change
On the inter-annual scale, U was the dominant climatic factor for the change of ET 0 , followed by T max , T mean , RH, Rn and T min , with contribution rates of 6.15%, 5.03%, 4.34%, −1.82%, and −1.41%, respectively (Table 3).The positive contribution rates of RH, T max and T mean to ET 0 have not been able to offset the negative contribution rates of U, Rn and T min ; therefore, the ET 0 showed a decreasing trend from 1960 to 2019.On the seasonal scale, U was the dominant factor for the decrease of ET 0 in spring and autumn, followed by T max and RH; the largest contribution rates to the change of ET 0 were −9.79% and −13.14%, respectively.In summer, Rn was the dominant factor for ET 0 change, the contribution rates of T max , T min , and T mean to ET 0 were less different, which were 2.57%, 2.1% and 2.09%,respectively, while in winter, Rn contributes the least to ET 0 , and T min became the dominant factor.
The spatial distribution of dominant meteorological factors is shown in Figure 7.In spring, the dominant factor for ET 0 was U in the study area except for Fujin (Figure 7a), In summer, Rn was the dominant factor in most regions of the study area, T min and RH were the dominant factors in the partial northern and eastern regions.In autumn, U was the dominant factor for 85% of the total 26 sites (Figure 7c); in winter, T min was the dominant factor in the northern and eastern regions, while U was the dominant factor in the western region.On the inter-annual scale, the dominant factor of all 26 sites was U (Figure 7e).summer, Rn was the dominant factor in most regions of the study area, Tmin and RH were the dominant factors in the partial northern and eastern regions.In autumn, U was the dominant factor for 85% of the total 26 sites (Figure 7c); in winter, Tmin was the dominant factor in the northern and eastern regions, while U was the dominant factor in the western region.On the inter-annual scale, the dominant factor of all 26 sites was U (Figure 7e).

Discussion
Under the general warming of the global climate, ET 0 was on the rise in most areas, such as the annual ET 0 for India, which, as a whole, has increased in the latter half of the 20th century [31]; Awash River basin, Ethiopia [32], and South Korea over the recent 100 years [33], which is consistent with the generally accepted trend that ET 0 increases with temperature.However, the phenomenon of the "evaporation paradox" appeared in many areas around the world [34], that is, with the continuous increase of temperature, the ET 0 showed a decreasing trend, such as the Canadian prairie region, the northern region of South America, Thailand, New Zealand [9,[35][36][37], and the northwest of India [38], the Lijiang watershed [39], as well as Alor Setar, Malaysia [40].In this study, ET 0 in Heilongjiang Province from 1960 to 2019 showed a decrease trend inter-annually and seasonally with the increasing temperature, proving that the "evaporation paradox" phenomenon also existed in our study area, which was consistent with the conclusion of Li et al. [41].At present, there is still no clear conclusion about the mechanism of the "evaporation paradox"; the existing related research studies are mostly qualitative analysis, and the reasons need to be further explored.
In China, the sensitivity coefficient and contribution rate of ET 0 to meteorological factors varied in different climatic regions.Even in the same climatic region, the sensitivity coefficient and contribution rate of ET 0 were different spatially.For example, in the subtropical monsoon climate zone, ET 0 in the Yangtze River basin [42], Sichuan basin [43], and Yunnan-Guizhou plateau [44] were most sensitive to the T max , RH, and sunshine hours, respectively.In the temperate continental climate zone, ET 0 in the middle and upper reaches of the Yellow River [45] and the Ebinur Lake basin in Xinjiang [46] were most sensitive to the actual water vapor pressure and U, respectively, with the highest contribution rate of U for both regions.In the plateau mountainous climate region, the meteorological factor with the largest sensitivity coefficient in the Qinghai Tibet Plateau was the actual water vapor pressure, and the sunshine hours had the largest contribution rate to ET 0 [47].Our study area is located in the temperate monsoon climate zone, where ET 0 was the most sensitive to U and had the largest contribution rate to ET 0 .However, in the Loess Plateau basin, which is also located in the same climate zone, ET 0 was most sensitive to actual water vapor pressure and T mean had the highest contribution rate [48].In the North China Plain of the temperate monsoon climate zone, sunshine hours had the largest contribution rate to ET 0 [21].In summary, the plateau area has strong solar radiation, long sunshine time, and low temperature, which may cause the ET 0 to be most affected by the sunshine hours; most of the basin areas have a dry climate, less rainfall, and lack of water resources, causing vegetation to be more sensitive to the water condition; therefore, RH might have a larger sensitivity coefficient to ET 0 .Heilongjiang Province is located in plain and mountainous areas, and affected by southeast monsoon in summer and controlled by northwest monsoon in winter [49].This might be the reason why U contributed the most to ET 0 change.
In this study, ET 0 in Heilongjiang Province showed a decrease trend on the inter-annual scale from 1960 to 2019, with U as the dominant factor, which was consistent with the trend of ET 0 in Jilin Province analyzed by Liu et al. [50] and in Liaoning Province analyzed by Cao et al. [51].This was probably because the above three provinces are geographically contiguous and have a similar type of climate characteristics.A research established by Xu et al. in 2003 predicted that the temperature in Heilongjiang Province would increase significantly by 2030 and 2050 with a higher increase in winter and a lower increase in summer by the CCCma (Canadian Center for Climate Modelling and analysis), CCSR (Center for Climate System Research), CSIRO (Commonwealth Scientific and Industrial Research Organization), GFDL (Geophysical Fluid Dynamics Laboratory) and Hadley climate models [52].Zhang et al. [53] pointed out that the climate of Heilongjiang Province would tend to be warm and humid in the next 41 years by the Hadley model.Wang et al. [54] believed that ET 0 was not only affected by various meteorological factors, but interactions between meteorological factors would also interfere with ET 0 changes.Moreover, ET 0 may also be affected by human activities including aerosol emissions, air pollution [55] and rice area expansion [56].Therefore, further exploration in these aspects will be needed for contribution analysis of ET 0 changes in Heilongjiang Province in the future, as well as in other regions of the world.

Conclusions
In this study, ET 0 was calculated by the Penman-Monteith formula, the sensitivitycontribution rate method was used to clarify the climatic factors affecting seasonal and inter-annual changing of ET 0. , The T max , T min , T mean showed an increasing trend and RH, U, Rn showed a decreasing trend from 1960 to 2019.ET 0 showed a decreasing trend with an average climate tendency rate of −8.69 mm/(10a).The annual average ET 0 change was negatively sensitive to T min and RH, and positively sensitive to T max , T mean , U and Rn.U was the dominant factor for the ET 0 decrease in spring and autumn, while RH and T min were the dominant factors in summer and winter, respectively.On the inter-annual scale, the sensitivity order of ET 0 change to each climatic factor was U > RH > T max > Rn > T min = T mean .These results indicated that the dominant factors for ET 0 changes were U in Heilongjiang Province.Moreover, the dominant factors for ET 0 change varies under different climatic and geographical conditions; impacts of human activities on ET 0 should also be considered in future studies.

Figure 1 .
Figure 1.Map of the study area and distribution of meteorological stations in Heilongjiang Province.

Figure 1 .
Figure 1.Map of the study area and distribution of meteorological stations in Heilongjiang Province.

Figure 4 .
Figure 4. (Row A) spatial distribution and (Row B) climate tendency of ET0 from 1960 to 2019.Figure 4. (Row A) spatial distribution and (Row B) climate tendency of ET 0 from 1960 to 2019.

Figure 4 .
Figure 4. (Row A) spatial distribution and (Row B) climate tendency of ET0 from 1960 to 2019.Figure 4. (Row A) spatial distribution and (Row B) climate tendency of ET 0 from 1960 to 2019.

Figure 7 .
Figure 7. Spatial distribution of dominant climatic factors in (a) spring, (b) summer, (c) autumn, (d) winter and (e) inter-annually in the study area from 1960 to 2019.

Figure 7 .
Figure 7. Spatial distribution of dominant climatic factors in (a) spring, (b) summer, (c) autumn, (d) winter and (e) inter-annually in the study area from 1960 to 2019.

Author Contributions:
Conceptualization, T.N.; methodology, T.N. and R.Y.; software, R.Y.; validation, R.Y., S.L. X.Z. and Y.L.; formal analysis, T.N. and R.Y.; data curation, P.C. and T.L.; writingoriginal draft preparation, T.N. and R.Y.; writing-review and editing, Z.G., C.D. (Chong Du), C.D. (Changlei Dai) and H.J.; visualization, R.Y.; supervision, T.N., Z.Z. and Z.G.; funding acquisition, T.N. and Z.Z.All authors have read and agreed to the published version of the manuscript.Funding: This work was fund by the Opening Project of Key Laboratory of Efficient Use of Agricultural Water Resources, Ministry of Agriculture and Rural Affairs of the People's Republic of China (number: AWR2021002), the Basic Scientific Research Fund of Heilongjiang Provincial Universities (number: 2021-KYYWF-0019), the National Key Research and Development Program of China (2021YFD1500802) and the National Natural Science Foundation Project of China (numbers: 51779046 & 52079028).

Table 1 .
Seasonal and inter-annual climate tendency rate of meteorological factors in Heilongjiang Province.

Table 2 .
Sensitivity coefficients of seasonal and inter-annual ET 0 change to climate factors.

Table 3 .
Contribution of seasonal and inter-annual climate factors to ET 0 changes from 1960 to 2019.