Spatiotemporal Variations in Reference Evapotranspiration and Its Contributing Climatic Variables at Various Spatial Scales across China for 1984–2019

: Reference evapotranspiration (ET 0 ) is of great signiﬁcance in studies of hydrological cycle, agricultural water resources, and hydrometeorology. The present study collected daily meteorological data at 536 meteorological stations in China from 1984 to 2019, calculated daily ET 0 using the FAO Penman–Monteith equation, analyzed the spatial distribution and temporal variation characteristics of ET 0 and meteorological variables at four different spatial scales (continental, regional, provincial, and local), and discussed the sensitivity of ET 0 to the meteorological variables and the contribution rates of the meteorological variables to the ET 0 variations. The results showed that ET 0 increased at 406 out of the 536 stations (75.7%), with the trends being signiﬁcant at 65 stations at the 5% signiﬁcance level, and 147 at the 1% signiﬁcance level. The slope calculated using Sen’s method and linear trend method showed that the annual ET 0 at the continental scale increased by approximately 12 mm/decade. Most of the stations showed decreasing trends in relative humidity (H m ), sunshine duration (S D ), and wind speed at 2 m height (U 2 ) while increasing trends in the maximum air temperature (T max ) and minimum air temperature (T min ). ET 0 was most sensitive to H m (sensitivity coefﬁcient, S t = − 0.66), followed by T max ( S t = 0.29), S D ( S t = 0.18), U 2 ( S t = 0.16), and T min ( S t = 0.07). Most of the stations showed increasing trends in S t for H m (56.16%), T max (72.95%), T min (87.31%), and U 2 (90.49%), and decreasing trends for S D (69.78%). The variations in H m , T max , and T min increased the ET 0 at most of the stations (82.28%, 98.13%, and 69.03%, respectively). The variations in S D and U 2 decreased ET 0 at most of the stations (66.04% and 56.34%, respectively). Some ET 0 characteristics in a few regions can be well described using a single spatial scale. However, most regions exhibited signiﬁcantly different ET 0 characteristics across spatial scales. The results of this project can provide reference for hydrological analysis and agricultural water management under climate change conditions and provide data and information for other hydrology-related applications.


Introduction
Evapotranspiration (ET) is known as a key parameter reflecting the variations in the climate and is a major component of the meteorological and hydrological cycles [1][2][3][4][5][6][7][8][9].With the risk of heatwave events and their persistence intensifying in recent past and increasing faster, the ET, precipitation, and other climatic variables have changed over the diverse regions [10,11].In many hydrodynamic and water quality modeling system, the ET can be directly used as an input variable, so it is important to better understand the spatiotemporal features of ET and the contributing factors in the context of ongoing climate change.Furthermore, a better understanding of ET is beneficial for analyzing the hydrodynamics and water quality, especially in terms of explaining the mechanisms for the variations in water quality, such as the variations in the balance of salinity and variations in contaminant concentrations.ET can be estimated by the reference evapotranspiration (ET 0 ).It has been widely accepted that ET 0 is an important factor in water balance and conversion and is thus its estimation and forecasting is important in water resources management [12][13][14][15][16][17].For example, estimating crop water demand is a necessary step in irrigation scheduling and the design of agricultural water conservancy projects, and the estimation of crop water demand is mainly based on the calculation of ET 0 [18,19].A clear understanding of the spatiotemporal variations in ET 0 is of importance for better planning and management of water resources.Contrary to intuitive expectation, ET 0 has been reported to decrease in many regions with increasing temperature [20][21][22], and this controversial behavior is known as the "evaporation paradox" [23][24][25].The "evaporation paradox" actually provides a straightforward indication that the ET 0 variations are products of changes in multiple variables rather than any one variable, and it has attracted many researchers to investigate the spatiotemporal variations in ET 0 and its contributing climatic factors.
Tang et al. [18] investigated ET 0 and its contributing meteorological variables in the Haihe River Basin of northern China and found that ET 0 exhibited a decreasing trend with a rate of −1 mm yr −2 ; the air temperature produced an increase in the derivative of ET over time (dET 0 /dt), while the net radiation and wind speed produced a decrease in (dET 0 /dt).Liu and McVicar [24] studied the potential evaporation and its contributing climatic factors in the Yellow River Basin in China and found that the potential evaporation in the water yielding region showed an increasing trend, while that in the water consuming region exhibited a decreasing trend during 1961-2010.Fully physics-based models should be employed to explicitly describe the contributing effects of the climatic factors on potential evaporation.Wang et al. [26] studied the spatial and temporal variations in ET 0 and the meteorological variables in China for the period 1961-2013 using a 0.5 • × 0.5 • gridded dataset; the study showed that there was a significant spatiotemporal variability in ET 0 , and the overall annual ET 0 declined at a rate of 6.84 mm/decade.The variation in ET 0 was mostly attributed to the wind speed (U w ), followed by the maximum air temperature (T max ), sunshine duration (S D ), relative humidity (H m ), and minimum air temperature (T min ).Jiang et al. [27] found that during the growing season in Southwest China ET 0 decreased at a rate of 10.25 mm/decade for the period 1961-1996 and increased at a rate of 8.12 mm/decade for the period 1997-2016; ET 0 was most sensitive to H m , followed by S D , T max , T min , and U w .Yan et al. [7] investigated the importance of the meteorological variables in predicting ET 0 at eight stations in China; the results indicated that the ET 0 in the arid region was mostly influenced by U w , followed by S D and H m , and ET 0 in the humid region was mostly influenced by S D , followed by H m and U w .
Although the existing studies have provided valuable understanding of the spatial and temporal variability in ET 0 and the relevant meteorological variables in China, a major drawback of all these works is that they have only examined the characteristics at a single or two spatial scales [28][29][30][31].The resulting findings of these previous studies cannot provide sufficient information for some large water conservancy projects that need to consider both global and local factors, and it is still unclear whether the global values derived from the data at discrete observation points can satisfactorily represent the local characteristics due to the lack of comparisons between observations across spatial scales.
Therefore, in contrast to the existing studies of ET 0 in China, there is a lack of comprehensive datasets of ET 0 distribution at different spatial scales across China.The primary purpose of the present study was to investigate the variability in ET 0 and its contributing climatic variables at various spatial scales.A total of four different spatial scales were considered, including continental, regional, provincial, and local.To specify, the present study was conducted toward achieving the following objectives: (1) to provide a new dataset of daily ET 0 at different spatial scales in China for 1984-2019, calculated using the FAO Penman-Monteith equation; this comprehensive dataset and the resulting findings do not currently exist in the literature and are important for projects that require both global and local considerations; (2) to analyze the spatial distribution and temporal variation characteristics of ET 0 and meteorological variables at different spatial scales adopting the Water 2022, 14, 2502 3 of 20 non-parametric Mann-Kendall test method; and (3) to quantitatively analyze the sensitivity of ET 0 to the meteorological variables and the contribution rates of the meteorological variables to ET 0 variations using the partial derivative sensitivity coefficients method, and determine the key climatic variables affecting ET 0 at different spatial scales.

Study Area and Climate Dataset
Located in the east of Asia, China is the third largest country in the world in land area, with an estimated land area of 9.6 million square kilometers and a population of over 1.4 billion.Its terrain is generally high in the West and low in the East.Most of the areas belong to temperate and subtropical monsoon climate, with rich and diverse climate types and geographical landscape.The issues surrounding in China are severe partially because it has uneven spatial and temporal distribution of water resources.The precipitation generally decreases from the southeast coast to northwest inland.The cultivated land in the Yangtze River Basin and south of the Yangtze River only accounts for 36% of China's total, while water resources account for 80% of China's total.In the Yellow River, Huaihe River and Haihe River basins, water resources only account for 8% of China's total, while arable land accounts for 40% of China's total.The annual and interannual variation of precipitation and runoff is large.These characteristics of water resources cause frequent floods and droughts in China, and agricultural production is also relatively unstable.
In the present study, a total of four spatial scales were considered: continental, regional, provincial, and local.The continental scale was equivalent to the national scale, which covered mainland China.At the regional scale, the study area has been divided into 7 subregions by the following geographic locations: North Eastern Region (NER), North China Region (NCR), East China Region (ECR), South China Region (SCR), Central China Region (CCR), North Western Region (NWR), and South Western Region (SWR).A large portion of the water resources activities in China are governed by provincial water management authorities, so the provincial scale was also considered.The local scale considered the stations that have sufficient quality data.
The daily climatic data employed in this study included daily H m , S D , T max , T min , and U w at 10 m height (U 10 ).The data of U 10 has been converted to U 2 (U w at 2 m height) using the transformation formula provided by Allen et al. [32].The data were primarily obtained from the China Meteorological Administration.The quality of the dataset has been carefully controlled, and the dataset has been widely used in previous studies [9,24,28].There were 665 stations in the original dataset, but 129 of the stations were excluded because they had extensive missing data; finally, 536 stations were selected.The study area and selected meteorological stations are presented in Figure 1.The subregions, the provinces, and the corresponding numbers of meteorological stations are provided in Table 1.

ET0 Calculation
The present study employed the FAO Penman-Monteith (P-M) equation to calculate the daily ET0.In March 1990, the Food and Agriculture Organization of the United Nations (FAO) recommended the use of this formula to calculate ET0 at the international symposium on crop evapotranspiration calculation methods held in Rome, Italy [19].The formula is the most widely used method to calculate ET0 in the world.The P-M formula, physically based on both the energy balance and aerodynamics, does not need a special regional calibration function, and it can calculate ET0 by using general meteorological data, which has high practical value and accuracy [33].The equation can be written as where 0 ET is the reference evapotranspiration (mm day -1 );  is the slope vapor pres- sure curve (kPa °C−1 ); n R is the net radiation at the crop surface (MJ m −2 day −1 ); G is

ET 0 Calculation
The present study employed the FAO Penman-Monteith (P-M) equation to calculate the daily ET 0 .In March 1990, the Food and Agriculture Organization of the United Nations (FAO) recommended the use of this formula to calculate ET 0 at the international symposium on crop evapotranspiration calculation methods held in Rome, Italy [19].The formula is the most widely used method to calculate ET 0 in the world.The P-M formula, physically based on both the energy balance and aerodynamics, does not need a special regional calibration function, and it can calculate ET 0 by using general meteorological data, which has high practical value and accuracy [33].The equation can be written as where ET 0 is the reference evapotranspiration (mm day -1 ); ∆ is the slope vapor pressure curve (kPa • C −1 ); R n is the net radiation at the crop surface (MJ m −2 day −1 ); G is the soil heat flux density (MJ m −2 day −1 ); γ is the psychrometric constant (kPa • C −1 ); T is the mean daily air temperature at 2 m height ( • C); U 2 is the U w at the 2m height (m/s); e s is the saturation vapor pressure (kPa); e a is the actual vapor pressure (kPa).The intermediate variables can be calculated from daily H m , S D , T max , T min , and U 2 , following the equations provided in the FAO56 report [32].

Trend Analysis
The present study used the Mann-Kendall trend test approach to detect significant trends in ET 0 and its contributing climatic variables.The Mann-Kendall (MK) test [34][35][36] is a nonparametric (distribution-free) test, and it was employed in this study because it has been widely used in assessing whether there is a monotonic upward or downward trend of hydrological variables over time [27,37].Its superiority exists in the ability of detecting the significance of the long time-series trend, no matter the sample series follow a certain type of distribution [38].In the MK test, S is a statistical value used to determine whether values obtained later in time tend to be larger than values obtained earlier, as: where n is the number of observations and sign(x i − x j ) can be expressed as: The variance of S can be computed as follows: where l is the number of tied groups, and t k is the number of observations in the tth group.
The MK test statistic, Z S , can be expressed as: Positive values of Z S indicate that the variables of interest tend to increase with time, and negative values of Z S indicate that the variables of interest tend to decrease with time.The null hypothesis of no trend is rejected when |Z S | > 1.96 at the α = 0.05 level and |Z S | > 2.576 at the α = 0.01 level.
To estimate the slope of the linear trend, the time series of a variable (e.g., ET 0 ) were fitted into a linear regression equation, such as where T L is slope of the liner trend.The slope of trend can also be estimated using a nonparametric procedure developed by Sen [39].In the procedures, the slopes of N pairs of observations are first estimated by where x j is the value at time j, while x k is the value at time k (j > k).Since there is only one datum in each time period in the present study, N can be computed as N = n(n−1) 2 and n is the number of time steps.The values of these slopes are then ranked from smallest to largest, and the median value is known as the Sen's estimator of slope.When N is odd, Sen's estimator of slope is Q (N+1)/2 ; when N is even, Sen's estimator of slope is . To verify whether the median slope is statistically different from zero, the confidence interval for Sen's estimator of slope, for a given level of significance α, should be computed [34].The confidence interval can be computed as follows: where Var(S) is computed by (4), and Z 1−α/2 is obtained from the standard normal distribution table.The approximate normal theoretical lower and upper confidence limit are respectively the M 1 th largest and the (M 2 +1)th largest of value of Q. M 1 and M 2 are defined by the Sen's estimator of slope is statistically different from zero if the two limits Q(M 1 ) and Q(M 2 ) have the same sign.

Sensitivity Coefficient and Contribution Rate
The sensitivity coefficient, St i , that quantifies the sensitivity of ET 0 to a climatic variable, c i , can be expressed as [27]: If St i > 0, ET 0 and c i will increase or decrease simultaneously, otherwise, they have a reverse relation.The larger the |St i |, the greater the effect of a climatic variable has on ET 0 variation.The sensitivity coefficient method is more convenient and accurate compared to other sensitivity analysis approaches [40], and it has been widely used in hydro-meteorological study.
The contribution rate of c i , C bi , to ET 0 variation, can be expressed as [27]: where n is the number of observations, Td i is the trend of the variable c i , and |Ac i | is the absolute value of the averaged c i .If the C bi is positive (negative), it denotes the variable had a positive (negative) contribution to ET 0 variation.The method can clearly reveal how the climatic variables contribute to ET 0 .

Spatiotemporal Variation of ET 0
Figure 2 shows the descriptive statistics of daily ET 0 at the 536 stations.The spatial distribution plots, which were obtained by interpolating the results using the Kriging interpolation approach, are provided primarily for the purpose of visualization.The mean daily ET 0 ranged from around 1 mm (mainly in the northeast) to around 4 mm (such as in the southern edge).The minimum daily ET 0 in the northern part was close to 0 mm, partially because the winter temperature in the north was quite low.The minimum daily ET 0 in the southern part was around 1 mm.Although the temperature increased from north to south, the maximum daily ET 0 showed a reverse trend, which may partially be attributed to the higher U 2 in the northern part, and revealed that ET 0 is the combined result of various climatic variables.The median values were close to the mean values.Generally, the standard deviation and variance decreased from north to south, and the kurtosis and skewness increased from northeast to southwest.
The results of the M-K trend analysis are shown in Figure 3. ET 0 increased at 406 of the 536 stations (75.7%), with the trends being significant at 65 stations at the 5% significance level, and 147 at the 1% significance level.ET 0 decreased at 128 out of the 536 stations (23.9%), with the trends being significant at 13 stations at the 5% significance level, and 13 at the 1% significance level.In general, most of the stations exhibited an increasing ET 0 , and the stations showing significant increasing trends were primarily located in NM, NX, GS, SC, CQ, and YN.The averaged slope calculated using Sen's method and linear trend method were 0.0033 and 0.0034, respectively.Therefore, the annual ET 0 at the continental scale increased by approximately 12 mm/decade in China.The results of the M-K trend analysis are shown in Figure 3. ET0 increased at 406 of the 536 stations (75.7%), with the trends being significant at 65 stations at the 5% significance level, and 147 at the 1% significance level.ET0 decreased at 128 out of the 536 stations (23.9%), with the trends being significant at 13 stations at the 5% significance level, and 13 at the 1% significance level.In general, most of the stations exhibited an increasing ET0, and the stations showing significant increasing trends were primarily located in NM, NX, method were 0.0033 and 0.0034, respectively.Therefore, the annual ET0 at the continental scale increased by approximately 12 mm/decade in China.

Spatiotemporal Variation of Climatic Factors
The descriptive statistics of the climatic variables were shown in Figure 4. Based on the mean data, Hm increased from the northwest (less than 40%) to the southeast (more than 80%); SD decreased from west (as high as 9 h) to east (as low as 3 h); Tmax increased from north (lower than 5 °C) to south (close to 30 °C); Tmin also increased from north (lower than −10 °C) to south (higher than 20 °C); and U2 decreased from north (close to 4 m/s) to south (as low as 1 m/s).The median values were close to the mean values.The standard deviations ranged from around 19% in the northwest to around 8% in southeast.Conversely, the standard deviations for SD increased from the northwest (approximately 4.5 h) to the southeast (approximately 2.5 h).The standard deviations for Tmax, Tmin, and U2 all decreased from north to east.These standard deviations can quantify how accurate it was to estimate climatic variables at a smaller spatial scale using characteristics at a larger spatial scale.

Spatiotemporal Variation of Climatic Factors
The descriptive statistics of the climatic variables were shown in Figure 4. Based on the mean data, H m increased from the northwest (less than 40%) to the southeast (more than 80%); S D decreased from west (as high as 9 h) to east (as low as 3 h); T max increased from north (lower than 5 • C) to south (close to 30 • C); T min also increased from north (lower than −10 • C) to south (higher than 20 • C); and U 2 decreased from north (close to 4 m/s) to south (as low as 1 m/s).The median values were close to the mean values.The standard deviations ranged from around 19% in the northwest to around 8% in southeast.Conversely, the standard deviations for S D increased from the northwest (approximately 4.5 h) to the southeast (approximately 2.5 h).The standard deviations for T max , T min , and U 2 all decreased from north to east.These standard deviations can quantify how accurate it was to estimate climatic variables at a smaller spatial scale using characteristics at a larger spatial scale.
Figure 5 shows the results of trend analysis for the climatic variables.H m decreased at 443 of the 536 stations (82.6%), and the averaged slopes calculated using Sen's method and linear trend method were −0.089 and −0.086, respectively.The majority of the stations (67.2%) exhibited decreasing S D , with an average trend slope of −0.006 (Sen's = −0.006;and linear = −0.006).Most of the stations showed increasing T max (532 stations; 99.3%), and the averaged slope was 0.039 (Sen's = 0.040; and linear = 0.039).Similarly, most of the stations (512 stations; 95.5%) had an increasing T min with the averaged slope being 0.042 (Sen's = 0.042; and linear = 0.042).The U 2 showed an overall decreasing trend, as 324 stations (60.4%) had a declining U 2 , and the averaged Sen's slope and linear slope were −0.007 and −0.005, respectively.Figure 5 shows the results of trend analysis for the climatic variables.Hm decreased at 443 of the 536 stations (82.6%), and the averaged slopes calculated using Sen's method and linear trend method were −0.089 and −0.086, respectively.The majority of the stations (67.2%) exhibited decreasing SD, with an average trend slope of −0.006 (Sen's = −0.006;and linear = −0.006).Most of the stations showed increasing Tmax (532 stations; 99.3%), and the averaged slope was 0.039 (Sen's = 0.040; and linear = 0.039).Similarly, most of the stations (512 stations; 95.5%) had an increasing Tmin with the averaged slope being 0.042 (Sen's = 0.042; and linear = 0.042).The U2 showed an overall decreasing trend, as 324 stations (60.4%) had a declining U2, and the averaged Sen's slope and linear slope were −0.007 and −0.005, respectively.

Variability of ET0 at Various Spatial Scales
The descriptive statistics of daily ET0 for each station (at the local have been presented in Figure 2 and the mean values were averaged to obtain the mean ET0 at the provincial, regional, and continental scales.This averaging approach was used for upscaling, instead of other approaches based on interpolations, to avoid the large uncertainties

Variability of ET 0 at Various Spatial Scales
The descriptive statistics of daily ET 0 for each station (at the local scale) have been presented in Figure 2 and the mean values were averaged to obtain the mean ET 0 at the provincial, regional, and continental scales.This averaging approach was used for upscal-ing, instead of other approaches based on interpolations, to avoid the large uncertainties caused by spatial interpolations.At the continental scale, the daily ET 0 across China for 1984-2019 was 2.72 mm/day.However, the standard deviation of the mean daily ET 0 at different stations was 0.40 mm/day, which was quite high, indicating that the ET 0 characteristics found at the continental scale cannot provide meaningful insight for local water management.Wang et al. [26] investigated the spatiotemporal variability of ET 0 during 1961-2013 using a gridded dataset and also revealed that the ET 0 in China differed significantly.Thus, the ET 0 variation characteristics across different spatial scales should be investigated.At the regional scale, the averaged daily ET 0 was the highest in South China (SCR), which was 3.25 mm/day, followed by NWR (2.80 mm/day), NCR (2.78 mm/day), ECR (2.76 mm/day), SWR (2.69 mm/day), and CCR (2.68 mm/day).The averaged ET 0 was the lowest in North East (NER), which was 2.24 mm/day.
Figure 6 illustrates the time series of annual mean daily ET 0 at the provincial scale.In the area NER, province LN generally had the highest ET 0 (mean ET 0 =2.51 mm/d), while HL had the lowest ET 0 (mean ET 0 = 2.06 mm/d).Province TJ had the highest ET 0 (mean ET 0 = 3.01 mm/d) in the subregion NCR, and NM had the lowest ET 0 (mean ET 0 = 2.61 mm/d).In ECR, the ET 0 was the highest in province SH (mean ET 0 = 2.83 mm/d) and the lowest in AH (mean ET 0 = 2.60 mm/d).In the subregion SCR, ET 0 ranged from 2.96 mm/d in GX to 3.60 mm/d in HI.In the subregion CCR, the difference in ET 0 was relatively small, with HA being the province with the highest ET 0 (mean ET 0 = 2.80 mm/d) and HB being the province with the lowest ET 0 (mean ET 0 = 2.61 mm/d).The province that had the maximum mean ET 0 in the subregion NWR was XJ (mean ET 0 =3.06 mm/d), and the one that had the minimum mean ET 0 was QH.In the subregion SWR, the mean ET 0 ranged from 2.47mm/d in province CQ to 3.02 mm/d in YN.In some areas, such as SCR, the differences at different locations were relatively small, and thus, it is reasonable to estimate ET 0 characteristics at a smaller spatial scale using ET 0 characteristics at a larger spatial scale.However, most of the regions exhibited significant differences of ET 0 at different locations; thus, the regional observations cannot well represent the provincial characteristics.
The trends in ET 0 variations at the local scale have been presented in Figure 3, and the trends at the provincial and regional scales can be analyzed based upon these trends.The statistics regarding the trends in ET 0 at the regional scale are summarized in Figure 7.In the subregion NCR, 34% of the stations showed a trend of "increase-significant α = 0.01", followed by "increase" (28%) and "increase-significant α = 0.05" (18%).In the area ECR, 39% of the stations showed a trend of "increase", followed by "decrease" (26%) and "increasesignificant α = 0.01" (20%).39% of the stations in the subregion SCR showed a trend of " increase", while 25% showed a trend of "decrease" and 22% showed a trend of "increasesignificant α = 0.05".In the area SWR, a large portion of the stations showed a trend of "increase" (45%) and "increase-significant α = 0.01" (31%), and 12% of the stations showed a trend of "decrease".There were 39%, 27%, and 13% of the stations in the subregion NWR that showed a trend of "increase-significant α = 0.01", "increase", and "decrease" (13%), respectively.The majority of the stations in CCR showed an increasing trend (34% "increase" and 30% "increase-significant α = 0.01" (30%)", and there were also 26% that showed the "decrease" trend.In the subregion NER, most of the stations (44%) showed a trend of "increase"; however, there were also 26% of the stations showed a trend of "decrease" (26%), exhibiting a significant diverse; the portion of stations that showed a trend of "increasesignificant α = 0.01" was 16%.At the continental scale, the annual ET 0 increased by approximately 12 mm/decade for 1984-2019.Most of the stations showed increasing trends, with 36% of the stations exhibiting a trend of "increase" and 27% exhibiting "increasesignificant α = 0.01"(27%); however, there were also 19% of the stations that showed a trend of "decrease".Zhang et al. [31] studied evapotranspiration (ET) in China for 1948-2018 and found that the ET in most regions (covered by 89.5% pixels) showed increasing trends, which was close to the current observation.Zhang et al. [31] reported that the fluctuation of ET was relatively larger in the north and west, and relatively weaker in the south and east.
The present study found the same characteristic, but it reported the standard deviations and variances of ET 0 variations at each station to quantitatively describe the fluctuations (Figure 2).The trends in ET0 variations at the local scale have been presented in Figure 3, the trends at the provincial and regional scales can be analyzed based upon these tre The statistics regarding the trends in ET0 at the regional scale are summarized in Figu In the subregion NCR, 34% of the stations showed a trend of "increase-significant 0.01", followed by "increase" (28%) and "increase-significant α = 0.05" (18%).In the ECR, 39% of the stations showed a trend of "increase", followed by "decrease" (26%) "increase-significant α = 0.01" (20%).39% of the stations in the subregion SCR show trend of " increase", while 25% showed a trend of "decrease" and 22% showed a tren

Variability of Climatic Variables at Various Spatial Scales
The statistics at the regional and continental (national) scales are illustrated in Figure 8, in which the trends and M-K Zs values are also indicated.The lower lines of the boxes indicated the lower quartiles, which were the medians of the lower half of the datasets, the medium lines correspond to the middle values of the entire datasets, and the upper lines indicate the upper quartiles, which were the medians of the upper half of the datasets.These lines of boxes can also provide insight for the variability of climatic variables at the provincial scale.
As can be seen, the Hm in ECR, SCR, and CCR was relatively higher, whereas that in NCR and NWR was relatively lower.All of the subregions showed decreasing trends in Hm, with NCR having the most significant trend (Zs = −2.39)and SCR having the least significant trend rate (Zs = −0.78).The subregions NCR, NWR, and NER had relatively greater values of SD than the other subregions.SD was declining in all of the areas, and the absolute value of Zs was the largest in NCR (Zs = −1.52)and smallest in SWR (Zs = −0.34).Yin et al. [15] has studied the potential evapotranspiration (PET) and its determining factors in China for 1971-2008 and has also found a decreasing trend of SD.Tmax in SCR, CCR, and

Variability of Climatic Variables at Various Spatial Scales
The statistics at the regional and continental (national) scales are illustrated in Figure 8, in which the trends and M-K Z s values are also indicated.The lower lines of the boxes indicated the lower quartiles, which were the medians of the lower half of the datasets, the medium lines correspond to the middle values of the entire datasets, and the upper lines indicate the upper quartiles, which were the medians of the upper half of the datasets.These lines of boxes can also provide insight for the variability of climatic variables at the provincial scale.

Sensitivity and Contribution
The variations in ET0 are the combined effects of the climatic variables; thus, it is important to quantify the sensitivity of ET0 to these factors for a better understanding of the driving mechanisms.The mean values of the sensitivity coefficients, the corresponding As can be seen, the H m in ECR, SCR, and CCR was relatively higher, whereas that in NCR and NWR was relatively lower.All of the subregions showed decreasing trends in H m , with NCR having the most significant trend (Z s = and SCR having the least significant trend rate (Z s = −0.78).The subregions NCR, NWR, and NER had relatively greater values of S D than the other subregions.S D was declining in all of the areas, and the absolute value of Z s was the largest in NCR (Z s = −1.52)and smallest in SWR (Z s = −0.34).Yin et al. [15] has studied the potential evapotranspiration (PET) and its determining factors in China for 1971-2008 and has also found a decreasing trend of SD.T max in SCR, CCR, and ECR was relatively higher than others in NCR and NER.All of the areas showed increasing trends, with the Z s values ranging from 2.18 in NER to 4.13 in ECR.Similarly, T min was also increasing in all of the subregions.T min in ECR, SCR, and NER was relatively higher, whereas that in NCR and NWR was relatively lower.The Z s values for T min ranged from 2.16 in NER to 4.59 in CCR.In contrast to T max and T min , U 2 showed a decreasing trend, with SCR having the most significant trend rate (Z s = −0.04)and ECR having the least significant (Z s = −2.68).The subregions with relatively high U 2 included NCR, ECR, and NER, and those subregions with low U 2 included SWR, SCR, and CCR.Yin et al. [15] has also found that U 2 showed a decreasing trend.In summary, the results at the continental scale can well represent the trend directions at the regional scale.In some regions, the findings at the regional scale can satisfactorily represent the characteristics at the provincial and local scales, such as H m in SCR and T max in SWR.However, in most regions, both the statistics and trend rates exhibited substantial variations across spatial scales, such as S D in ECR and T min in NER.

Sensitivity and Contribution
The variations in ET 0 are the combined effects of the climatic variables; thus, it is important to quantify the sensitivity of ET 0 to these factors for a better understanding of the driving mechanisms.The mean values of the sensitivity coefficients, the corresponding trends (Z s values), and trend rates (Sen's slope) at the local scale are presented in Figure 9.
Figure 9 shows that ET 0 at most stations was negatively correlated to the variations in H m , while ET 0 at most stations was positively correlated to the variations in the other climatic variables.At the continental scale, ET 0 was most sensitive to H m (average S t = −0.66),followed by T max (S t = 0.29), S D (S t = 0.18), U 2 (S t = 0.16), and T min (S t = 0.07).Yin et al. (2010) also revealed that H m was the most sensitive variable for PET.
Consistent with the observations of Fan and Thomas [29], who studied ET 0 in China for the period of 1960-2011, the importance of the climatic variables to ET 0 were not stable.Most of the stations showed increasing trends in S t for H m (56.16%), with a Z s value of 0.37 and a slope of 0.00092 at the continental scale.Fan and Thomas [29] has also found that the importance of H m increased in China.In terms of S D , 69.78% of the stations showed decreasing trends, and the Z s value and slope at the continental scale were −1.15 and −0.00027, respectively.Liu et al. [25] found that the sensitivity of T max exhibited declining trends while that of T min showed increasing trends in China for 1960-2007.However, based on the present study, which covered a longer time period, the majority of the stations showed increasing trends in S t for T max (72.95%) and T min (87.31).The Z s values at the continental scale corresponding to T max and T min were 1.25 and 2.16, respectively, and the slopes at the continental scale were 0.00040 and 0.00032, respectively.Interestingly, 90.49% of the stations showed increasing trends in S t for U 2 , implying that U 2 was becoming more important in affecting ET 0 .Liu et al. [25] also found that the sensitivity of U 2 showed an increasing trend in China for 1960-2007.The Z s value and slope at the continental scale corresponding to U 2 was 2.38 and 0.00106, respectively.The contribution of each climatic variable to the variation in ET0 was quantified in Figure 10, in which the spatial distribution of Zs for ET0 was also shown for a better interpretation of the climatic effects.As can be seen from Figure 10, the variations in Hm increased ET0 at most of the stations (441 stations; 82.28%), and the Cb at the continental scale was 2.83%.Similarly, the variations in Tmax and Tmin increase ET0 at 526 stations (98.13%) and 370 stations (69.03%), respectively, and the Cb values at the continental scale were 2.32% and 0.07%, respectively.However, for most stations, the variations in SD and U2 decreased ET0 at most of the stations (354 stations and 302 stations, respectively), and the Cb values at the continental scale were −0.67% and −0.70%, respectively.Yao et al. [30] studied the PET across China for 1982-2010 and found that PET was mostly attributed to The contribution of each climatic variable to the variation in ET 0 was quantified in Figure 10, in which the spatial distribution of Z s for ET 0 was also shown for a better interpretation of the climatic effects.As can be seen from Figure 10, the variations in H m increased ET 0 at most of the stations (441 stations; 82.28%), and the C b at the continental scale was 2.83%.Similarly, the variations in T max and T min increase ET 0 at 526 stations (98.13%) and 370 stations (69.03%), respectively, and the C b values at the continental scale were 2.32% and 0.07%, respectively.However, for most stations, the variations in S D and U 2 decreased ET 0 at most of the stations (354 stations and 302 stations, respectively), and the C b values at the continental scale were −0.67% and −0.70%, respectively.Yao et al. [30] studied the PET across China for 1982-2010 and found that PET was mostly attributed to the incident solar radiation, followed by pressure deficit, air temperature, and wind speed, which were close to the current observations.Ullah et al. [40] and Hina et al. [41] investigated the climatic variables related to drought and possible cycles over Pakistan.The analysis of atmospheric circulation patterns revealed that large-scale changes in wind speed, air temperature, relative humidity, and so on are the likely drivers of droughts in the region, which also confirmed the current observations were reasonable.Sein et al. [42] assessed the spatiotemporal variation of summer monsoon precipitation and its potential drivers in Myanmar and found that the spatial trends for monthly and seasonal precipitation vary from station to station and region to region due to a fluctuated shift of climatic dynamics, which demonstrated the precipitation and ET 0 variability is associated with regional scale atmospheric circulation and some climatic variables.Based on the values of C b obtained in this study using the latest dataset, it can be concluded that the variations in ET 0 were mostly attributed to H m , T max , U 2 , S D , and T min .
Water 2022, 14, x FOR PEER REVIEW 20 of 23 the incident solar radiation, followed by pressure deficit, air temperature, and wind speed, which were close to the current observations.Ullah et al. [40] and Hina et al. [41] investigated the climatic variables related to drought and possible cycles over Pakistan.The analysis of atmospheric circulation patterns revealed that large-scale changes in wind speed, air temperature, relative humidity, and so on are the likely drivers of droughts in the region, which also confirmed the current observations were reasonable.Sein et al. [42] assessed the spatiotemporal variation of summer monsoon precipitation and its potential drivers in Myanmar and found that the spatial trends for monthly and seasonal precipitation vary from station to station and region to region due to a fluctuated shift of climatic dynamics, which demonstrated the precipitation and ET0 variability is associated with regional scale atmospheric circulation and some climatic variables.Based on the values of Cb obtained in this study using the latest dataset, it can be concluded that the variations in ET0 were mostly attributed to Hm, Tmax, U2, SD, and Tmin.

Conclusion
Daily ET0 at 536 meteorological stations in China for the period 1984-2019 were calculated, the trends in ET0 and its contributing climatic variables were analyzed, and the sensitivity of ET0 to the meteorological variables and the contribution rates of the meteorological variables to the ET0 variations were discussed.A key novelty of this study was that the variability of ET0 and the contributing meteorological variables at various spatial scales were investigated.Generally, ET0 decreased from south to north.The core of high ET0 was primarily located near the southern edge.Generally, the standard deviation and variance decreased from north to south, and the kurtosis and skewness increased from northeast to southwest.The ET0 increased at 406 of the 536 stations, and the averaged slope showed that the annual ET0 increased by approximately 12 mm/decade in China.Most of the stations showed decreasing trends in Hm (Zs = −1.74),SD (Zs = −0.87),and U2 (Zs = −1.07)and increasing trends in Tmax (Zs = 3.30) and Tmin (Zs = 3.85).
ET0 was negatively correlated to Hm, and the average sensitivity coefficient was −0.66.In terms of positively correlated variables, ET0 was most sensitive to Tmax (St = 0.29), followed by SD (St = 0.18), U2 (St = 0.16), and Tmin (St = 0.07).The St for Hm at most stations (56.16%) were increasing while the SD at most stations (69.78%) were decreasing.The majority of the stations showed increasing trends in St for Tmax (72.95%) and Tmin (87.31%).Interestingly, 90.49% of the stations showed increasing trends in St for U2, implying that U2 was becoming more important in affecting ET0.The variations in Hm, Tmax, and Tmin

Conclusions
Daily ET 0 at 536 meteorological stations in China for the period 1984-2019 were calculated, the trends in ET 0 and its contributing climatic variables were analyzed, and the sensitivity of ET 0 to the meteorological variables and the contribution rates of the meteorological variables to the ET 0 variations were discussed.A key novelty of this study was that the variability of ET 0 and the contributing meteorological variables at various spatial scales were investigated.Generally, ET 0 decreased from south to north.The core of high ET 0 was primarily located near the southern edge.Generally, the standard deviation and variance decreased from north to south, and the kurtosis and skewness increased from northeast to southwest.The ET 0 increased at 406 of the 536 stations, and the averaged slope showed that the annual ET 0 increased by approximately 12 mm/decade in China.Most of the stations showed decreasing trends in H m (Z s = −1.74),S D (Z s = −0.87),and U 2 (Z s = −1.07)and increasing trends in T max (Z s = 3.30) and T min (Z s = 3.85).
ET 0 was negatively correlated to H m , and the average sensitivity coefficient was −0.66.In terms of positively correlated variables, ET0 was most sensitive to T max (S t = 0.29), followed by S D (S t = 0.18), U 2 (S t = 0.16), and T min (S t = 0.07).The S t for H m at most stations (56.16%) were increasing while the S D at most stations (69.78%) were decreasing.The majority of the stations showed increasing trends in S t for T max (72.95%) and T min (87.31%).Interestingly, 90.49% of the stations showed increasing trends in S t for U 2 , implying that U 2 was becoming more important in affecting ET 0 .The variations in H m , T max , and T min increased the ET 0 at most of the stations (the mean contribution rate was 2.83%, 2.32%, and 0.07%, respectively).However, at most stations, the variations in S D and U 2 decreased ET 0 , and the mean C b values were −0.67% and −0.70%, respectively.
The present study reported an analysis of trends in ET 0 and its contributing climatic variables at four different spatial scales.Compared to analyses of a single or two spatial scales, the present analysis can better illustrate the variation trends and aid in future studies in water resources, agriculture, and hydrology in China.This study is also meaningful for modeling of hydrodynamics and water quality in inland and coastal water: first, the data can be directly used as inputs of the models; and second, the studied ET 0 characteristics can provide a better understanding of the mechanics and reasons for the variations in hydrodynamics and water quality.The present study found that the variations in ET 0 , climatic variables, sensitivity coefficients, and contribution rates were significant at some stations, so it is important to regularly update the corresponding forecasts.The significant variations also emphasize that it is important to develop sustainable water resources management and water-saving irrigation techniques under the condition of climate change.This study used the averaging approach for upscaling, it deserves a separate study to perform other upscaling techniques to study ET 0 characteristics at different scales.

Figure 1 .
Figure 1.The study area and locations of the meteorological stations (the dots).

Figure 1 .
Figure 1.The study area and locations of the meteorological stations (the dots).

Figure 4 .
Figure 4. Descriptive statistics of the daily climatic variables: (a) mean; (b) median; and (c) standard deviation.

Figure 4 .
Figure 4. Descriptive statistics of the daily climatic variables: (a) mean; (b) median; and (c) standard deviation.

Figure 5 .
Figure 5. Spatial distribution maps for the trend analysis of the climatic variables: (a) M-K ZS, (b) Sen's slope, and (c) linear trend.

Figure 5 .
Figure 5. Spatial distribution for the trend analysis of the climatic variables: (a) M-K Z S , (b) Sen's slope, and (c) linear trend.

Figure 8 .
Figure 8. Statistics of the climatic factors in different subregions.The lower line = lower quartile; the medium line = the middle values; the upper line = the upper quartile; the arrow = the trend; and the numbers = the values of Zs.

Figure 8 .
Figure 8. Statistics of the climatic factors in different subregions.The lower line = the lower quartile; the medium line = the middle values; the upper line = the upper quartile; the arrow = the trend; and the numbers = the values of Z s .

Figure 10 .
Figure 10.Spatial maps for the contribution analysis: (a) the M-K Zs value of ET0; (b) the Cb of RH; (c) the Cb of SD; (d) the Cb of Tmax; (e) the Cb of Tmin; and (f) the Cb of U2.

Figure 10 .
Figure 10.Spatial maps for the contribution analysis: (a) the M-K Zs value of ET 0 ; (b) the C b of RH; (c) the C b of S D ; (d) the C b of T max ; (e) the C b of T min ; and (f) the C b of U 2 .

Table 1 .
The study area, province, and corresponding number of stations.

Table 1 .
The study area, province, and corresponding number of stations.