Spatiotemporal Characteristics of Actual Evapotranspiration Changes and Their Climatic Causes in China

Abstract

As the main expenditure item in water balance, evapotranspiration has an important impact on the surface ecosystem. Assessing the impact of changes in meteorological elements on evapotranspiration is essential to identify the spatiotemporal heterogeneity of hydrographic responses to climate changes. Based on the actual evapotranspiration (ET_a) product (GPR-ET) generated by Gaussian process regression (GPR), as well as temperature and precipitation datasets, our study employed various statistical analysis methods, including geographic detector, the center of gravity migration model, spatial variation coefficients, and partial differential models, to investigate the spatiotemporal variation in ET_a in China from 2000 to 2018. The analysis covered future trends in ET_a changes and the contribution of meteorological factors. Our results showed that the ET_a in northwest China had stronger spatial heterogeneity and the mean value was generally lower than that in the southeast. But the center of gravity of ET_a was shifting towards the northwest. In most areas, the future trend was expected to be inconsistent with the current stage. ET_a in the regions of north and west was mainly driven by precipitation, while its increase in southeast China was largely attributed to temperature. In addition to spatial variations, the joint enhancement effect of temperature and precipitation on ET_a exists. According to the contribution analysis, precipitation contributed more to the change in ET_a than temperature. These findings have enhanced our comprehension of the contribution of climate variability to ET_a changes, providing scientific proof for the optimization apportion of future water resources.

1. Introduction

The transmission process of water to the atmosphere via soil evaporation and vegetation transpiration, collectively called evapotranspiration (ET), is crucial to the heat and mass exchange processes that take place in the Earth’s atmosphere [1,2]. Accurately calculating ET is of utmost importance for formulating reasonable irrigation systems and improving water resource utilization efficiency. The estimation models of ET can be divided into three categories: physical models, semi-physically based models, and black-box models. (1) Physical models establish ET formulas based on theoretical assumptions including energy conservation, complementary relationships, Budyko’s assumption, similarity of gradient fluxes, etc. The classical ET models include the Penman-Monteith model [3,4], surface energy balance systems (SEBS) [5], and Budyko–coupled water-energy balance equations [6], the Priestley–Taylor equation [7], etc. (2) Semi-physically based models have the advantages of both physical laws and mathematical statistics. While maintaining the basic physical laws, the strong approximation performance of machine learning algorithms is also adopted to compute the critical parameters and poorly described subroutines. Physically constrained machine learning models [8] and hybrid models [9] are two common methods among them. (3) The properties of the black-box models make them useful for mapping the relation among inputs and outputs even without a relevant physical basis. Its utilization in ET evaluation can effectively substitute the data-intensive empirical model with relatively poor adaptability. The commonly used methods include multilayer perceptron [10,11], fuzzy models [12,13], and neural networks [14], etc. At present, the Penman–Monteith formula recommended by the Food and Agriculture Organization of the United Nations (FAO) is a recognized standardized approach for counting reference crop evapotranspiration (ET_o), which is sufficiently accurate in terms of calculation and also a sound theoretical basis [15]. The Penman–Monteith formula has a strict physical foundation, comprehensively considering the effects of plant stomatal conductance and meteorological factors on ET; the formula can be applied directly and does not require parameter correction for different climatic conditions. It has been shown to be a high-precision method for estimating ET_o in different regions and environments [16,17]. The crop ET under standardized conditions is noted as ET_c. The crop coefficient method proposed by FAO indicates that multiplying the crop coefficient K_c with ET_o can acquire ET_c. Due to its clear mechanism, strong practicality, and reliable estimation accuracy, the crop coefficient method has been widely used globally [18,19]. It can be divided into the single and the dual crop coefficient method, the latter method divides the crop coefficient K_c into the basic crop coefficient K_cb and the soil evaporation coefficient K_e [15]. Previous research has shown that the dual crop coefficient method can more accurately gauge field ET_c due to its ability to estimate soil evaporation and vegetation transpiration [20].

Except for estimating ET accurately, considering the dramatic transformation in land use/cover and climate of China over the last 50 years, several recent works have also evolved to investigate the trend of ET patterns as well as ET sensitivity to meteorological factors. Li et al. [21] analyzed the space distribution of annual ET changes in China between 1980 and 2010 using four ET datasets. The results showed that ET increased massively in southeastern China while decreasing in northeastern China. Precipitation, net radiation, and vapor pressure difference were the primary climate variables influencing ET changes. Ma et al. [22] investigated water resource transformation in the Loess Plateau utilizing remote-sensing ET_a data and found that wind velocity and vegetation cover fraction (f_veg) had the greatest influence on ET_a. This was followed by barometric pressure, air dampness, precipitation, duration of sunshine, and temperature. However, there are still two inadequacies in such literature on China. First, extant analyses of the spatiotemporal patterns of ET and their climatic causes in China often rely heavily on datasets constructed before the 21st century. Since 2000, China has undergone booming innovation in agriculture, industry, and other fields, and the influence of agricultural practices on land-use changes is especially evident. The information extracted from last century’s data may be of limited value for current policymakers, resource managers, and other stakeholders. Second, previous studies have certain limitations in that they are mostly concentrated on specific regions of China, such as the Tibetan Plateau [23], Yanhe River Basin [6], and areas prone to desertification [22], primarily due to restricted availability of climate data and/or diverse research objectives. These studies do not fully address the needs of the whole country for water resource assessment and planning. Thus, one contribution of this paper is to investigate the spatiotemporal characteristics of ET_a changes and their climatic causes for all of China in the context of climatic changes since the 21st century. The spatial heterogeneity of the drivers of ET_a across the diverse regions of China is also discussed.

Previous research has mostly concentrated on the causes and influencing factors of changes of potential ET (ET_p) and/or ET_o, with relatively little attention paid to the spatiotemporal patterns and determining factors of ET_a. As stated by Luo et al. [1] and Nie et al. [24], ET_p is defined as the highest possible level of ET under specific meteorological conditions on a stable subsurface that has a plentiful water supply, while ET_a reflects ET levels under the real conditions of the surface. In other words, the water supply conditions of one district are crucial links between ET_p and ET_a. Unfortunately, despite its importance as an indicator of moisture status, soil water content is often difficult to measure accurately at national or even global scales. This has hindered researchers from accurately identifying spatiotemporal patterns of ET_a. Therefore, another contribution of this paper is the first attempt to use the Chinese Land ET_a Dataset (GPR-ET) generated through the Gaussian Process Regression method, published by Yin et al., to explore the changes and determining factors of ET_a [25]. The magnitude of ET_a relies on two main components: energy conditions and water supply conditions [26]. The energy conditions include factors such as temperature, humidity, and radiation. Our study selected easily obtainable temperature observations as proxies to reflect the impact of energy conditions on changes in ET_a. In addition, because of the difficulty in accurately measuring the surface moisture situation across China, we chose the most closely relevant variable, precipitation data, to describe the surface water supply conditions [27].

Due to the large span of China’s territory, there are significant climate and terrain differences among different regions. When exploring ET_a, in addition to the overall commonness of the whole country, the research on the characteristics of each region is equally important. Therefore, by using the geographic detector, the center of gravity migration model, and other methods, we adopted a new geographical division method. Using observed datasets from 496 meteorological sites, we aim to address the following questions: (1) What is the spatial and temporal distribution of ET_a across China during the initial decades of the 21st century and how will this distribution evolve in the future? (2) How can analysis methods such as correlation, sensitivity, and contribution analysis be utilized to quantify the influence of different meteorological factors (including precipitation and temperature) on variations in ET_a? (3) Are the dominant drivers of the change in ET_a the same in different regions? If not, how do they vary spatially? The aim of the current study is to offer quantitative guidance for the optimal allocation of water resources in the future. Section 2 presents the principles of various statistical analysis methods. The source and site descriptions of meteorological data and ET_a data are also summarized in this section. Section 3 presents the results, and further details are discussed in Section 4. Finally, our conclusions are presented in Section 5.

2. Materials and Methods

2.1. Materials

The meteorological data in the present study were collected from the National Earth System Science Data Center (http://www.geodata.cn/index.html, accessed on 29 April 2022), which included the precipitation and temperature (2000–2018) data. The spatial resolution of the above datasets was 1 km, and the geographic spatial range of the data sets was the nation of the mainland in China (excluding the islands and reefs in the South China Sea).

The GPR-ET dataset was used for our research. These data were generated by integrating the Gaussian process regression (GPR) algorithm and five process-based ET algorithms, with 36 maps per year and a resolution of 1 km. Yin et al. [25] demonstrated that the accuracy of this dataset is higher than the single calculation results of five process-based ET_a calculation methods and the existing eight high-resolution products. The five process-based ET_a calculation methods included Semi-empirical Penman–Monteith ET algorithm (SEMI-PM), Remote-sensing-based Penman–Monteith ET algorithm (RS-PM), Revised remote-sensing-based Penman–Monteith ET algorithm (RRS-PM), MODIS ET product algorithm (MOD16), and Penman–Monteith-Leuning version 2 (PMLv2). The eight products with high resolution are GLASS, SSEBop, BESS, MODIS, MTE, CR, and SEBS.

Table 1. Summary of information for the datasets used in this study.

Datasets	Time Range	Temporal Resolution	Spatial Resolution	Data Sources
Temperature (T)	1901–2022	monthly	1 km	National Earth System Science Data Center
Precipitation (P)
ETa	2000–2018	10 days	1 km	(https://doi.org/10.6084/m9.figshare.12278684.v5, accessed on 29 April 2022)

showed a summary of information for the dataset used in this study.

The existing seven geographical regions of China are divided into administrative units, which is convenient for the statistical of social and economic indicators. However, this zoning method leads to great differences in climate and geomorphological characteristics, which is not conducive to the analysis of regional heterogeneity. Therefore, in this research, we used the similarity of regional climate characteristics combined with geomorphological characteristics to divide China into seven new areas (Figure 1) according to the method devised by Wang et al. [28]. Accordingly, the zones comprised the northern region (N), the northeast region (NE), the southeast region (SE), the eastern area of the northwest (ENW), the western area of the northwest (WNW), the western area of the southwest (WSW), and the southwest region (SW).

2.2. Methods

In our study, we first studied ET_a spatiotemporally. In terms of time, the interannual ET_a trends were analyzed via simple linear regression [29]. Spatially, the Theil–Sen estimator and Mann–Kendall test were used to test the tendency for ET_a and its significance, respectively. Then, the center of gravity migration model was utilized to delve into the spatiotemporal changes of ET_a, and the coefficient of variation (CV) was also utilized to analyze relative volatility. Subsequently, to forecast the future tendency for ET_a, the Hurst index was calculated using rescaled range analysis (R/S). Finally, to explore the relevance between ET_a and meteorological elements, we adopted three correlative analysis methods (including simple, partial, and multiple correlative analysis) to explore the correlation between the two. In addition, the spatial heterogeneity was detected by the geographical detector, and the sensitivity coefficient and contribution rate were calculated.

2.2.1. Spatio-Temporal Characteristics of ET_a Changes

Theil–Sen Median Estimator and Mann–Kendall Test

The Theil–Sen median estimator and Mann–Kendall test are often adopted in meteorological and hydrological calculations to ascertain the tendency of long-term sequence data [30]. The data changes passed the 95% test (p < 0.05), indicating a significant change trend. Referring to Li et al. [31], and considering the range size of the ET_a trend, the trend can be divided into three categories: slope < −0.01 mm/year is a decreasing trend, slope > 0.01 mm/year is an increasing trend, and values between the two are considered unchanged.

Center of Gravity Migration Model

The gravity center is the location where, at any given time, a geographical factor achieves equilibrium in the space plane. It can reflect the movement and spatial aggregation of the factors and assess the development condition and trend of the factors [32]. Many scholars have used the gravity center migration model to calculate the changes in vegetation, meteorological hazards, and other factors [33]. In this study, grid point samples with a resolution of 20 m were generated within the study area. The following formula was used to determine the gravity center of ET_a over 19 years:

$$X_t=\frac{{\displaystyle \sum _{i=1}^n(E{T}_{ai}\times X_{ti})}}{{\displaystyle \sum _{i=1}^{n}E{T}_{ai}}}$$

(1)

$$Y_t=\frac{{\displaystyle \sum _{i=1}^n(E{T}_{ai}\times Y_{ti})}}{{\displaystyle \sum _{i=1}^{n}E{T}_{ai}}}$$

(2)

where n is the prime number of the ET_a value point and x_i, y_i, and ET_ai are the longitude, latitude, and ET_a value of point_i, respectively.

Coefficient of Variation

In addition to the variation trend values, the degree of variation in the data of long-term series should also be considered. The coefficient of variation (CV) is an indicator used to measure the degree of data dispersion, which represents the relative fluctuation of ET_a, and the calculation formula is [34]:

$$Cv=\frac{SD}{\overline{{ET}_a}}$$

(3)

where SD is the standard deviation of multi-year ET_a and $\overline{{ET}_a}$ is the multi-year mean value of ET_a.

2.2.2. Prediction of Spatiotemporal ET_a Trends

The Hurst index [35] obtained using the rescaled range (R/S) analysis method can effectively utilize the secular dependence in the time serial message to predict the future development trends of these data [35,36,37]:

For time series, $X=\left\{x_1,x_2,\cdots x_n\right\}$, the mean sequence is denoted as:

$${\overline{X}}_{(i)}=\frac{1}{i}\sum _{t=1}^i{x}_{(t)} i=1,2,\dots ,n$$

(4)

The cumulative deviation $X_{(t,i)}$ is denoted as:

$$X_{(t,i)}={\displaystyle \sum _{u=1}^i(x_{(u)}-x_{(i)})}\hspace{1em}1\le t\le i$$

(5)

The extreme difference $R_{(i)}$ is denoted as:

$$R_{(i)}=\underset{1\le t\le i}{\max } x_{(t,i)}-\underset{1\le t\le i}{\min } x_{(t,i)}$$

(6)

The standard deviation $S_{(i)}$ is denoted as:

$$S_{(i)}={\left[\frac{1}{i}{\displaystyle \sum _{t=1}^i{(x_{(t)}-x_{(i)})}^2}\right]}^{\frac{1}{2}}$$

(7)

For the ratio $R_{(i)}/S_{(i)}\cong R/S$, if $R/S\subset {\tau }^H$, it indicates a “Hurst” situation in the time series; the H value is the Hurst index. The specific value must be obtained through the least squares algorithm using the dual-logarithm coordinates system $\left(lni, ln R/S\right)$.

2.2.3. Attribution Analysis of ET_a

Correlation Analysis

Simple correlation analysis, partial correlation analysis, and complex correlation analysis can provide different levels of correlation information [38]. Simple correlation analysis can quantify the strength and orientation of the linear relationship between two variables, but it cannot control the effect of other variables. Partial correlation analysis is based on simple correlation analysis, which can eliminate the effect of other variables and provide a pure relationship between two variables [39]. Complex correlation analysis can comprehensively consider the common impact of multiple variables on the same goal. Therefore, using these three methods simultaneously can provide more accurate correlation information, thereby helping researchers to have a more comprehensive understanding of the relationships between variables [40].

• Simple correlative analysis: this method calculates the correlative coefficient between ET_a and climatic factors (precipitation and temperature [41]).

$$r_{xy}=\frac{\sum _{i=1}^n\left[\left(x_i-\overline{x}\right)\left(y_i-\overline{y}\right)\right]}{\sqrt{\sum _{i=1}^n{\left(x_i-\overline{x}\right)}^2}\sqrt{\sum _{i=1}^n{\left(y_i-\overline{y}\right)}^2}}$$

(8)

where $r_{xy}$ is the correlation coefficient of ET_a and temperature, $x_i$ and $y_i$ are the ET_a value and the corresponding value of temperature, respectively, and $\overline{x}$ and $\overline{y}$ are the averages of ET_a and temperature, respectively.

2.

Partial correlative analysis: this method can effectively obviate the interference of the third factor and describe the correlation between the first two more accurately [42,43]. The correlation coefficient and partial correlation coefficient <−0.2 represent negative correlations, those >0.2 represent positive correlations, and values between the two represent weak correlations.

$$r_{xy,z}=\frac{r_{xy}-r_{xz}{r}_{yz}}{\sqrt{\left(1-r_{xz}^2\right)\left(1-r_{yz}^2\right)}}$$

(9)

where $r_{xy,z}$ is the partial correlation coefficient between ET_a and temperature after fixed precipitation, and $r_{xy}$, $r_{xz}$, and $r_{yz}$ are the correlation coefficients between ET_a, temperature, and precipitation respectively.

The partial correlation coefficient’s significance can be evaluated using the T-test approach. The following is the formula:

$$t=\frac{r_{xy,z}}{\sqrt{1-r_{xy,z}^2}}\sqrt{n-m-1}$$

(10)

where m is the degree of freedom.

3.

Multiple correlation analysis: the multiple correlation coefficient can be used to explore the correlation of multiple factors [6,8].

$$r_{x,yz}=\sqrt{1-\left(1-r_{xy}^2\right)\left(1-r_{xz,y}^2\right)}$$

(11)

where the $r_{x,yz}$ is the multiple correlation coefficient.

The multiple correlation coefficient’s significance can be evaluated using the F-test approach:

$$F=\frac{r_{x,yz}^2}{1-r_{x,yz}^2}\times \frac{n-k-1}{k}$$

(12)

where the k is the number of independent variables.

Geographical Detector

A group of statistical techniques known as geographic detectors were used to identify spatial heterogeneity and the elements contributing to it. The fundamental premise of the theory is that the spatial pattern of the independent and dependent variables ought to be analogous if an independent variable (X) significantly influences a dependent variable (Y). To compensate for the limited number of samples, grid samples measuring 20 m were generated within the study area. The factors were discretized automatically using five categorized ways. The parameter for the geographical detector analysis was selected using the spatial scale with the highest q-value. The main categories of geographic detectors include ecological, risk, interaction, differentiation, and factor detection. Our research primarily utilized factor and interaction detection [44]. The q-value was determined using the following formulas:

$$q=1-\frac{{\displaystyle \sum _{h=1}^L{N}_h{\sigma }_h^2}}{N{\sigma }^2}=1-\frac{SSW}{SST}$$

(13)

$$SSW={\displaystyle \sum _{h=1}^L{N}_h}{\sigma }_h^2,SS{\rm T}={\rm N}{\sigma }^2$$

(14)

where h = 1, …, L is the classification of variable Y or factor X, $N_h$ and $N$ are the number of units of stratum h and the entire area, respectively, ${\sigma }_h^2$ and ${\sigma }^2$ are the variances of the variable Y for stratum h and the entire area, respectively, and $SSW$ and $SST$ are the sum of squares of the intralayer and the total sum of squares, respectively. The q value has a range of [0, 1]. The spatial heterogeneity of Y became more noticeable as the q value increased. The q value was greater when the independent variable X possessed more explanatory capacity on the characteristic Y if the stratification was caused by it.

Interaction detection: it evaluated whether the influence of variables X₁ and X₂ on Y were independent of one another or whether they functioned in concert to raise or reduce the dependent variable Y’s explanatory power. The relationships between the two factors can be categorized as follows (

Table 2. The type of interaction between two independent variables and the dependent variable.

Judge Theorem	Interaction
$q\left(X_1\cap X_2\right)<Min\left(q-X_1,q-X_2\right)$	Nonlinear weakening
$Min\left(q-X_1,q-X_2\right)<q\left(X_1\cap X_2\right)<Max\left(q-X_1,q-X_2\right)$	Single-factor nonlinear attenuation
$q\left(X_1\cap X_2\right)>Max\left(q-X_1,q-X_2\right)$	Double-factor enhancement
$q\left(X_1\cap X_2\right)=\left(q-X_1\right)+\left(q-X_2\right)$	Independence
$q\left(X_1\cap X_2\right)>\left(q-X_1\right)+\left(q-X_2\right)$	Nonlinear enhancement

):

Sensitivity Analysis

The partial differential model method was used to quantify the precise sensitivity connection between the meteorological elements and ET_a. The dimensionless sensitivity coefficient calculated by using this method can analyze the degree of impact of a single climate factor change on ET_a change [45]:

$$S{v}_i=\lim (\frac{\Delta E{T}_a/E{T}_a}{\Delta v_i/v})=\frac{\partial E{T}_a}{\partial v_i}\times \frac{v_i}{E{T}_a}$$

(15)

where ${Sv}_i$ is the sensitivity coefficient and the dimension is 1, and $v_i$ represents the independent variables of different climate factors.

Contribution Analysis

Contribution analysis can be used to evaluate the degree of impact of factors on overall changes and identify the main influencing factors. The specific contribution rate of corresponding meteorological elements to ET_a can be expressed by multiplying the sensitivity coefficient with the relative variation ratio of meteorological elements over the years [46]:

$$C{v}_i=S{v}_i\times R{v}_i$$

(16)

$$R{v}_i=\frac{n\cdot T{v}_i}{\left|\overline{v}\right|}\times 100\%$$

(17)

where $C{v}_i$ is the contribution rate (%) of the climatic factor $v_i$ and its change to ET_a, $R{v}_i$ is the multi-year relative variation ratio for the climatic factor $v_i$, $T{\nu }_i$ is the annual average linear change rate for the climatic factor, and $\overline{v}$ is the average of the climatic factor. The positive (negative) value of $C{v}_i$ indicates that the role of climate factors in ET_a change is positive (negative), and this role becomes stronger with the increase in the $C{v}_i$ absolute value.

3. Results

3.1. Temporal-Spatial Variations in ET_a

Figure 2 shows the linear trends for ET_a, temperature, and precipitation in China from 2000 to 2018. Annual ET_a ranged between 381.047 and 420.224 mm, with an average of 397.602 mm. From 2000 to 2018, the ET_a trend increase was 1.470 mm/year (p < 0.01). Temperature and precipitation also showed an increasing trend, with values of 0.014 °C/year and 2.583 mm/year (p < 0.10), respectively. Figure 3 shows the trend and significance of ET_a in China from 2000 to 2018. The rate of change in ET_a from 2000 to 2018 ranged from −78.471 mm/year to 153.291 mm/year, and negative values indicate a decreasing trend in ET_a in this region during the study period, with an increasing trend in 69.395% of the study region. The areas showing an increasing trend were concentrated in the southeast of WSW and the southern part of SW. The WNW region showed the largest proportion of area with an increase in ET_a at 88.208%. The regions with the largest proportions of area with a significant increase in ET_a were the N and ENW regions, with 37.843% and 33.669%, respectively, mostly found in the Loess Plateau area. Overall, the results showed that among the seven geographical divisions, the proportion of areas with a non-significant increase in ET_a is the highest, which further indicated that the main trend of ET_a in the country and various regions is increasing.

The coordinates of the center of gravity of ET_a were calculated for the study area and the seven geographic sub-regions by using the center of gravity migration model. Accordingly, the change in each center of gravity from 2000 to 2018 was plotted, and the results are shown in Figure 4. The center of gravity of ET_a generally migrated to the northwest throughout the study area. This was closely associated with the rise in precipitation in the Northwest since the 21st century. There were differences in the transformation of the center of gravity of ET_a within each climate subregion. The center of gravity of ET_a in WSW, N, SE, and NE showed a migration to the southwest in 2011. The same happened in ENW in 2012. The center of gravity of ET_a in the SW suddenly moved to the northeast in 2013. The reason was that the severe drought that occurred in SW in 2013 was concentrated in its southwestern region, which greatly reduced ET_a in SW and shifted the center of gravity to the northeast.

Figure 5 shows the spatial distribution of the coefficients of variation and stability for ET_a in China from 2000 to 2018. The coefficients of variation ranged between 0.001 and 4.243, indicating that ET_a variation possessed obvious spatial heterogeneity; the CV values varied significantly in different regions. The high CV region was predominantly centralized in the WSW region and its western and northern regions, as well as most of the western region of the Inner Mongolia Plateau. Among the five stability classifications, the region with low fluctuating stability was the largest, at 39.319%. The region with the lowest volatility stability accounted for only 5.553%, and most of them were located in the region surrounded by low volatility stability. The CV change in western China was relatively stronger than those in the east, and the WNW region, ENW region, and WSW region had the largest area share of the highest volatility changes, with area shares of 68.007%, 42.544%, and 30.865%, respectively. This phenomenon may be related to the significant increase in ET in Xinjiang and other regions due to global warming. In the other four regions, the areas with the highest proportion of area share were all low-fluctuation stable areas, and this fluctuation type had a share of over 60% in both the SE and NE regions.

Figure 6 shows the spatial distribution of the ET_a Hurst index (a) and future trends (b) in China, from 2000 to 2018. The average ET_a Hurst index was 0.472, and 60.831% of the region had a Hurst index of less than 0.5, suggesting that, in most regions, the future trends for ET_a were inconsistent with the current stage. The high-value area of the Hurst index (H > 0.60) accounted for 9.661%, and it was scattered throughout the country. The future change trends for ET_a in these areas were highly consistent with previous change trends. The median area of the Hurst index (0.40 < H ≤ 0.60) accounted for 66.362%, and the future trends for ET_a had the characteristics of weak co-direction and weak reverse.

Following Wu et al. [47], the future trends for ET_a can be obtained by overlaying the spatial change trends with the Hurst index classified as shown in

Table 3. The future trends of ET_a.

Number	β	Hurst	Future Trend of ETa	Meaning
I	<−0.01	<0.5	Anti-persistent increase	The declining ETa in the past will increase in the future
II	−0.01~0.01	<0.5	The trend of change cannot be determined	/
III	>0.01	<0.5	Anti-persistent decrease	The increasing ETa in the past will decline in the future
IV	<−0.01	>0.5	Persistent decrease	The declining ETa in the past will persistently decrease in the future
V	0.01~0.01	>0.5	Basically unchanged	The degree of change in ETa is relatively small
VI	>0.01	>0.5	Persistent increase	The increasing ETa in the past will persistently increase in the future

. Figure 3 and Figure 6a were applied to the spatial superposition analysis to form the distribution of future tendency for ET_a (Figure 6b). Future trends of the ET_a for 0.671% of the study area were undetermined. From 2000 to 2018, the ET_a in these regions fluctuated greatly, and the complex pattern of change made the trends for ET_a in these regions difficult to predict. The percentages of areas showing an anti-persistent increase and anti-persistent decrease in the WSW region were close (33.323% and 32.402%, respectively). The type with the highest proportion in the SW area had a persistent increase (44.174%), while in other regions, this type had an anti-persistent decrease (Figure 7). The results above demonstrated that, in the future, there will be more areas with decreasing trends (52.981%) than those with increasing trends. Areas that showed a continuous decreasing trend were predominantly centralized in Taiwan Province, the central-eastern area of WSW, the mountain range of northern NE, and the hilly area of SE.

3.2. Correlation between ET_a and Meteorological Factors

Figure 8 shows the significance of the correlation and partial correlation between ET_a and temperature (a) and precipitation (b). The percentages of the places where ET_a and temperature were positively and negatively correlated were 34.917% and 16.309%, respectively, while ET_a and precipitation were positively and negatively correlated for 48.325% and 23.884% of places, respectively. The results showed that ET_a correlated positively with both meteorological factors in a spatially wider area, and the precipitation correlated more positively than temperature with ET_a. The area with a positive correlation between ET_a and temperature was 10.348%, primarily concentrated in the central region. This correlation type occupied the largest area (36.422%) in the SE region. In the N region and NE region, the large continuous distribution of ET_a showed a highly significant positive correlation with precipitation, and in the ENW region, the area of significant positive correlation even exceeded half of the entire space, which was related to the low precipitation in these regions. In addition, the areas where ET_a demonstrated the same correlation with temperature and precipitation were largely situated in the southwestern area of ENW, the southeastern part of SW, and other areas. This phenomenon further suggested that precipitation and temperature were the primary meteorological elements that influenced regional ET_a and led to the divergence of spatial and temporal patterns, while temperature and precipitation played a common role in ET_a in these regions.

The significance of the partial correlation between ET_a and temperature and precipitation is shown in Figure 8c,d. After eliminating their reciprocal influences, there was no significant negative correlation between temperature or precipitation and ET_a. The percentages of areas that correlated positively and negatively with ET_a and temperature were 37.585% and 15.784%, respectively, while 49.150% and 23.730% of areas correlated similarly with ET_a and precipitation, respectively. The pattern of partial correlation significance for ET_a and temperature and precipitation was similar to that seen in the distribution of correlation significance map for the positive correlation area, the difference mainly appeared in the significant negative correlation area. In the NE, SW, and SE regions in particular, the areas with the non-significant negative correlation type increased significantly once the influence of temperature was eliminated. In general, compared to ET_a and temperature, there was a larger partial correlation between ET_a and precipitation.

Figure 9 shows the multiple correlation coefficient and its significance for ET_a and temperature and precipitation. The multiple correlation coefficient values for each region cover a large span, and the spatial average was 0.470. Only 9.291% of the regions had a weak multiple correlation, indicating that in most regions, temperature and precipitation may have a complex combined effect on ET_a. From the p value distribution map, it can be seen that the areas passing the significance test (p< 0.05) were primarily in the WNW, ENW, N, and SE regions, and were centered in the Inner Mongolia Plateau area, accounting for 34.867% of the entire space.

Based on the research of previous scholars on the driving factors of NPP [48], NDVI [49], and ET_a [50],

Table 4. Driving factor division criteria.

Driving Factors	α = 0.05
	Partial Correlation Coefficient (ETa-P)	Partial Correlation Coefficient (ETa-T)	Multiple Correlation Coefficient
Other-factors-driven			F < F0.05
Precipitation-driven	|t| > t0.05		F > F0.05
Temperature-driven		|t| > t0.05	F > F0.05
Joint precipitation- and temperature-driven	|t| < t0.05	|t| < t0.05	F > F0.05

shows the zoning rules of climate factors on ET_a.

The driving force partition was determined (Figure 10) following the division standard of Table 4. The figure shows that the driving force partition showed strong spatial heterogeneity. Precipitation-driven regions (II) were predominantly centralized in the WNW, ENW, and N areas. ET_a in these regions was mainly limited by water supply, so ET_a was driven by precipitation, having a total area of 18.701%. The area of the temperature-driven regions (III) accounted for only 8.050%, indicating that the temperature-driven zones are not significant compared to other types in the country. The temperature-driven regions were mostly found in the SE area. Water storage and precipitation in this region were sufficient, and temperature mainly affected evaporation capacity. The joint precipitation- and temperature-driven regions (IV) were discretely distributed throughout China, accounting for 8.131% of the total area. Other-factors-driven regions (I) accounted for the largest proportion. In these areas, changes in ET_a were not primarily driven by temperature or precipitation. More factors are required to further explore the mechanism.

3.3. Time Variation of Factors Detected by Geographic Detectors

Figure 11 analyzed the effect of temperature and precipitation separately and jointly on evapotranspiration using q -value, and the larger the q-value, the stronger the explanatory power of the corresponding element of the changes in ET_a. The figure shows that the q-values range from 0.036–0.586 for temperature, 0.027–0.815 for precipitation, and 0.093–0.847 for temperature and precipitation interaction.

Table 5. Annual average q value.

Geographical Division	q-T	q-P	q-T&P
The whole study area	0.407	0.739	0.779
WSW	0.214	0.461	0.499
WNW	0.191	0.374	0.458
ENW	0.166	0.780	0.818
SW	0.272	0.215	0.402
N	0.468	0.617	0.683
SE	0.258	0.196	0.370
NE	0.128	0.371	0.500

shows the multi-year average values of the corresponding q-values, and the results indicate that the interaction between two factors on ET_a changes is stronger than that of a single factor in all regions. The results of the analyses for each year in the remaining regions indicated that precipitation explains the spatial heterogeneity of ET_a more strongly than temperature, except for the SW and SE regions. The outcomes of combining temperature and precipitation detection on ET_a in the SW region (2007, 2013, and 2017) and the SE region (2007) were nonlinearly augmented, while the results of the rest of the regions and years were double factor enhancement. This result further suggested that there is a nonlinear enhancement of the influence of temperature and precipitation on ET_a, in addition to spatial differences.

3.4. Sensitivity Analysis of ET_a and Meteorological Factors

The spatial distribution of the ET_a sensitivity type to temperature is shown in Figure 12a. In general, the sensitivity coefficient of ET_a to temperature was greater in the southeast than northwest. With regard to the sensitivity level, the area with the largest proportion had high co-directional sensitivity, at 35.135%. This sensitivity type accounted for 70.248% of the area in the SE region and was also widely distributed in the SW and N regions. ET_a in these areas was highly sensitive to temperature, and both ET_a and temperature changed in the same direction. In the NE region, the ET_a was also sensitive to temperature, but the direction of change of the two was the opposite.

The spatial distribution of ET_a sensitivity types to precipitation is shown in Figure 12b, and its average value was 0.15. The sensitivity coefficient of ET_a to precipitation showed a continuous high-value area in the northern and western regions of China (including WSW, WNW, ENW, N, NE, and SW). High-value concentration areas also appeared in eastern Sichuan and northern Chongqing in the SW region. Unlike temperature, the coefficient of ET_a sensitivity to precipitation was greater in the northwest than in the southeast. The high co-directional sensitive area accounted for the largest proportion of all types of sensitive levels, at about 48.845%, followed by the high reverse sensitive area, which accounted for more than one-fifth of the area. The regions with medium co-directional sensitivity and negligible sensitivity were distributed around the regions with high co-directional sensitivity and accounted for a small proportion of the regions.

3.5. Contribution Analysis of Climate Factors to ET_a

The rates at which climatic conditions contributed to ET_a are displayed in Figure 13a,b. Overall, precipitation and temperature contributed to 3.015% and 0.240% of ET_a, respectively. Those results indicated that the overall upward trends for ET_a were caused by both temperature and precipitation, and the upward trends for precipitation were greater than those for temperature. The temperature contribution rate to ET_a was between −3% and −3% in 87.728% of the area. The high contribution rate was concentrated in the area near Hulun Lake in Inner Mongolia in the WNW, N, and ENW regions. The means and trends for temperature in the southern and eastern parts of WSW showed complex regularity due to their particular environmental conditions, resulting in a large difference in the temperature contribution rate to ET_a in these areas and surrounding areas. For 87.762% of the region, the precipitation contribution rate to ET_a was concentrated between −6% and 16%. The high-value regions were primarily concentrated on the Inner Mongolia Plateau and were also distributed in central and southern Xinjiang and western WSW. The contribution rate for the northern part of the NE, as well as the south of the Yangtze River’s middle and lower reaches in the SE region, was low. In the original area of rivers in the WSW, large fluctuations in elevation led to significant temperature differences, causing their contributions to ET to vary in a complex manner. In the east area of the 110°E meridian and below 500 m above sea level, precipitation contributed negatively to ET, while temperature mainly contributed positively to ET, and the absolute contribution rate of precipitation was even greater.

4. Discussion

4.1. The Spatiotemporal Characteristics of Evapotranspiration

In this research, the mean and temporospatial changes in ET_a from 2000 to 2018 in China were analyzed. We found that the average ET_a was 397.602 mm/year, showing an increasing trend with a rate of 1.470 mm/year. However, due to inconsistent data sources and the selection of start and end years, these results were slightly different from other previous studies. Cheng et al. [51] utilized MODIS products to calculate the average ET_a value in China from 2001 to 2018. They calculated the value as 359.61 mm/year, with a linear variable rate of 2.95 mm/year, and noted that MODIS products underestimated the ET_a value in China. In addition, unlike the downward trend for both ET_p [52] and ET_a [53] in the mid to late 20th century, our study demonstrated that the average annual ET_a showed a fluctuating and increasing trend after entering the 21st century [54,55]. The main reason for this discrepancy is that, in the early 21st century, China made significant resource investments in ecological restoration, which included the implementation of returning farmland to forests and grasslands, returning grazing to grasslands, natural forest protection, and protective forest system engineering [56]. These projects have achieved good results with improved vegetation conditions. Accordingly, the areas with significant increases in ET_a mainly occur in areas with increased vegetation coverage [57].

At the regional scale, the ET_a tendency in China was different. The increase in ET_a mainly occurred in the North China Region [58], the Loess Plateau [59], and the Huai River Basin [60], while the significant decrease in ET_a mainly occurred in the Tibetan Plateau [23], the Hai River Basin [61], and the northern Xinjiang [57]. Our research results also showed that in most of the Qinghai Tibetan Plateau, ET_a displayed a declining pattern, while there was a marked growth in the Loess Plateau region. In our research, the annual average ET_a in China decreased from the southeast seaboard to the northwest outback, which coincides with the conclusions of He and Shao [57], Li et al. [62], and Cheng et al. [51]. Moreover, the north and west had higher CV values in ET_a than the south and east. The high fluctuation changes in ET_a were concentrated in the northwest, while ET_a fluctuation was the most stable in most of the eastern regions. This was consistent with previous research results [51,63].

Our research indicates that the correlation between precipitation and ET_a is stronger than that of temperature, both with and without the influence of the other variable. In a study conducted in the upper and intermediate reaches of the Heihe River Basin, Qiu et al. [64] also discovered a high partial correlation between precipitation and ET. In addition, in our research, ET_a in the western and northern areas of China was primarily driven by precipitation, while in the southeastern regions, it was mainly driven by temperature. The combined driving force of precipitation and temperature was distributed in discrete patches across various regions of the country. ET_a was more sensitive to precipitation than to temperature, and precipitation affected it in more regions. In the area near Hulun Lake in Inner Mongolia and the border area between western Heilongjiang and Inner Mongolia, the temperature contribution rate to ET_a was higher. The high precipitation contribution rate to ET_a was mainly concentrated in the Inner Mongolian Plateau. Fu et al. [54] also found that temperature and precipitation are the considerable factors causing variations in ET, with precipitation being the largest driving factor for ET changes. Furthermore, precipitation and ET_a have a specific relationship that is influenced by factors such as the intensity and duration of precipitation [65], soil type [66], vegetation type [67], etc., resulting in a lag effect. However, this lag effect often persists in the short term [68]. Unfortunately, the precipitation and ET_a data used in this study are both monthly averages. Due to the limitations of available data, the lag effect of precipitation and ET_a may not be significant on a monthly scale. The impact of the lag effect between precipitation and ET_a on their correlation will be explored in our future research.

By studying the changes in ET_a, we can better understand the water use of vegetation and soil, thereby guiding the rational allocation and management of agricultural water resources. Based on the monitoring results of ET_a, the irrigation water amount can be adjusted promptly to avoid excessive irrigation or water resource waste, thus, improving the efficiency of water resource utilization and promoting sustainable development of agricultural production. We found that more than half of the study areas have inconsistent trends in future ET_a changes compared to the current situation, and the water resource allocation in these areas needs to be adjusted promptly. In the anti-persistent decrease/increase area, ET_a will change the current pattern and show a decrease/increase phenomenon. In this area, it is necessary to moderately control the irrigation water amount to avoid waste/crop water shortage.

Because we only considered two meteorological factors in this study, there were some normal discrepancies between our results and those of previous studies. Most previous studies showed that the variation in ET_a was significantly influenced by precipitation [54,69], but temperature-driven [70,71] and jointly driven [62] factors were also found. For example, Su et al. [53] found that the decrease in ET_a in Northeast China was led by changes in precipitation, while the decrease in ET_a in the eastern Qinghai Tibet Plateau was caused by changes in surface conditions.

4.2. Limitations and Future Prospects

Our study only analyzed two driving factors for ET_a, namely temperature and precipitation. Although previous studies have confirmed that these two factors exert a significant impact on ET_a [72,73], the changes in ET are also influenced by other factors. These include other meteorological factors, terrain factors, vegetation coverage, and human activities. In our study, the area driven by other factors was the largest, which also confirms this observation.

The meteorological factors that dominate ET changes vary in different regions and at different scales. Apart from temperature and precipitation, wind velocity and solar radiation have also been shown to be considerable driving forces for ET changes [54,62]. Nie et al. [24] used meteorological data from Heilongjiang Province from 1960 to 2019 to analyze the climate factors affecting ET_o using the sensitivity contribution rate approach. Their results indicated that the highest temperature, average temperature, wind velocity, and solar radiation positively drove ET_o, with wind speed being the primary element. Due to the lack of research on wind speed in our study, other factors had the largest effects on the area in the Heilongjiang region.

Altitude is one of the terrain factors that has a remarkable impact on ET. Liu et al. [74] analyzed the differences in ET_o changes and climate-driven mode at different altitudes and found that the annual ET_o sequence of the Qinghai Tibet Plateau presented an upward tendency with this elevation increase. The increase is more significant at an elevation below 5000 m. Ma et al. [75] explored soil moisture, solar radiation, and ET_o data at elevations from 3300 to 4700 m in the Qinghai Lake Basin. They found that, as the altitude increases, the main control factor for ET_o shifts from water conditions to energy conditions. According to Figure 1 and Figure 13, the positive contribution of precipitation to ET_a dominates at altitudes between 1000 and 2500 m. As the altitude continues to increase, the impact of temperature on ET_a gradually becomes stronger than precipitation. The results above show that the influence of altitude on ET_a can be indirectly achieved by affecting the climate.

Different substrates are suitable for different types of vegetation growth; moreover, the presence of vegetation can alter the surface’s albedo and roughness, as well as intercepting and transpiring water. This makes vegetation an important factor affecting ET [76,77]. In our study, ET_a in the Loess Plateau region was driven by other factors, except for some areas affected by precipitation. Similar to previous studies, it was shown that these areas were greatly affected by vegetation cover. Jin et al. [77] isolated the influence of climatic variation and forest restoration on ET by adjusting variables. They discovered that the primary factor causing ET changes on the Loess Plateau was vegetation restoration.

The Hurst index analysis used in our study can quantitatively describe the long-range dependence of time series, and this method can determine the future trend of time series change [78]. However, Hurst index analysis is based on historical data, revealing long-term correlations while ignoring future situational changes. When estimating the future trends of ET_a, these scenario changes include changes in land use policies, the occurrence of extreme weather, etc. Our research is based on further analysis of existing remote sensing ET_a datasets, and many scholars have made contributions in directly using remote sensing technology to estimate ET. Yang et al. [79] used thermal imaging satellite remote sensing and Wang et al. [80] used satellite microwave remote sensing to effectively estimate ET. At the same time, with the development of machine learning algorithms and scene simulation technology, these two techniques have also been applied by many scholars to estimate ET [81,82].

In the past few years, human activities, including agricultural irrigation, soil and water conservation projects, and so on, have exerted an increasing influence on ET. However, human activity’s effects on ET are difficult to measure accurately due to the lack of systematic and complete indicators. Given that interactions can occur between climate variation, vegetation, and human endeavors, future studies should consider the interactions between various factors. Further research on these factors and the contribution of ET_a are required to reduce the uncertainty of such relationships. The results of this paper can be used as a theoretical foundation for the management and distribution of water resources in China and as a reference for a follow-up study of ET_a climate feedback mechanisms. It should be noted, however, that our study only used two meteorological factors, temperature and precipitation; therefore, the ET_a change mechanism could not be explained in more detail. These deficiencies should be addressed further to improve future work in this area.

5. Conclusions

Based on GPR-ET data and temperature and precipitation interpolation data, our research analyzed the spatiotemporal features of ET_a in China from 2000 to 2018. On this basis, the connection between ET_a and climate change was studied, and the climatic variation contribution rate to ET_a variation was quantitatively analyzed. The conclusions are as follows:

• From 2000 to 2018, the multi-year spatial mean for ET_a was 397.602 mm, which showed a significant upward trend. The variation characteristics had significant spatial differences, while the southeast and northwest have high and low general distributions, respectively. The center of gravity of ET_a generally shifts towards the northwest direction. The spatial distribution of the ET_a variation coefficient decreased from northwest to southeast. The enforcement of public ecological restoration measures and the acceleration of urbanization may be the main reasons for ET_a fluctuations.
• In the future, the tendency of ET_a changes in most regions will be inconsistent with the current stage. A total of 52.981% of the regional ET_a will decrease to varying degrees, and the future tendency of ET_a in the total area of 0.671% is still uncertain. Combined with the results of the coefficient of variation, it can be seen that the ET_a in these areas fluctuated greatly from 2000 to 2018, and the complex changes make the ET_a trend in this part of the region difficult to predict.
• Spatially, the correlation between ET_a and precipitation was stronger than that of temperature, and precipitation had a positive correlation with ET_a in more regions. In the division of ET_a driving factors, the precipitation-driven region was predominantly spread in the western and northern regions, and the regional proportion was larger than the temperature-driven region. It is worth noting that in the national and seven regional studies, the driving force of the interaction between temperature and precipitation is stronger than a single factor.
• ET_a was more sensitive to temperature in the southeast, while it was more sensitive to precipitation in the northwest. The coefficient for ET_a sensitivity to precipitation was higher in more regions. The positive contribution area of temperature and precipitation was greater for ET_a than the negative contribution area. In the southeast of China, temperature was the major cause of the ET_a increase, whereas in the north and northwest, the role of precipitation was more significant.

To date, researchers have investigated the impact of multiple factors on ET_a, but few studies systematically analyzed the relationships between various factors. Future research should consider the overall study of complex relationships and the isolation of single factors. In addition, although existing models have developed rapidly, there is insufficient research to demonstrate the reliability of any single method in ET_a prediction; therefore, additional research is still required in this field.