Spatiotemporal Changes of Reference Evapotranspiration in the Highest-Latitude Region of China

Reference evapotranspiration (ET0) is often used to make management decisions for crop irrigation scheduling and production. In this study, the spatial and temporal trends of ET0 in China’s most northern province as well as the country’s largest agricultural region were analyzed for the period from 1964 to 2013. ET0 was calculated with the Penman-Monteith of Food and Agriculture Organization of the United Nations irrigation and drainage paper NO.56 (FAO-56) using climatic data collected from 27 stations. Inverse distance weighting (IDW) was used for the spatial interpolation of the estimated ET0. A Modified Mann–Kendall test (MMK) was applied to test the spatiotemporal trends of ET0, while Pearson’s correlation coefficient and cross-wavelet analysis were employed to assess the factors affecting the spatiotemporal variability at different elevations. The results from this study showed a clear decreasing trend for annual ET0 from the low elevation plain area to the high elevation mountainous area. Over the past five decades, ET0 in Heilongjiang Province decreased in all seasons, except for the winter months, during which a steady increase in temperature was found. Elevation played an important role in estimating ET0 in this higher-latitude region, while relative humidity was the most relevant meteorological factor that affected the spatiotemporal variation of ET0 in the province. Overall, the findings from the study suggest that winter ET0 in a high altitude region will continue to increase in the future as climate change persists, which could worsen spring droughts and irrigation management for semi-arid areas in the province.


Introduction
Reference evapotranspiration (ET 0 ) is an important parameter in the water and energy balance of terrestrial ecosystems [1,2].It has been used for determining crop water demand and irrigation scheduling [3,4], as it is closely related to hydrological processes such as soil moisture transfer to the atmosphere, groundwater recharge, and runoff generation [5,6].ET 0 can be also considered an indicator for the ecological water demand of natural landscapes.Therefore assessing ET 0 spatial distribution can be useful for determining regional water consumption and developing effective water management strategies.
Although global temperatures have increased in the past 50 years, ET 0 has been found to respond differently in different geographical regions.For instance, ET 0 has been reported to increase in the Mediterranean countries and Central Asia [7,8].However, ET 0 has been found to decline in Australia [9], the United States [10], Canada [11], India [12], and China [13,14].Despite a number of large-scale studies on ET 0 , few studies have looked at its spatiotemporal trend in high-altitude regions that are vulnerable to temperature changes in the future.
In the context of a hydrological cycle, any change in meteorological factors could affect ET 0 because evapotranspiration is directly or indirectly related to these factors [15].Previous studies have analyzed the sensitivity of ET 0 to various meteorological factors such as air temperature, humidity, wind speed, or solar radiation [16][17][18][19].Roderick and Farquhar [20] found that the decrease in ET 0 was consistent with the decreases in sunlight resulting from increasing cloud coverage and aerosol concentration.Qian et al. [21] reported that increased air pollution may decrease solar radiation and result in less energy being available for evapotranspiration.Differing from most opinions, McVicar et al. [22] considered that change in wind speed was one of the most important factors contributing to the decline of the ET 0 on the global scale.In addition, it may also be relevant that solar radiation was a vital factor likely to effect ET 0 .However, Brutsaert et al. [23] argued that, as actual evaporation increases, atmospheric humidity increases, thus reducing the apparent potential evaporation, which is reflected in reduced pan evaporation.In short, the cause of the change in ET 0 is complex and remains a scientific problem, especially in high-altitude regions, due to these regions experiencing the fastest warming trends [9,24].
As one of the most important grain production bases in China, Heilongjiang is the country's northernmost province (43)(44)(45)(46)(47)(48)(49)(50)(51)(52)(53) • N) and is also the world's northernmost rice-growing area.Several studies have found a significant increase in temperature in the region during the past several decades [13,14,25].In addition, a few previous studies had qualitatively indicated that annual ET 0 had a decreasing trend over the whole Northeast of China [18,19,[26][27][28].However, few studies quantitatively analyze ET 0 variation against elevation in high altitude regions like northeast China.Furthermore, it is not clear whether and how the temperature increase affects ET 0 variability in this province.The variation of ET 0 is a key element, which is closely related to water demands for irrigation.As global temperatures are projected to increase in the coming century [24], it is vital to understand the spatiotemporal variation of ET 0 in high latitude regions.This is especially true for analyzing future water demands for irrigation in these regions as many of the regions suffer from periodic droughts [3,29,30].
This study aims to assess the past five decades of spatiotemporal change in reference to evapotranspiration in China's highest latitude region.Specifically, the study aims to (1) quantify ET 0 at different elevations in Heilongjiang Province for the period from 1964 to 2013; (2) map ET 0 spatially across the province to discern the regional heterogeneity; and (3) analyze the sensitivity of ET 0 to key meteorological variables.

Study Area
Heilongjiang Province (43 • 22 -53 • 24 N and 121 • 13 -135 • 05 E) is located in northeast China, covering a total land area of approximately 473,000 km 2 (Figure 1a).About one third of the province (i.e., 37%) consists of low-lying plain areas with extensive wetlands, many of which have been converted to rice cultivation in the past half century [31].Currently, 28% of the province's area (i.e., 133,000 km 2 ) is agricultural land (Figure 1b) that annually produces a total of 0.65 × 10 8 t grains, ranking first in China (Heilongjiang Bureau of Statistics, 2015) [32].The province is designed by the Chinese Government to be the most important grain production base in the future.and 50 × 10 8 J/m 2 , 55% to 60% of which is accumulated during the growing season from May to September, with a decreasing trend from South to North.The total duration of sunshine in the province is about 2400-2800 h per year.The monthly average wind speed is 2.7 m/s, varying from 2.1 m/s to 3.8 m/s.Normally, strong winds occur in Spring, causing lower humidity in the season.

Meteorological Data Collection
Daily meteorological data from 27 stations in Heilongjiang Province were obtained for the period from 1964 to 2013 from the China Meteorological Administration (http://data.cma.cn/)(Table 1).Six meteorological variables were used to calculate ET0 for each station, and they included maximum air temperature (Tmax, °C), minimum air temperature (Tmin, °C), mean air temperature (Tmean, °C), mean relative humidity (RH, %), wind speed at 10 m above the ground (WS, m/s), and sunshine duration (SD, h/d).All these meteorological data were qualitatively checked before use.

Calculation of Reference Evapotranspiration
In this study, the Penman-Monteith equation, which is applied widely all over world [34] as recommended by the United Nations Food and Agriculture Organization (FAO), was used for calculating the ET0 due to lack of measured pan evapotranspiration: The climate of this high-latitude region can be characterized as temperate, humid, and semi-humid continental monsoon with a long-term annual average temperature of 1.9 • C, which ranges from −5 to 5 • C, and a long-term annual average precipitation of about 550 mm, most of which falls during the summer season [33].The total annual solar radiation in the province is between 44 × 10 8 and 50 × 10 8 J/m 2 , 55% to 60% of which is accumulated during the growing season from May to September, with a decreasing trend from South to North.The total duration of sunshine in the province is about 2400-2800 h per year.The monthly average wind speed is 2.7 m/s, varying from 2.1 m/s to 3.8 m/s.Normally, strong winds occur in Spring, causing lower humidity in the season.

Meteorological Data Collection
Daily meteorological data from 27 stations in Heilongjiang Province were obtained for the period from 1964 to 2013 from the China Meteorological Administration (http://data.cma.cn/)(Table 1).Six meteorological variables were used to calculate ET 0 for each station, and they included maximum air temperature (T max , • C), minimum air temperature (T min , • C), mean air temperature (T mean , • C), mean relative humidity (RH, %), wind speed at 10 m above the ground (WS, m/s), and sunshine duration (SD, h/d).All these meteorological data were qualitatively checked before use.

Calculation of Reference Evapotranspiration
In this study, the Penman-Monteith equation, which is applied widely all over world [34] as recommended by the United Nations Food and Agriculture Organization (FAO), was used for calculating the ET 0 due to lack of measured pan evapotranspiration: where ET 0 is the reference evapotranspiration (mm d −1 ), R n is net radiation (MJ m −2 d −1 ); G is the soil heat flux density (MJ m −2 d −1 ), T mean is the mean air temperature ( • C), e s is the saturation vapor pressure (kPa), e a is the actual vapor pressure (kPa), ∆ is the slope of the saturation vapor pressure function (kPa • C −1 ), γ is the psychometric constant (kPa • C −1 ), and U 2 is the wind speed at 2 m height (m s −1 ).Due to the lack of U 2 data, this variable was estimated from wind speed at 10 m height, using the logarithmic vertical wind speed profile [34]: where z is the height of measurement above ground surface (m) and U z is the measured wind speed at z m above ground surface (m s −1 ).The detailed description used in estimating ET 0 was given in Chapter 3 of FAO-56 [34].The seasonal and annual values of ET 0 were calculated using the daily ET 0 data.

Spatial Interpolation
In recent years, based on meteorological station data, some spatial interpolation methods have been applied to study the spatial distributions of ET 0 and the climate variations.Due to the fact that the inverse distance weighting (IDW) method gave the lowest mean error, the method is applied widely all over world [3,27,35,36].Thus, in this study, the IDW method is selected to analyze the spatial variation of ET 0 .
The weights of the IDW interpolation method are solely a function of the distance between the interest point and the sampling points for i = 1, 2, . . ., n. Considering the distance D i between these two points, the value of an interest point is in the form [37]: where M is the interpolated value of a interest point, M i is the value of the sampling point I, and D i is the distance to the sampling point.

Data Analysis
In this study, all analyses of the spatiotemporal variation of ET 0 were performed two times; once with annual data and once with seasonal data.Meanwhile, a season is defined in the standard climatological way; Spring is considered as occurring from March to May, Summer from June to August, Autumn from September to November, and Winter from December to February.
The Modified Mann-Kendall test (the Mann-Kendall test with trend-free pre-whitening) was applied to analyze the long-term trend of ET 0 throughout the whole study period.The traditional Mann-Kendall test assumes that the series is independent and not robust against autocorrelation.However, some meteorological time series may frequently show statistically significant serial correlation.This may lead to a disproportionate rejection of the null hypothesis without the trend, Water 2017, 9, 493 5 of 17 whereas the null hypothesis is actually true.Therefore, the effect of serial correlation is a main source of uncertainty in testing and interpreting trends.To correct the serial correlation in the data, a trend-free pre-whitening procedure is used to remove the lag-1 serial correlation (r 1 ) from the time series.Recently, this method has been applied in more and more studies [25,38,39].To determine whether the observed data set is serially correlated, the significance of the lag-1 serial correlation (r 1 ) should be tested at the 0.10 significance level.r 1 is calculated using the following formula [40]: , the time series is assumed to be independent at the 0.10 significance level and can be subjected to the original Mann-Kendall test.Otherwise, the effect of serial correlation should be removed from the time series by pre-whitening prior to applying the Mann-Kendall test.The Mann-Kendall test was then used to detect a trend in the residual series.The new time series was obtained as follows [41]: The r 1 value of this new time data set was calculated and used to determine the residual series as: The value of β × i was added again to the residual data set as follows: The y i series was subjected to trend analysis.
Based on the Modified Mann-Kendall test, a normalized test statistic Z was calculated for all the meteorological stations.When the significance levels are set at 0.01 and 0.05, |Z a | is 2.58 and 1.96, respectively.At a certain significance level, if |Z| > |Z a |, the null hypothesis H 0 is rejected.In other words, the trend is significant at the set level of significance.Otherwise, no significant trend exists [25].
Pearson's correlation coefficient was applied to analyze the correlativity between the ET 0 and key meteorological factors, and r value was calculated via SPSS 16.0 software.
The sensitivity coefficients of the ET 0 to each meteorological factor were simulated as follows [42,43]: where C s is the sensitivity coefficient, CH ET0 is the change in ET 0 with respect to change in the climate variable, and CH CV is the change (increase or decrease) in the climate variable, which was changed within a range between −20% and 20% at an interval of 1%.
In order to analyze the temporal evolution of the different components at a number of scales (i.e., scale analysis) that are present in the time series, wavelet coherence (WTC), which is widely used in hydrological and meteorological research all over the world, was applied to determine the intensity of the covariance of two time series, x t and y t , in time-frequency domains [44].Ranging from 0 to 1, with 1 representing the highest covariance between x t and y t , the values of the wavelet squared coherency R were calculated as [45]: Water 2017, 9, 493 6 of 17 where W x t and W y t are the wavelet transforms of time series x and y, respectively, according to time t, and W xy t is their cross-wavelet spectrum.S represents a smoothing operator and can be expressed as: where S scale and S time represent the smoothing along the wavelet scale and in time, respectively.For the Morlet wavelet, a suitable smoothing operator is given by Torrence and Webster as follows: where c 1 and c 2 are normalization constants and Π is the rectangle function.The factor of 0.6 is the empirically determined scale decorrelation length for the Merlet wavelet [44].
In addition, the Lilliefors test and the Kolmogorov-Smirnov goodness-of-fit test were applied to check the probability distribution of each time series.Furthermore, it was worth noting that a non-normal time series was transformed into a record of percentiles in terms of its cumulative distribution function in order to avoid the presence of any oulier [45].The WTC procedures were applied after the transformation.In this study, the wavelet analysis was performed by a Matlab toolbox developed by Grinsted et al. [45], which is available at http://noc.ac.uk/using-science/crosswaveletwavelet-coherence.

Spatial Variation of ET 0
The long-term (1964-2013) ET 0 in Heilongjiang Province averaged 729 mm/year, ranging from 542 mm/year to 960 mm/year.The annual average ET 0 showed a clear decreasing trend from the low elevation plain area to the high elevation mountainous area (Figure 2).The annual average ET 0 was 772.5 mm/year in the elevation from 0 to 200 m, 705.4 mm/year in the elevation from 200 to 400 m, and 668.1 mm/year in the elevation above 400 m (Table 2).The ET 0 in the Southwest was higher than that in other areas.The annual average ET 0 was highest (960 mm) at the Tailai station in the Southwest and lowest (542 mm) at the Mohe station in the Northwest.
Water 2017, 9, 493 6 of 16 where Sscale and Stime represent the smoothing along the wavelet scale and in time, respectively.For the Morlet wavelet, a suitable smoothing operator is given by Torrence and Webster as follows: where c1 and c2 are normalization constants and Π is the rectangle function.The factor of 0.6 is the empirically determined scale decorrelation length for the Merlet wavelet [44].
In addition, the Lilliefors test and the Kolmogorov-Smirnov goodness-of-fit test were applied to check the probability distribution of each time series.Furthermore, it was worth noting that a nonnormal time series was transformed into a record of percentiles in terms of its cumulative distribution function in order to avoid the presence of any oulier [45].The WTC procedures were applied after the transformation.In this study, the wavelet analysis was performed by a Matlab toolbox developed by Grinsted et al. [45], which is available at http://noc.ac.uk/using-science/crosswavelet-waveletcoherence.

Spatial Variation of ET0
The long-term (1964-2013) ET0 in Heilongjiang Province averaged 729 mm/year, ranging from 542 mm/year to 960 mm/year.The annual average ET0 showed a clear decreasing trend from the low elevation plain area to the high elevation mountainous area (Figure 2).The annual average ET0 was 772.5 mm/year in the elevation from 0 to 200 m, 705.4 mm/year in the elevation from 200 to 400 m, and 668.1 mm/year in the elevation above 400 m (Table 2).The ET0 in the Southwest was higher than that in other areas.The annual average ET0 was highest (960 mm) at the Tailai station in the Southwest and lowest (542 mm) at the Mohe station in the Northwest.Over the past five decades, the ET0 in Heilongjiang Province averaged 248 mm, 312 mm, 134 mm, and 35 mm in Spring, Summer, Autumn, and Winter, respectively.(Figure 3).The ET0 during summer was clearly higher in the low elevation plain area than in the high elevation mountainous area.The spring, summer, and autumn ET0 showed a decreasing trend with increasing elevation.Over the past five decades, the ET 0 in Heilongjiang Province averaged 248 mm, 312 mm, 134 mm, and 35 mm in Spring, Summer, Autumn, and Winter, respectively.(Figure 3).The ET 0 during summer was clearly higher in the low elevation plain area than in the high elevation mountainous area.The spring, summer, and autumn ET 0 showed a decreasing trend with increasing elevation.However, the winter ET 0 showed a tendency to increase first and decrease thereafter as elevation increased.
Water 2017, 9, 493 7 of 16 However, the winter ET0 showed a tendency to increase first and decrease thereafter as elevation increased.Note: a.s.l.means above sea level.

Temporal Variation of ET0
During 1964 to 2013, the average annual ET0 in Heilongjiang Province varied from 722 mm to 633 mm, with a decreasing trend of 7.4 mm ever 10 years (Table 2).Furthermore, different decreasing trends were shown for different elevations (Figure 4a).The trend decreased trend with an increase of elevation (Table 2).The decreasing trend of ET0 was very significant (p < 0.01 or p < 0.05) for most regions from 0 to 400 m and not significant for the northern region above 400 m aor for a few local regions (Figure 4b).

Temporal Variation of ET 0
During 1964 to 2013, the average annual ET 0 in Heilongjiang Province varied from 722 mm to 633 mm, with a decreasing trend of 7.4 mm ever 10 years (Table 2).Furthermore, different decreasing trends were shown for different elevations (Figure 4a).The trend decreased trend with an increase of elevation (Table 2).The decreasing trend of ET 0 was very significant (p < 0.01 or p < 0.05) for most regions from 0 to 400 m and not significant for the northern region above 400 m aor for a few local regions (Figure 4b).However, the winter ET0 showed a tendency to increase first and decrease thereafter as elevation increased.Note: a.s.l.means above sea level.

Temporal Variation of ET0
During 1964 to 2013, the average annual ET0 in Heilongjiang Province varied from 722 mm to 633 mm, with a decreasing trend of 7.4 mm ever 10 years (Table 2).Furthermore, different decreasing trends were shown for different elevations (Figure 4a).The trend decreased trend with an increase of elevation (Table 2).The decreasing trend of ET0 was very significant (p < 0.01 or p < 0.05) for most regions from 0 to 400 m and not significant for the northern region above 400 m aor for a few local regions (Figure 4b).Seasonally, the average ET0 varied from 264 mm to 221 mm in Spring, from 332 mm to 280 mm in Summer, from 144 mm to119 mm in Autumn, and from 40 mm to 31 mm in Winter, with a changing trend of -4.2 mm every 10 years, -2.7 mm every 10 years, and -1.1 mm every 10 years, and 0.6 mm every 10 years during the period from 1964 to 2013 in Heilongjiang Province, respectively (Table 3).Seasonally, the average ET 0 varied from 264 mm to 221 mm in Spring, from 332 mm to 280 mm in Summer, from 144 mm to119 mm in Autumn, and from 40 mm to 31 mm in Winter, with a changing trend of -4.2 mm every 10 years, -2.7 mm every 10 years, and -1.1 mm every 10 years, and 0.6 mm every 10 years during the period from 1964 to 2013 in Heilongjiang Province, respectively (Table 3).
Furthermore, the same trend as for the annual ET 0 occured for the ET 0 during Spring and Summer, as shown for different elevations in Figures 5a and 6a.The trend is declining with the increase of elevation (Table 3).The decreasing trend of ET 0 is very significant (p < 0.01 or p < 0.05) for most regions from 0 to 400 m and not significant for a few local regions (Figures 5b and 6b).A decreasing trend of ET 0 in Autumn was shown from 0 to 400 m from 1964 to 2013.However, an increasing trend is shown above 400 m (Figure 7a).The decreasing trend of ET 0 is significant (p < 0.05) for most regions from 0 to 400 m and very significant (p < 0.01) in a few local regions.However, no significant level is shown (Figure 7b).A tendency for ET 0 to increase in Winter is shown for different elevations from 1964 to 2013 in Figure 8a.In addition, it is worth noting that the trend is most obvious for the region above 400 m, with an increasing trend 0.9 mm every 10 years (Figure 8b).Furthermore, the same trend as for the annual ET0 occured for the ET0 during Spring and Summer, as shown for different elevations in Figures 5a and 6a.The trend is declining with the increase of elevation (Table 3).The decreasing trend of ET0 is very significant (p < 0.01 or p < 0.05) for most regions from 0 to 400 m and not significant for a few local regions (Figures 5b and 6b).A decreasing trend of ET0 in Autumn was shown from 0 to 400 m from 1964 to 2013.However, an increasing trend is shown above 400 m (Figure 7a).The decreasing trend of ET0 is significant (p < 0.05) for most regions from 0 to 400 m and very significant (p < 0.01) in a few local regions.However, no significant level is shown (Figure 7b).A tendency for ET0 to increase in Winter is shown for different elevations from 1964 to 2013 in Figure 8a.In addition, it is worth noting that the trend is most obvious for the region above 400 m, with an increasing trend 0.9 mm every 10 years (Figure 8b).Furthermore, the same trend as for the annual ET0 occured for the ET0 during Spring and Summer, as shown for different elevations in Figures 5a and 6a.The trend is declining with the increase of elevation (Table 3).The decreasing trend of ET0 is very significant (p < 0.01 or p < 0.05) for most regions from 0 to 400 m and not significant for a few local regions (Figures 5b and 6b).A decreasing trend of ET0 in Autumn was shown from 0 to 400 m from 1964 to 2013.However, an increasing trend is shown above 400 m (Figure 7a).The decreasing trend of ET0 is significant (p < 0.05) for most regions from 0 to 400 m and very significant (p < 0.01) in a few local regions.However, no significant level is shown (Figure 7b).A tendency for ET0 to increase in Winter is shown for different elevations from 1964 to 2013 in Figure 8a.In addition, it is worth noting that the trend is most obvious for the region above 400 m, with an increasing trend 0.9 mm every 10 years (Figure 8b).

Sensitivity Analysis of ET0 with Climatic Factors
Climatic factors, including wind speed (WS), maximum and minimum temperature (Tmax and Tmin), sunshine duration (SD), and relative humidity (RH), in Heilongjiang Province showed a clear association with elevation change.The average WS, Tmax, Tmin, and SD all had a decreasing trend from the low elevation plain area to the high elevation mountainous area in the province across the four seasons (Table 4).By contrast, the average relative humidity (RH) had a reverse relation with  Note: ** indicates a 1% significance level; Z is the statistical value by the MMK method.

Sensitivity Analysis of ET0 with Climatic Factors
Climatic factors, including wind speed (WS), maximum and minimum temperature (Tmax and Tmin), sunshine duration (SD), and relative humidity (RH), in Heilongjiang Province showed a clear association with elevation change.The average WS, Tmax, Tmin, and SD all had a decreasing trend from the low elevation plain area to the high elevation mountainous area in the province across the four seasons (Table 4).By contrast, the average relative humidity (RH) had a reverse relation with

Sensitivity Analysis of ET 0 with Climatic Factors
Climatic factors, including wind speed (WS), maximum and minimum temperature (T max and T min ), sunshine duration (SD), and relative humidity (RH), in Heilongjiang Province showed a clear association with elevation change.The average WS, T max , T min , and SD all had a decreasing trend from the low elevation plain area to the high elevation mountainous area in the province across the four Water 2017, 9, 493 10 of 17 seasons (Table 4).By contrast, the average relative humidity (RH) had a reverse relation with elevation.With respect to the seasonal patterns at different elevations, the variation of the meteorological factors in Winter was complex.The average T max and SD in Winter showed a decreasing trend as the elevation increased.The average WS and T min in Winter showed a tendency to decrease first and then to increase with increasing elevation.However, the average RH in Winter showed a tendency to increase first and then to decrease (Table 4).The variation of interannual and seasonal WS, T max , T min , SD, and RH showed a complex trend from 1964 to 2013 at different elevations in Heilongjiang Province.The change of T min was most significant, and the change of SD was not significant at different elevations (Table 5).The correlation coefficients between ET 0 and the other meteorological factors (WS, T max , T min , SD, and RH) was calculated in the interannual category, Spring, Summer, Autumn, and Winter from 1964 to 2013 at different elevations in Heilongjiang Province, northeast China (Table 6).The most relevant meteorological factor to ET 0 for Heilongjiang Province was RH annually and in Spring, Summer, and Autumn and T max in Winter.In addition, the most relevant meteorological factor changed for different elevations, i.e., on the annual scale the most sensitive meteorological factor was RH in the elevation range of 0 m to 200 m and above 400 m but SD in the elevation range of 200 m to 400 m.On an annual scale, the most sensitive variable was RH, followed by T max at different elevations in the whole of Heilongjiang Province.On the seasonal scale, the most sensitive variable was RH in Spring, Summer, Autumn and T max in Winter.In addition, WS also had a relatively large sensitivity coefficient on ET 0 , except in Winter.SD has a relatively large sensitivity coefficient on ET 0 in Summer.There was no obvious pattern of change at different elevations, except that ET 0 's sensitivity coefficients to T max and T min in Winter increased with increasing elevation (Table 7).

Discussion
In this study, we found a clear spatial pattern of ET 0 in China's far northern province.Across the high-latitude region, there was a decreasing trend of annual ET 0 from the low elevation plain area to the high elevation mountainous area, which is more apparent in Spring, Summer, and Autumn.This may be attributed to a combined effect of temperature and precipitation variations since they are closely related to the topography and monsoon circulation in the region [18,35,46].Similar to two recent studies [13,14], our study found that there was a significantly decreasing trend of annual ET 0 in the typically high altitude region.The same trend was found in Spring, Summer, and Autumn.In addition, some studies that reported ET 0 was significantly reduced in high altitude regions such as the Taoer River Basin [47], the northwest of China [48], and so on.However, we found an increasing trend of ET 0 in Winter.Furthermore, the higher the elevation, the more obvious the increasing trend (Figure 8a).In addition, as discussed in the introduction, ET 0 has been found to respond differently in different geographical regions.One of reasons is that ET 0 had different sensitivities to meteorological variables in different regions [43,49,50].
Our study showed that ET 0 in Heilongjiang was sensitive to the key meteorological variables and the sensitivity varies depending on the season.ET 0 in Spring, Summer, and Autumn was most sensitive to relative humidity in high altitude region.However, ET 0 was most sensitive to T max in the wintertime, followed by T min .These findings differ from previous studies [13,51] on evapotranspiration across China, which reported that ET 0 was only sensitive to relative humidity in northeast China.There was a significantly increased trend (p < 0.05) in winter T max (Figure 9a) and a very significantly increased trend in winter T min , which was almost double the T max (p < 0.01) (Figure 9b) at different elevations in Heilongjiang Province.In order to determine how winter temperatures affect ET 0 specifically, we applied wavelet coherence (WTC).ET 0 has a correlation of a 0.05 significance level with T max (Figure 10a) and T min (Figure 10b) in Winter on a time scale of up to nine years or three to nine years for the areas below 400 m a.s.l.from 1964 to 2013.The time scale of significance correlation was complex above 400 m a.s.l.As shown in the Figure 10a,b, there were correlations of a 0.05 significance level in a time scale of up to five years from 1970 to 1985, six to eight years from 1976 to 1990, and up to seven years from 1990 to 2005.In addition, all these areas showed an in-phase (the arrow pointed to the right), which indicated that T max and T min had a strong effect on winter ET 0 at this time scale.This explained why ET 0 in winter has had an increasing trend in the past five decades, and we expect that the trend will continue in the future as climate change persists, which could worsen spring droughts and irrigation management for semi-arid areas in the province.Some uncertainty may exist in estimating reference evapotranspiration rates with the Penman-Monteith equation.However, due to the lack of actual pan evapotranspiration measurements in this region, we were unable to quantify possible error estimates.Despite the uncertainties and limitations, our study sheds light on the spatiotemporal variation of ET 0 and its sensitivity to meteorological variables at different elevations in the highest latitude region of China.Future studies are needed to constrain the uncertainties and refine ET 0 estimation for different land use types and crop cultivations.

Summary and Conclusions
This study analyzed the spatiotemporal variation of seasonal and interannual reference evapotranspiration in China's far northern region, Heilongjiang Province, at different elevations from 1964 to 2013.Based on the findings, the following conclusions can be drawn: (1) Reference evapotranspiration in Heilongjiang Province has decreased in the past five decades during all seasons, except Winter.This is mainly a result of the steadily-increasing winter temperature in this higher-latitude region.(2) Elevation plays an important role in estimating reference evapotranspiration in Heilongjiang because of the greater sensitivity of transpiration to temperature in the higher-altitude region.(3) Over the past five decades, the greatest changes in Heilongjiang's climate have occurred in relative humidity and temperatures (especially minimum temperature) during the winter months.
Relative humidity is the most important meteorological factor that has affected the spatiotemporal variation of reference evapotranspiration in this province.

Figure 1 .
Figure 1.(a) Geographical locations and (b) land use types of the meteorological stations used in this study for Heilongjiang Province, northeast China.

Figure 1 .
Figure 1.(a) Geographical locations and (b) land use types of the meteorological stations used in this study for Heilongjiang Province, northeast China.

Figure 2 .
Figure 2. Spatial variation of annual average ET0 from 1964 to 2013 in Heilongjiang Province, northeast China.

Figure 2 .
Figure 2. Spatial variation of annual average ET 0 from 1964 to 2013 in Heilongjiang Province, northeast China.

Figure 4 .
Figure 4. (a) Temporal trends of interannual ET0 for different elevations during 1964 to 2013 and (b) the spatial pattern of the Modified Mann-Kendall (MMK) test in Heilongjiang Province, northeast China.Down ** indicates the decreasing trend at the 0.01 significance level, Down * indicates the decreasing trend at 0.05 significance level, and not significant means no significantly decreasing trend.(Green thick and solid line is the linear trend of ET0 at 0-200 m, the Red thick and solid line is the linear trend of ET0 at 200-400 m, Blue thick and solid line is the linear trend of ET0 at >400 m).Seasonally, the average ET0 varied from 264 mm to 221 mm in Spring, from 332 mm to 280 mm in Summer, from 144 mm to119 mm in Autumn, and from 40 mm to 31 mm in Winter, with a changing trend of -4.2 mm every 10 years, -2.7 mm every 10 years, and -1.1 mm every 10 years, and 0.6 mm every 10 years during the period from 1964 to 2013 in Heilongjiang Province, respectively (Table3).

Figure 4 .
Figure 4. (a) Temporal trends of interannual ET0 for different elevations during 1964 to 2013 and (b) the spatial pattern of the Modified Mann-Kendall (MMK) test in Heilongjiang Province, northeast China.Down ** indicates the decreasing trend at the 0.01 significance level, Down * indicates the decreasing trend at 0.05 significance level, and not significant means no significantly decreasing trend.(Green thick and solid line is the linear trend of ET0 at 0-200 m, the Red thick and solid line is the linear trend of ET0 at 200-400 m, Blue thick and solid line is the linear trend of ET0 at >400 m).

Figure 4 .
Figure 4. (a) Temporal trends of interannual ET 0 for different elevations during 1964 to 2013 and (b) the spatial pattern of the Modified Mann-Kendall (MMK) test in Heilongjiang Province, northeast China.Down ** indicates the decreasing trend at the 0.01 significance level, Down * indicates the decreasing trend at 0.05 significance level, and not significant means no significantly decreasing trend.(Green thick and solid line is the linear trend of ET 0 at 0-200 m, the Red thick and solid line is the linear trend of ET 0 at 200-400 m, Blue thick and solid line is the linear trend of ET 0 at >400 m).

Figure 5 .Figure 6 .
Figure 5. (a) Temporal trends of spring ET0 for different elevations from 1964 to 2013 and (b) spatial pattern of MMK test in Heilongjiang Province, northeast China.Down ** indicates the decreasing trend at the 0.01 significance level, Down * indicates the decreasing trend at 0.05 significance level, and not significant means no significantly decreasing trend.(Green thick and solid line is the linear trend of ET0 at 0-200 m, the Red thick and solid line is the linear trend of ET0 at 200-400 m, Blue thick and solid line is the linear trend of ET0 at >400 m).

Figure 5 .
Figure 5. (a) Temporal trends of spring ET 0 for different elevations from 1964 to 2013 and (b) spatial pattern of MMK test in Heilongjiang Province, northeast China.Down ** indicates the decreasing trend at the 0.01 significance level, Down * indicates the decreasing trend at 0.05 significance level, and not significant means no significantly decreasing trend.(Green thick and solid line is the linear trend of ET 0 at 0-200 m, the Red thick and solid line is the linear trend of ET 0 at 200-400 m, Blue thick and solid line is the linear trend of ET 0 at >400 m).

Figure 5 .Figure 6 .
Figure 5. (a) Temporal trends of spring ET0 for different elevations from 1964 to 2013 and (b) spatial pattern of MMK test in Heilongjiang Province, northeast China.Down ** indicates the decreasing trend at the 0.01 significance level, Down * indicates the decreasing trend at 0.05 significance level, and not significant means no significantly decreasing trend.(Green thick and solid line is the linear trend of ET0 at 0-200 m, the Red thick and solid line is the linear trend of ET0 at 200-400 m, Blue thick and solid line is the linear trend of ET0 at >400 m).

Figure 6 .
Figure 6.(a) Temporal trends of summer ET 0 for different elevations from 1964 to 2013 and (b) spatial pattern of MMK test in Heilongjiang Province, northeast China.Down ** indicates the decreasing trend at the 0.01 significance level, Down * indicates the decreasing trend at 0.05 significance level, and not significant means no significantly decreasing trend.(Green thick and solid line is the linear trend of ET 0 at 0-200 m, the Red thick and solid line is the linear trend of ET 0 at 200-400 m, Blue thick and solid line is the linear trend of ET 0 at >400 m).

Figure 7 .Figure 8 .
Figure 7. (a) Temporal trends of autumn ET0 for different elevations from 1964 to 2013 and (b) spatial pattern of MMK test in Heilongjiang Province, northeast China.Down ** indicates the decreasing trend at the 0.01 significance level, Down * indicates the decreasing trend at 0.05 significance level, and not significant means no significantly decreasing trend.(Green thick and solid line is the linear trend of ET0 at 0-200 m, the Red thick and solid line is the linear trend of ET0 at 200-400 m, Blue thick and solid line is the linear trend of ET0 at >400 m).

Figure 7 .Figure 7 .Figure 8 .
Figure 7. (a) Temporal trends of autumn ET 0 for different elevations from 1964 to 2013 and (b) spatial pattern of MMK test in Heilongjiang Province, northeast China.Down ** indicates the decreasing trend at the 0.01 significance level, Down * indicates the decreasing trend at 0.05 significance level, and not significant means no significantly decreasing trend.(Green thick and solid line is the linear trend of ET 0 at 0-200 m, the Red thick and solid line is the linear trend of ET 0 at 200-400 m, Blue thick and solid line is the linear trend of ET 0 at >400 m).

Figure 8 .
Figure 8.(a) Temporal trends of winter ET 0 for different elevations from 1964 to 2013 and (b) spatial pattern of MMK test in Heilongjiang Province, northeast China.Down ** indicates the decreasing trend at the 0.01 significance level, Down * indicates the decreasing trend at 0.05 significance level, and not significant means no significantly decreasing trend.(Green thick and solid line is the linear trend of ET 0 at 0-200 m, the Red thick and solid line is the linear trend of ET 0 at 200-400 m, Blue thick and solid line is the linear trend of ET 0 at >400 m).

Figure 9 .
Figure 9. Trend of (a) winter maximum temperature (Tmax), and (b) winter minimum temperature (Tmin) at different elevations from 1964 to 2013 in Heilongjiang Province, northeast China.As the change slopes show, winter Tmin increased much more than winter Tmax in China's northernmost region.(Green thick and solid line is the linear trend of ET0 at 0-200 m, the Red thick and solid line is the linear trend of ET0 at 200-400 m, Blue thick and solid line is the linear trend of ET0 at >400 m).

Figure 10 .
Figure 10.The wavelet coherence (a) between winter ET0 and Tmax and (b) between winter ET0 and Tmin in Heilongjiang Province.A significance level of 5% against red noise is exhibited as a thick contour, and the relative phase relationship is denoted as arrows (with anti-phase pointing left and in-phase pointing right).

Figure 9 .
Figure 9. Trend of (a) winter maximum temperature (T max ), and (b) winter minimum temperature (T min ) at different elevations from 1964 to 2013 in Heilongjiang Province, northeast China.As the change slopes show, winter T min increased much more than winter T max in China's northernmost region.(Green thick and solid line is the linear trend of ET 0 at 0-200 m, the Red thick and solid line is the linear trend of ET 0 at 200-400 m, Blue thick and solid line is the linear trend of ET 0 at >400 m).

Figure 9 .
Figure 9. Trend of (a) winter maximum temperature (Tmax), and (b) winter minimum temperature (Tmin) at different elevations from 1964 to 2013 in Heilongjiang Province, northeast China.As the change slopes show, winter Tmin increased much more than winter Tmax in China's northernmost region.(Green thick and solid line is the linear trend of ET0 at 0-200 m, the Red thick and solid line is the linear trend of ET0 at 200-400 m, Blue thick and solid line is the linear trend of ET0 at >400 m).

Figure 10 .
Figure 10.The wavelet coherence (a) between winter ET0 and Tmax and (b) between winter ET0 and Tmin in Heilongjiang Province.A significance level of 5% against red noise is exhibited as a thick contour, and the relative phase relationship is denoted as arrows (with anti-phase pointing left and in-phase pointing right).

Figure 10 .
Figure 10.The wavelet coherence (a) between winter ET 0 and T max and (b) between winter ET 0 and T min in Heilongjiang Province.A significance level of 5% against red noise is exhibited as a thick contour, and the relative phase relationship is denoted as arrows (with anti-phase pointing left and in-phase pointing right).

Table 1 .
Geographical location and elevation of the 27 meteorological stations used in this study for estimating reference evapotranspiration for Heilongjiang Province, northeast China.

Table 1 .
Geographical location and elevation of the 27 meteorological stations used in this study for estimating reference evapotranspiration for Heilongjiang Province, northeast China.

Table 2 .
Average ET 0 (±std) of different seasons at different elevations in Heilongjiang Province, northeast China.

Table 2 .
Average ET0 (±std) of different seasons at different elevations in Heilongjiang Province, northeast China.

Table 3 .
Seasonal and interannual variation of ET0 from 1964 to 2013 at different elevations in Heilongjiang Province, Northeast China.

Table 3 .
Seasonal and interannual variation of ET0 from 1964 to 2013 at different elevations in Heilongjiang Province, Northeast China.

Table 3 .
Seasonal and interannual variation of ET 0 from 1964 to 2013 at different elevations in Heilongjiang Province, Northeast China.

Table 4 .
Average wind speed (WS), maximum and minimum temperature (T max , T min) ), sunshine duration (SD), and relative humidity (RH) during the four seasons at different elevations in Heilongjiang Province, northeast China.

Table 5 .
Decadal changes in wind speed (WS), maximum and minimum temperature (T max , T min ), subnshine duration (SD), and relative humidity (RH) from 1964 to 2013 in Heilongjiang Province, northeast China.

Table 6 .
Correlation coefficients between ET 0 and meteorological parameters, including wind speed (WS), maximum and minimum temperature (T max , T min ), sunshine duration (SD), and relative humidity (RH) from 1964 to 2013 at different elevations in Heilongjiang Province, northeast China.

Table 7 .
Seasonal and annual average sensitivity coefficients for individual climatic factors, including wind speed (WS), maximum and minimum temperature (T max , T min ), sunshine duration (SD), and relative humidity (RH) to ET 0 from 1964 to 2013 at different elevations in Heilongjiang Province, northeast China.
Note: Numbers in bold indicate the greatest sensitivity coefficients.