Quantifying the Spatiotemporal Variation of Evapotranspiration of Different Land Cover Types and the Contribution of Its Associated Factors in the Xiliao River Plain

: Evapotranspiration (ET) is a vital constituent of the hydrologic cycle. Researching changes in ET is necessary for understanding variability in the hydrologic cycle. Although some studies have clariﬁed the changes and inﬂuencing factors of ET on a regional or global scale, these variables are still unclear for different land cover types due to the range of possible water evaporation mechanisms and conditions. In this study, we ﬁrst investigated spatiotemporal trends of ET in different land cover types in the Xiliao River Plain from 2000 to 2019. The correlation between meteorological, NDVI, groundwater depth, and topographic factors and ET was compared through spatial superposition analysis. We then applied the ridge regression model to calculate the contribution rate of each inﬂuencing factor to ET for different land cover types. The results revealed that ET in the Xiliao River Plain has shown a continuously increasing trend, most signiﬁcantly in cropland (CRO). The correlation between ET and inﬂuencing factors differed considerably for different land cover types, even showing an opposite result between regions with and without vegetation. Only precipitation (PRCP) and NDVI had a positive impact on ET in all land cover types. In addition, we found that vegetation can deepen the limited depth of land absorbing groundwater, and the inﬂuence of topographic conditions may be mainly reﬂected in the water condition difference caused by surface runoff. The ridge regression model eliminates multicollinearity among inﬂuencing factors; R 2 in all land cover types was over 0.6, indicating that it could be used to effectively quantify the contribution of various inﬂuencing factors to ET. According to the results of our model calculations, NDVI had the greatest impact on ET in grass (GRA), cropland (CRO), paddy (PAD), forest (FOR), and swamp (SWA), while PRCP was the main inﬂuencing factor in bare land (BAR) and sand (SAN). These ﬁndings imply that we should apply targeted measures for water resources management in different land cover types. This study emphasizes the importance of comprehensively considering differences among various hydrologic cycles according to land cover type in order to assess the contributions of inﬂuencing factors to ET. m to 1200 m. The annual average temperature is 5.5 ◦ C to 7.50 ◦ C, sunshine hours are 2800–3100 h, and the RH is 45–58%, with a moderately temperate, semi-arid monsoon climate. The spatial and temporal distribution of PRCP is uneven; it spatially increases from northwest to southeast and is heavily concentrated in the period from June to September. The average annual precipitation is between 325 mm and 729 mm [37,52].


Introduction
Evapotranspiration (ET) refers to the process of water entering the atmosphere in gaseous form from the underlying surface, including from evaporation from the soil and the transpiration of water through plant stomata [1][2][3]. ET is the main pathway for energy and water exchange in the soil-vegetation-atmosphere system [4,5]. About 60% of atmospheric precipitation on the Earth's surface is returned to the atmosphere through ET; in arid The Xiliao River Plain (XRP) is an important development base for agriculture and animal husbandry in China. Due to the dual impact of climate change and human activities in recent years, the runoff of the Xiliao River has decreased year by year. Especially since 2000, the downstream section of the Xiliao River has dried up due to water diversions. The shortage of surface water resources has aggravated land desertification and deteriorated the ecological environment in the XRP [51]. Consequently, water loss through ET in the XRP has attracted the attention of many scholars. However, as the factors affecting the changes in ET are complicated, qualitative analyses of the contribution of climate change to ET remain limited. In this study, in addition to meteorological factors, we selected more comprehensive influencing factors such as groundwater depth and vegetation coverage factors to quantitatively evaluate the contribution of each influencing factor to ET for different land cover types. Specifically, this study is designed to (1) detect differences in temporal and spatial changes of ET in different land cover types from 2000 to 2019, (2) compare the sensitivity of ET to influencing factors in different land cover types, and (3) quantify the contribution of each influencing factor to ET and determine the main driving factors in different cover types. The results will not only help us gain a better understanding of water cycle processes in different land cover types, but also help policymakers to devise targeted measures for water resources management for different land cover types.

Study Area
The study area (42 • 18 -44 • 30 N, 119 • 14 -123 • 42 E) was located.0. in the south-west of the Northeast Plain of China, with a total area of about 96,000 km 2 (Figure 1), consisting of the alluvial plain and residual sandy ancient alluvial plain formed by the Xiliao River and its tributaries. Hilly grassland and loess hills are the main geomorphic types in the study area. The terrain of the study area slowly inclines from west to east and presents wavy undulation, with elevation ranging from 38 m to 1200 m. The annual average temperature is 5.5 • C to 7.50 • C, sunshine hours are 2800-3100 h, and the RH is 45-58%, with a moderately temperate, semi-arid monsoon climate. The spatial and temporal distribution of PRCP is uneven; it spatially increases from northwest to southeast and is heavily concentrated in the period from June to September. The average annual precipitation is between 325 mm and 729 mm [37,52].
land cover types in order to clarify the water and energy interactions between land and atmosphere covered by various vegetation types.
The Xiliao River Plain (XRP) is an important development base for agriculture and animal husbandry in China. Due to the dual impact of climate change and human activities in recent years, the runoff of the Xiliao River has decreased year by year. Especially since 2000, the downstream section of the Xiliao River has dried up due to water diversions. The shortage of surface water resources has aggravated land desertification and deteriorated the ecological environment in the XRP [51]. Consequently, water loss through ET in the XRP has attracted the attention of many scholars. However, as the factors affecting the changes in ET are complicated, qualitative analyses of the contribution of climate change to ET remain limited. In this study, in addition to meteorological factors, we selected more comprehensive influencing factors such as groundwater depth and vegetation coverage factors to quantitatively evaluate the contribution of each influencing factor to ET for different land cover types. Specifically, this study is designed to (1) detect differences in temporal and spatial changes of ET in different land cover types from 2000 to 2019, (2) compare the sensitivity of ET to influencing factors in different land cover types, and (3) quantify the contribution of each influencing factor to ET and determine the main driving factors in different cover types. The results will not only help us gain a better understanding of water cycle processes in different land cover types, but also help policymakers to devise targeted measures for water resources management for different land cover types.

Study Area
The study area (42°18′-44°30′N, 119°14′-123°42′E) was located.0. in the south-west of the Northeast Plain of China, with a total area of about 96,000 km 2 (Figure 1), consisting of the alluvial plain and residual sandy ancient alluvial plain formed by the Xiliao River and its tributaries. Hilly grassland and loess hills are the main geomorphic types in the study area. The terrain of the study area slowly inclines from west to east and presents wavy undulation, with elevation ranging from 38 m to 1200 m. The annual average temperature is 5.5 °C to 7.50 °C, sunshine hours are 2800-3100 h, and the RH is 45-58%, with a moderately temperate, semi-arid monsoon climate. The spatial and temporal distribution of PRCP is uneven; it spatially increases from northwest to southeast and is heavily concentrated in the period from June to September. The average annual precipitation is between 325 mm and 729 mm [37,52].   As the ET data source of the study, the MOD16 global terrestrial ET product data set acquired from the United States Geological Survey (https://earthexplorer.usgs.gov/ (accessed on 5 September 2021)) includes ET, potential evapotranspiration (PET), latent heat flux (LE), and potential latent heat flux (PLE). The data set has a spatial resolution of 500 m, and temporal resolutions are 8-day and annual syntheses [53,54]. MOD16 is based on the P-M equation to calculate wet canopy evaporation (E wet , mm·day −1 ), bare soil evaporation (E s , mm·day −1 ), and vegetation transpiration (E T , mm·day −1 ) as ET (E, mm·day −1 ) by using remote sensing data (reflectance, vegetation coverage rate, and land cover type) and meteorological data (RH, TEMP, and air pressure) from measuring stations [25,55]. The calculation formulas are as follows: where ρ (kg·m −3 ) is the air density, C p (MJ·kg −1 · • C −1 ) is the specific heat capacity of air, f c is fractional vegetation cover, VPD (kPA) is the atmospheric vapor-pressure deficit, f wet is the relative surface wetness, Pa (kPA) is the atmospheric pressure, γ cv (s·m −1 ), γ hrc (s·m −1 ) represents the surface resistance and aerodynamic resistance to evaporated water on the wet canopy surface, ε is the ratio of the molecular weight of water to dry air (i.e., 0.622), γ tot (s·m −1 ) is the total aerodynamic resistance, γ as (s·m −1 ) is the aerodynamic resistance at the soil surface, f sm is the soil water constrain, γ a (s·m −1 ), γ s (s·m −1 ) represents the canopy surface resistance and the aerodynamic resistance between the mean canopy height and the air above the canopy, and ∆ is the derivative of saturated vapor pressure to temperature [55,56].

Land Cover Data
Landsat remote sensing data were used to extract information regarding the spatial distribution of the various land cover types in XRP. The remote sensing data were derived from the United States Geological Survey, which contained Landsat7-ETM data from 2000 to 2013 and Landsat8-OLI data from 2014 to 2019. Images with no noise or strips and less than 5% cloud cover were selected from June to September every year. This is the period of vigorous plant growth; as such, the color and texture of the images are richer, and the characteristics of different land cover types are more obvious. Land cover types could be accurately identified according to their color, texture and other features. The selected remote sensing images were processed by ortho-correction, layer stacking, image sharpening, etc. According to the classification system of land cover data in China and the actual distribution of land cover types in the XRP, the land cover types in the study area were divided into nine categories: bare land (BAR), sand (SAN), grass (GRA), swamp (SWA), forest (FOR), cropland (CRO), paddy (PAD), water area, and construction land. According to the color, shape, texture, etc., of the remote sensing images, interpretation marks for different land cover types were established. With the support of the GIS platform, the human-machine interactive interpretation method was used to obtain the spatial distribution data of land cover types in the study area ( Figure 1b).

Influencing Factors for Analysis Data
Wind direction speed (WDSP), PRCP, TEMP, RH, NDVI, groundwater depth, elevation, and slope were selected as factors influencing ET. Meteorological data were obtained from the China National Meteorological Information Center (http://www.nmic.cn/ (accessed on 10 September 2021)), which comprised the WDSP (m·s −1 ), TEMP ( • C), RH (%), and PRCP (mm) of 22 meteorological stations in and around XRP from 2000 to 2019. NDVI data were derived from the MOD13 products of the United States Geological Survey. The original data were first subjected to bidirectional reflectance atmospheric correction to remove the influence of water, cloud, aerosols, and cloud shadows. The spatial resolution was 500 m, and the temporal resolution was both monthly and annual. From the Geological Cloud Information Service Platform of the China Geological Survey, data from 2019 from 138 groundwater depth monitoring stations distributed throughout the XRP and surrounding areas were used to determine groundwater depths. Elevation and slope data were derived from NASA Earthdata (https://search.earthdata.nasa.gov/ (accessed on 10 September 2021)), which provided 30 m spatial resolution DEM data covering global land.
Based on the GIS platform, the MOD16 data set was processed by ET data extraction and monthly scale data synthesis. After preprocessing by projection transformation, spatial interpolation, splicing, and clipping, the data of each influencing factor were uniformly resampled to 500 m based on the spatial resolution of the ET data. Since the MOD16 ET product used the leaf area index to indirectly reflect soil moisture content, there were no data for water and construction areas. Combined with the land cover data in the study area, we selected BAR, SAN, GRA, SWA, FOR, CRO, and PAD to analyze the temporal and spatial variation characteristics and influencing factors of ET.

Trend Analysis
In previous studies, scholars mostly adopted Sen's slope to reflect the linear trends of ET changes [57]. However, without a large number of sample points, the method is easily affected by outliers. This study used Theil-Sen Median to analyze the change trends of ET in each pixel. By taking the median as the function, the Theil-Sen Median trend can reduce the impact of noise due to its insensitivity to errors and outliers [58]. The equation is as follows: where j and i are ordinal numbers of time, ET j and ET i represent the ET of each pixel in year j and year i, β represents the trend of change, β > 0 indicates an upward trend, and β < 0 a downward trend.

Significance Test
The Mann-Kendall method was used to conduct a significance analysis of the Theil-Sen trend in each pixel. This method does not require sequence data to obey normal distribution or linear trends, and it is not significantly affected by a few missing values or outliers. This method has been widely used in trend significance tests of long-term sequence data. The statistical test method is as follows: for time series ET = (ET 1 , ET 2 , . . . , ET n ), the normalized test flux Z is defined as: Remote Sens. 2022, 14, 252 Var(s) = [n(n − 1)(2n + 5)]/18 (9) where n is the research period (20 years from 2000 to 2019), ET j and ET i are the values for each pixel sequence, α is defined as the significance test level, and Z (1−α)/2 is defined as a standard normal variance. When |Z| ≥ Z (1−α)/2 , the null hypothesis is rejected. We set significance levels at α = 0.01, 0.05, and 0.1 [59,60].

Correlation Analysis
Pearson correlation coefficients (PCCs) were used to calculate the correlation coefficient between ET and various meteorological factors and NDVI in each pixel [61]. The formula is as follows: where n is the research period (i.e., from 2000 to 2019), i is the number of years, x i and y i are the values of two factors x, y in the i-year, and x, y are the n-years average of two factors.

Ridge Regression
Ridge regression is a biased estimation regression method for multiple collinearity data analyses. At the cost of losing part of the information or the fitting accuracy, ridge regression abandons the lack of bias of the least-square method to obtain more stable estimated results and a smaller mean square error. The ridge regression equation is as follows:β (k) = X X + kI −1 X y (11) where X X is the coefficient matrix. When the independent variable has multiple collinearities, i.e., X X ≈ 0 mathematically, a constant matrix kI (k > 0) is added to significantly reduce the possibility of X X + kI being singular. Therefore, the mean square error is smaller. k is the ridge parameter whose value is not unique; theoretically, an optimal value exists. The judgment criterion was to select the minimum k as the ridge parameter when the ridge estimation of each regression coefficient was stable by observing the ridge trace graph [62]. The equation for ridge regression is as follows: If we standardize y i and x i , then z y = y i −ŷ σ y , z xj =

The normalized equation is as follows:
z y =b 1 z x1 +b 2 z x2 + · · · +b k z xk (13) The relative contribution rate of each influencing factor to ET is calculated by Equation (14).
where C h is the relative contribution rate of a certain influencing factor to ET, andb k is the regression coefficient of each influencing factor in the standardized equation. The absolute value ofb k represents the degree of influence of each factor on ET, and the symbol represents the positive or negative effects thereon. Using MOD16 ET data for XRP from 2000 to 2019, the spatial distribution characteristics of the average ET for that period were statistically analyzed. As shown in Figure 2, the spatial distribution of ET in the study area was significantly different, i.e., it was high in the east and low in the west. The mean values on a spatial scale were 183.26-743.36 mm. ET range was divided into five levels at equal intervals according to ET value distribution, and the proportion and spatial distribution of different levels were statistically analyzed. The low-value areas of ET (183-250 mm) accounted for 13.68% of the total area; these were mostly concentrated in the central and western regions of the study area. The high-value areas (460-743 mm) accounted for only 3.27% of the total area, and were mostly concentrated in the southeast of the study area.
value of ˆk b represents the degree of influence of each factor on ET, and the symb sents the positive or negative effects thereon.

Interannual Change Feature
Using MOD16 ET data for XRP from 2000 to 2019, the spatial distribution ch istics of the average ET for that period were statistically analyzed. As shown in F the spatial distribution of ET in the study area was significantly different, i.e., it w in the east and low in the west. The mean values on a spatial scale were 183.26-743 ET range was divided into five levels at equal intervals according to ET value dist and the proportion and spatial distribution of different levels were statistically a The low-value areas of ET (183-250 mm) accounted for 13.68% of the total area; th mostly concentrated in the central and western regions of the study area. The hig areas (460-743 mm) accounted for only 3.27% of the total area, and were mostly trated in the southeast of the study area. To study the temporal variation characteristics of ET in different land cov we combined ET with data regarding land cover types to statistically analyze the nual variation trends of ET in BAR, SAN, GRA, SWA, FOR, CRO, and PAD from 2019 ( Figure 3). Among the seven land cover types, the annual average ET varied cantly. In PAD it was the largest (421.15 mm), followed by CRO (344.22 mm), FOR mm), SWA (278.89 mm), GRA (271.89 mm), and SAN (246.80 mm). The annual av ET was the lowest in BAR, i.e., only 240.58 mm. The temporal variation characte To study the temporal variation characteristics of ET in different land cover types, we combined ET with data regarding land cover types to statistically analyze the interannual variation trends of ET in BAR, SAN, GRA, SWA, FOR, CRO, and PAD from 2000 to 2019 ( Figure 3). Among the seven land cover types, the annual average ET varied significantly. In PAD it was the largest (421.15 mm), followed by CRO (344.22 mm), FOR (332.84 mm), SWA (278.89 mm), GRA (271.89 mm), and SAN (246.80 mm). The annual average of ET was the lowest in BAR, i.e., only 240.58 mm. The temporal variation characteristics of ET in different land cover types were relatively consistent, showing an overall increasing trend. The peak value appeared in 2013. Among different land cover types, the increasing trend in CRO was the most significant, with a growth rate of 7.60 mm/year, and least significant in BAR, with the growth rate of 4.79 mm/year. ET in different land cover types were relatively consistent, showing an overall increasing trend. The peak value appeared in 2013. Among different land cover types, the increasing trend in CRO was the most significant, with a growth rate of 7.60 mm/year, and least significant in BAR, with the growth rate of 4.79 mm/year. In order to study spatial variation, we analyzed the trend of ET by Theil-Sen Median based on MOD16 ET data from 2000 to 2019 in the XRP, and used the Mann-Kendall trend to test the significance of the observed trends. α was defined as the level of a significance test; α = 0.1, 0.05, 0.01 corresponded to |Z| ≥ 1.645, 1.96, 2.576, respectively, indicating that the study sequence passed a 90%, 95%, and 99% reliability of significance test, respectively. |Z| ≤ 1.645 indicated that the study sequence did not pass the significance test with a reliability of 90%. Regarding spatial distribution from 2000 to 2019 (Figure 4), ET showed an increasing trend in 99.42% of the study area. Only 0.58% of the area showed a decreasing trend; this was concentrated in the central and southeast of the study area. Notably, 93.67% of the area showed an increasing trend that surpassed α = 0.01 significance. In order to study spatial variation, we analyzed the trend of ET by Theil-Sen Median based on MOD16 ET data from 2000 to 2019 in the XRP, and used the Mann-Kendall trend to test the significance of the observed trends. α was defined as the level of a significance test; α = 0.1, 0.05, 0.01 corresponded to |Z| ≥ 1.645, 1.96, 2.576, respectively, indicating that the study sequence passed a 90%, 95%, and 99% reliability of significance test, respectively. |Z| ≤ 1.645 indicated that the study sequence did not pass the significance test with a reliability of 90%. Regarding spatial distribution from 2000 to 2019 (Figure 4), ET showed an increasing trend in 99.42% of the study area. Only 0.58% of the area showed a decreasing trend; this was concentrated in the central and southeast of the study area. Notably, 93.67% of the area showed an increasing trend that surpassed α = 0.01 significance.

Monthly Change Characteristics
Our statistical analysis of monthly ET distribution in the study area was based on MOD16 ET data from 2000 to 2019 ( Figure 5). The monthly mean of ET ranged from 5.8 mm to 148.88 mm, and the spatial distribution had obvious temporal differences. The monthly ET decreased from south to north in January, February, and December. From

Monthly Change Characteristics
Our statistical analysis of monthly ET distribution in the study area was based on MOD16 ET data from 2000 to 2019 ( Figure 5). The monthly mean of ET ranged from 5.8 mm to 148.88 mm, and the spatial distribution had obvious temporal differences. The monthly ET decreased from south to north in January, February, and December. From March to November, the it showed strong distribution in the southeast and low distribution in the northwest. The spatial distribution difference from June to September was particularly obvious. The highest monthly ET that appeared in the southeast of the study area in July was 135.04 mm higher than the lowest monthly ET, which was observed in the north of the study area in the same period.  The monthly ET trends of different land cover types from 2000 to 2019 in the XRP were analyzed. Figure 6 shows a single peak change of ET for all land cover types, as well as a peak appearing in July and August. By comparing the monthly ET, we found that the ET differences between different land cover types in the study area mainly occurred in July and August. The monthly ET of PAD was the highest (94.81 mm) in July, i.e., 59.29 mm higher than that of BAR, the lowest in the same period. Except for July and August, there was no significant difference in the monthly ET between different land cover types in the same month.

Correlation with Meteorological Factors
As an important environmental factor, meteorological elemen all aspects of the ecosystem. WDSP, TEMP, RH, and PRCP are imp factors affecting the inland water cycle [38,63]. In order to explore orological factors on ET, the PCCs between annual ET and WDSP, R calculated using MOD16 ET data and annual raster data regarding in the XRP based on the GIS platform ( Figure 7). The spatial distribu ET and WDSP was significantly different. The areas with significa (PCCs > 0.7), concentrated in the central and southeast of the stud 31.06% of the total. By contrast, in the west of the study area, 15.68

Correlation with Meteorological Factors
As an important environmental factor, meteorological elements are irreplaceable in all aspects of the ecosystem. WDSP, TEMP, RH, and PRCP are important meteorological factors affecting the inland water cycle [38,63]. In order to explore the influence of meteorological factors on ET, the PCCs between annual ET and WDSP, RH, TEMP, PRCP were calculated using MOD16 ET data and annual raster data regarding meteorological factors in the XRP based on the GIS platform ( Figure 7). The spatial distribution of PCCs between ET and WDSP was significantly different. The areas with significant positive correlation (PCCs > 0.7), concentrated in the central and southeast of the study area, accounted for 31.06% of the total. By contrast, in the west of the study area, 15.68% of the total showed negative correlation. ET and RH were mainly negatively correlated, and the correspondent areas with significant negative correlation (PCCs < −0.7) accounted for 12.71% of the study area; these were concentrated in the central and northern areas. The absolute value of the PCCs between ET and TEMP was less than 0.5, showing a weak correlation. The positively correlated areas accounted for 41.64%; these were mainly concentrated in the northwest and east of the study area. ET and PRCP were positively correlated in the study area, i.e., more significantly in the northwest than in the southeast. The areas with significant positive correlation (PCCs > 0.7) accounted for 54.48% of the study area. ent areas with significant negative correlation (PCCs < −0.7) accounted for 12.71% of the study area; these were concentrated in the central and northern areas. The absolute value of the PCCs between ET and TEMP was less than 0.5, showing a weak correlation. The positively correlated areas accounted for 41.64%; these were mainly concentrated in the northwest and east of the study area. ET and PRCP were positively correlated in the study area, i.e., more significantly in the northwest than in the southeast. The areas with significant positive correlation (PCCs > 0.7) accounted for 54.48% of the study area. In order to study the correlation between ET and meteorological factors for different land cover types, the PCCs between the two elements were calculated (Figure 8). The results revealed that the PCCs between ET and WDSP in GRA, CRO, PAD, SWA ranged from 0.06 to 0.97, showing a relatively strong positive correlation. Among them, the positive correlation of PAD was the most significant, with an average PCC of 0.67. However, ET and WDSP showed a negative correlation in more than 75% of BAR areas. The PCCs between ET and RH were negative in GRA, CRO, and SWA, indicating a strong negative correlation. The negative correlation of GRA was the most obvious, with an average PCC of −0.52. The absolute values of the average PCCs between ET and TEMP in all land cover types were less than 0.2, showing no correlation. ET and PRCP showed a significant positive correlation in all land cover types. Except for PAD, the average PCCs between ET and PRCP were above 0.7. In order to study the correlation between ET and meteorological factors for different land cover types, the PCCs between the two elements were calculated (Figure 8). The results revealed that the PCCs between ET and WDSP in GRA, CRO, PAD, SWA ranged from 0.06 to 0.97, showing a relatively strong positive correlation. Among them, the positive correlation of PAD was the most significant, with an average PCC of 0.67. However, ET and WDSP showed a negative correlation in more than 75% of BAR areas. The PCCs between ET and RH were negative in GRA, CRO, and SWA, indicating a strong negative correlation. The negative correlation of GRA was the most obvious, with an average PCC of −0.52.
The absolute values of the average PCCs between ET and TEMP in all land cover types were less than 0.2, showing no correlation. ET and PRCP showed a significant positive correlation in all land cover types. Except for PAD, the average PCCs between ET and PRCP were above 0.7.

Correlation with NDVI
Vegetation transpiration is the main mechanism of ET. It is mainly affected by vegetation species, root depth, and plant moisture content, and varies with vegetation growth and coverage. NDVI is the best indicator of growth status and spatial distribution density of vegetation. It is linearly correlated with vegetation distribution density and can effec-

Correlation with NDVI
Vegetation transpiration is the main mechanism of ET. It is mainly affected by vegetation species, root depth, and plant moisture content, and varies with vegetation growth and coverage. NDVI is the best indicator of growth status and spatial distribution density of vegetation. It is linearly correlated with vegetation distribution density and can effectively reflect the growth status [64,65]. In order to explore the influence of vegetation growth status and coverage on ET, the PCCs between annual ET and NDVI were calculated based on the GIS platform in each pixel. As shown in Figure 9, ET and NDVI were mainly positively correlated in the study area, with an average of 0.71. The areas with significant positive correlation (PCCs > 0.7) account for 76.34% of the study area. The areas with weak correlation (−0.5 < PCCs < 0.5) account for 12.14%, and were mainly concentrated in the south-east and middle-west of the study area. Based on the land cover data of XRP, the PCCs between ET and NDVI in different land cover types were statistically analyzed. As shown in Figure 10, the positive correlation between ET and NDVI was significant in CRO and FOR with higher vegetation density. The PCCs of more than 75% of the areas were over 0.6; the FOR average was the largest, i.e., up to 0.81. However, in BAR and SAN, with lower vegetation density, the average PCCs between ET and NDVI were only about 0.3, showing a relatively weak correlation. In GRA, PAD, and SWA, though the PCCs between ET and NDVI had a wide range, the average values were more than 0.6.

Correlation with Groundwater Depth
ET can be divided into soil moisture evaporation and groundwater evaporation from evaporated water sources. Groundwater evaporation is the process in which groundwater is transported to the gas belt and the atmosphere through soil evaporation and plant transpiration. It is the main mechanism for the transformation of groundwater to soil moisture and atmospheric water [66]; its extent is mainly affected by groundwater depth. According to the situation in the XRP, the groundwater depth was divided into 14 levels at intervals of 1 m. As shown in Figure 11a, the groundwater depths in the study area ranged

Correlation with Groundwater Depth
ET can be divided into soil moisture evaporation and groundwater evaporation from evaporated water sources. Groundwater evaporation is the process in which groundwater is transported to the gas belt and the atmosphere through soil evaporation and plant transpiration. It is the main mechanism for the transformation of groundwater to soil moisture and atmospheric water [66]; its extent is mainly affected by groundwater depth. According to the situation in the XRP, the groundwater depth was divided into 14 levels at inter-

Correlation with Groundwater Depth
ET can be divided into soil moisture evaporation and groundwater evaporation from evaporated water sources. Groundwater evaporation is the process in which groundwater is transported to the gas belt and the atmosphere through soil evaporation and plant transpiration. It is the main mechanism for the transformation of groundwater to soil moisture and atmospheric water [66]; its extent is mainly affected by groundwater depth. According to the situation in the XRP, the groundwater depth was divided into 14 levels at intervals of 1 m. As shown in Figure 11a, the groundwater depths in the study area ranged from 1.14 m to 14.89 m, and were generally deeper in the west and shallower in the east. The areas with groundwater depths of over 6.0 m were mainly concentrated in the northwest and middle-west of the study area. A spatial overlay analysis of ET, land cover types, and groundwater depth data was carried out to calculate the mean ET corresponding to the different groundwater depths in each land cover type. Figure 11b illustrates that ET decreases significantly with the increase in groundwater depth in GRA, CRO, BRA, SAN, and SWA. In FOR and PAD, although ET also decreased with the increase of groundwater depth, the curves had significantly greater fluctuation, indicating the minor effect of groundwater depth on ET in FOR and PAD. In SWA, GRA, and CRO, the influence of groundwater depth on ET was relatively consistent. When the groundwater depth was less than 3.0 m, the curves had no obvious change trend, indicating that the evaporated water in the surface soil could be fully replenished by groundwater through capillary action. Therefore, ET was not affected by groundwater depth. When the groundwater depth was more than 3.0 m, the curves showed a downward trend with the increase of groundwater depth, indicating that with the increase of water transport distance, the water supplied to the surface by capillary action, was reduced, thereby also decreasing ET. When the groundwater depth was 6.0 m, ET was lowest and tended to be stable. Meanwhile, ET was mainly dominated by soil moisture evaporation and vegetation transpiration. Therefore, the limit depths of groundwater evaporation of SWA, GRA, and CRO in the XRP could be determined as 6.0 m. In BAR and SAN, ET also decreased with the increase of groundwater depth and remained steady when the depth reached 4.0 m, indicating that the limit depth of groundwater evaporation in BAR and SAN in the XRP was 4.0 m. These findings imply that vegetation can increase the depth of the groundwater absorbed by land due to the water potential difference caused by transpiration allowing plant roots to absorb water.
Remote Sens. 2022, 13, x 14 of 24 the groundwater depth was 6.0 m, ET was lowest and tended to be stable. Meanwhile, ET was mainly dominated by soil moisture evaporation and vegetation transpiration. Therefore, the limit depths of groundwater evaporation of SWA, GRA, and CRO in the XRP could be determined as 6.0 m. In BAR and SAN, ET also decreased with the increase of groundwater depth and remained steady when the depth reached 4.0 m, indicating that the limit depth of groundwater evaporation in BAR and SAN in the XRP was 4.0 m. These findings imply that vegetation can increase the depth of the groundwater absorbed by land due to the water potential difference caused by transpiration allowing plant roots to absorb water.

Correlation with Topography Factors
Elevation and slope are basic topographical characteristics and the main factors affecting species, growth states, and the spatial distribution of vegetation. Their change alters soil moisture and surface runoff, thereby influencing ET. An elevation distribution map of XRP was made using ASTER DEM data. According to their distribution character-

Correlation with Topography Factors
Elevation and slope are basic topographical characteristics and the main factors affecting species, growth states, and the spatial distribution of vegetation. Their change alters soil moisture and surface runoff, thereby influencing ET. An elevation distribution map of XRP was made using ASTER DEM data. According to their distribution characteristics, the elevations were divided into 10 categories (0-100 m, 100-200 m, 200-300 m, 300-400 m, 400-500 m, 500-600 m, 600-700 m, 700-800 m, 800-900 m, and 900-1000 m). As shown in Figure 12a, the geography of the XRP is low and flat, tilting from the southwest to the northeast, with an average elevation of 253.94 m above sea level. A spatial overlay analysis of ET, land cover types, and elevation data in the study area was carried out to calculate the mean ET corresponding to different elevation levels for each land cover type. As shown in Figure 12b, except for PAD, the elevation effects on ET were relatively consistent. ET decreased obviously with the increase of elevation and tended to be stable when the elevation reached 300 m, indicating that elevation has a negative influence on ET when below 300 m. Although the ET of PAD showed a decreasing trend with the increase of elevation, the fluctuations in the curve indicate only slight effects. These phenomena may be related to surface runoff. After rainfall, rainwater is more likely to collect at low altitude areas, providing more water sources for ET. Affected by artificial irrigation, PAD has sufficient water, so it is not affected. Based on ASTER DEM data, a slope distribution map of the XRP (Figure 13a) was made. According to the classification standard of landform slopes by the International Geographical Union, the slope of the XRP was divided into six categories (0-0.5°, 0.5-2°, 2-5°, 5-15°, 15-35°, and 35-55°). The average ET of the different slope levels for each land cover type was calculated. As shown in Figure 13b, except for PAD, the slope effects on ET were relatively consistent. When the slope was less than 15°, ET changed slightly, and remained stable with the increase of slope. When the slope was greater than 15°, ET increased obviously with the increase of slope. This may have been caused by the relatively fast flow rate of rainwater in areas with large slopes. The larger the slope, the faster the runoff, and thus, the greater the ET. However, due to the low adhesion and loose texture of sand, SAN had a strong absorption capacity for rainwater, preventing it from forming surface runoff, even in areas with large slopes. Therefore, slope had no significant effect on ET in SAN. This implies that the influence of topographic conditions on ET may be mainly reflected in the differences in water conditions caused by surface runoff. Based on ASTER DEM data, a slope distribution map of the XRP (Figure 13a) was made. According to the classification standard of landform slopes by the International Geographical Union, the slope of the XRP was divided into six categories (0-0.5 • , 0.5-2 • , 2-5 • , 5-15 • , 15-35 • , and 35-55 • ). The average ET of the different slope levels for each land cover type was calculated. As shown in Figure 13b, except for PAD, the slope effects on ET were relatively consistent. When the slope was less than 15 • , ET changed slightly, and remained stable with the increase of slope. When the slope was greater than 15 • , ET increased obviously with the increase of slope. This may have been caused by the relatively fast flow rate of rainwater in areas with large slopes. The larger the slope, the faster the runoff, and thus, the greater the ET. However, due to the low adhesion and loose texture of sand, SAN had a strong absorption capacity for rainwater, preventing it from forming surface runoff, even in areas with large slopes. Therefore, slope had no significant effect on ET in SAN. This implies that the influence of topographic conditions on ET may be mainly reflected in the differences in water conditions caused by surface runoff.

The Relative Contribution Rate of Influencing Factors to ET
ET and its influencing factors showed different correlations for different land cover types. In this study, we used the raster data of annual average ET (mm), WDSP (m·s −1 ), TEMP (°C), PRCP (mm), RH (%), NDVI, groundwater depth (m), meters above sea level (m), and slope (°) in the study area from 2000 to 2019 to create a multisource spatial data set. Using this, we conducted spatial balance sampling in different land cover types at 5km intervals. Then, we took the WDSP, TEMP, PRCP, RH, NDVI, groundwater depth, elevation, and slope as the inputs (Xi) for the training sample and the annual ET value as the output (Yi) to construct a training sample set (Xi, Yi) in order to construct the statistical regression analysis model to quantitatively analyze the influence of various factors on ET. Independent variables showed multicollinearity because of their correlations. Linear regression can lead to overfitting of the regression equation, causing the independent variable coefficient to fail to conform to the actual meaning of the influencing factors. Therefore, it was necessary to diagnose the collinearity of each influencing factor before constructing the regression model. The variance inflation factor (VIF) was generally used as the standard of collinearity diagnoses. If VIF was greater than 10, it indicated that the model had serious multicollinearity. The collinearity diagnostic results (Figure 14a) showed that the VIFs of PRCP, RH and, NDVI were greater than 10, indicating significant multicollinearity. In order to weaken the influence of multicollinearity in independent variables on the regression coefficient of the model, this study used ridge regression analysis to achieve a more realistic regression result at the expense of rejecting some information and reducing the accuracy of the model [62].

The Relative Contribution Rate of Influencing Factors to ET
ET and its influencing factors showed different correlations for different land cover types. In this study, we used the raster data of annual average ET (mm), WDSP (m·s −1 ), TEMP ( • C), PRCP (mm), RH (%), NDVI, groundwater depth (m), meters above sea level (m), and slope ( • ) in the study area from 2000 to 2019 to create a multisource spatial data set. Using this, we conducted spatial balance sampling in different land cover types at 5-km intervals. Then, we took the WDSP, TEMP, PRCP, RH, NDVI, groundwater depth, elevation, and slope as the inputs (Xi) for the training sample and the annual ET value as the output (Yi) to construct a training sample set (Xi, Yi) in order to construct the statistical regression analysis model to quantitatively analyze the influence of various factors on ET. Independent variables showed multicollinearity because of their correlations. Linear regression can lead to overfitting of the regression equation, causing the independent variable coefficient to fail to conform to the actual meaning of the influencing factors. Therefore, it was necessary to diagnose the collinearity of each influencing factor before constructing the regression model. The variance inflation factor (VIF) was generally used as the standard of collinearity diagnoses. If VIF was greater than 10, it indicated that the model had serious multicollinearity. The collinearity diagnostic results (Figure 14a) showed that the VIFs of PRCP, RH and, NDVI were greater than 10, indicating significant multicollinearity. In order to weaken the influence of multicollinearity in independent variables on the regression coefficient of the model, this study used ridge regression analysis to achieve a more realistic regression result at the expense of rejecting some information and reducing the accuracy of the model [62]. From the results of our ridge regression analysis (Figure 14b, Table 1), the multicollinearity among the influencing factors was shown to be significantly weakened (VIF < 10), and the model determination coefficient (R 2 ) reached more than 0.6 for all land cover types. PRCP, NDVI, and slope had positive effects on ET for all land cover types, while RH, groundwater depth, and elevation had negative effects. It is noteworthy that TEMP had a negative effect on ET in GRA, BAR, and SAN, while WDSP had a negative effect on ET in BAR and SAN, and TEMP and WDSP had a positive effect on ET in other land cover types. We calculated the relative contribution rate of each influencing factor to ET in different land cover types according to the ridge regression analysis results. In Figure 15, the NDVI, RH, PRCP, and WDSP represent a cumulative contribution rate of over 60% to ET in all land cover types, thereby serving as the main influencing factors. The NDVI, RH, PRCP, WDSP, and groundwater depth in GRA, CRO, and SWA represent a cumulative contribution rate of over 80%, comprising the best explanation for changes in ET. In FOR, NDVI shows the largest relative contribution rate, i.e., 35.09%, making it the largest affecting factor, followed by RH, PRCP, WDSP, slope, and groundwater depth. NDVI is also the leading determinant of ET in PAD, although the influence of WDSP is greater than that of PRCP and RH. In BAR and SAN, PRCP is the leading determinant of ET, followed by NDVI, WDSP, and RH. From the results of our ridge regression analysis (Figure 14b, Table 1), the multicollinearity among the influencing factors was shown to be significantly weakened (VIF < 10), and the model determination coefficient (R 2 ) reached more than 0.6 for all land cover types. PRCP, NDVI, and slope had positive effects on ET for all land cover types, while RH, groundwater depth, and elevation had negative effects. It is noteworthy that TEMP had a negative effect on ET in GRA, BAR, and SAN, while WDSP had a negative effect on ET in BAR and SAN, and TEMP and WDSP had a positive effect on ET in other land cover types. We calculated the relative contribution rate of each influencing factor to ET in different land cover types according to the ridge regression analysis results. In Figure 15, the NDVI, RH, PRCP, and WDSP represent a cumulative contribution rate of over 60% to ET in all land cover types, thereby serving as the main influencing factors. The NDVI, RH, PRCP, WDSP, and groundwater depth in GRA, CRO, and SWA represent a cumulative contribution rate of over 80%, comprising the best explanation for changes in ET. In FOR, NDVI shows the largest relative contribution rate, i.e., 35.09%, making it the largest affecting factor, followed by RH, PRCP, WDSP, slope, and groundwater depth. NDVI is also the leading determinant of ET in PAD, although the influence of WDSP is greater than that of PRCP and RH. In BAR and SAN, PRCP is the leading determinant of ET, followed by NDVI, WDSP, and RH.

Reason Analysis of Temporal Changes in ET
This study found that among the various factors, NDVI, RH, PRCP, and WDSP wer the leading determinants of interannual ET changes for all land cover types. Through comparative analysis of the interannual trends of ET, NDVI, RH, PRCP, and WDSP in th XRP from 2000 to 2019 (Figure 16), ET changes were divided into three periods: the firs was a rising period from 2000 to 2005, the second another rising period from 2006 to 2013 and the third a stable period from 2014 to 2019. The first period represented the mos significant ET increase, with a 7.06% average growth rate. During this period, only RH showed a decreasing trend, while NDVI, WDSP, and PRCP increased. These changes ben efited the increase of ET. The increase of NDVI was the most significant, with an averag annual growth rate of 3.57%; this may have been the main reason for the increase in ET In the second period, ET showed a relatively smooth increase, with an average growth rate of 3.86%. NDVI represented a relatively smooth increasing trend, while RH basically showed a consistent increase relative to the previous period. The significant increase o PRCP and WDSP may have been the main driving forces of ET increase during this period In 2014-2019, NDVI, RH, PRCP, and WDSP tended to be stable, and the annual averag ET remained unchanged at about 370 mm. The results showed that the change of climati factors had a significant influence on ET, consistent with results for most basins in China such as the Pearl River and Tao River Basins. However, the difference lay in the fact tha TEMP is the main climatic factor affecting ET in the aforementioned basins [43,67]. In thi study, RH, PRCP, and WDSP were found to have a large influence on ET. This may hav been related to climate variability and the distribution of water resources in different re gions [68][69][70][71]. Combined with the results of this study, TEMP and other energy condition are the main drivers for ET change in humid or subhumid climate areas where water re sources are relatively abundant. In contrast, in arid climates with relatively scarce wate resources, water conditions such as vegetation and PRCP can further determine change in ET. This further explains why changes in ET are subject to multiple influencing factors and why the main factors vary in different periods or under different climatic condition [41][42][43][44]72].

Reason Analysis of Temporal Changes in ET
This study found that among the various factors, NDVI, RH, PRCP, and WDSP were the leading determinants of interannual ET changes for all land cover types. Through a comparative analysis of the interannual trends of ET, NDVI, RH, PRCP, and WDSP in the XRP from 2000 to 2019 (Figure 16), ET changes were divided into three periods: the first was a rising period from 2000 to 2005, the second another rising period from 2006 to 2013, and the third a stable period from 2014 to 2019. The first period represented the most significant ET increase, with a 7.06% average growth rate. During this period, only RH showed a decreasing trend, while NDVI, WDSP, and PRCP increased. These changes benefited the increase of ET. The increase of NDVI was the most significant, with an average annual growth rate of 3.57%; this may have been the main reason for the increase in ET. In the second period, ET showed a relatively smooth increase, with an average growth rate of 3.86%. NDVI represented a relatively smooth increasing trend, while RH basically showed a consistent increase relative to the previous period. The significant increase of PRCP and WDSP may have been the main driving forces of ET increase during this period. In 2014-2019, NDVI, RH, PRCP, and WDSP tended to be stable, and the annual average ET remained unchanged at about 370 mm. The results showed that the change of climatic factors had a significant influence on ET, consistent with results for most basins in China, such as the Pearl River and Tao River Basins. However, the difference lay in the fact that TEMP is the main climatic factor affecting ET in the aforementioned basins [43,67]. In this study, RH, PRCP, and WDSP were found to have a large influence on ET. This may have been related to climate variability and the distribution of water resources in different regions [68][69][70][71]. Combined with the results of this study, TEMP and other energy conditions are the main drivers for ET change in humid or subhumid climate areas where water resources are relatively abundant. In contrast, in arid climates with relatively scarce water resources, water conditions such as vegetation and PRCP can further determine changes in ET. This further explains why changes in ET are subject to multiple influencing factors, and why the main factors vary in different periods or under different climatic conditions [41][42][43][44]72].

Analysis of ET Differences for Different Land Cover Types
The factors affecting ET can be divided into three categories: water conditions conditions, and dynamic conditions. The water conditions describe the amount output on the evaporation surface. The energy conditions indicate the speed of wa version from liquid to gas, and the dynamic conditions indicate the speed of wate transport over the evaporation surface. As the main parameter to reflect pressure o resources, the water stress index (WSI) is closely related to water resource quant In order to further explore the reasons behind the obvious differences in ET, we u to analyze the influence of water conditions on ET for different land cover ty shown in Figure 17, the WSI for different land cover types varied significantly; ET decreased with the increase of WSI. In PAD and CRO, artificial irrigation ensu ficient water resources, i.e., a smaller WSI corresponds to greater ET. For other lan types, PRCP and capillary action (which transports groundwater to the surface) w main sources of surface water. Meanwhile, different species and densities of ve were key factors influencing the absorption of groundwater and the maintenanc water content in surface soil. The ultimate depths after absorbing groundwater SWA, and GRA were greater than those of SAN and BAR. The groundwater tran more water to the surface through capillary action. Moreover, the vegetation cov FOR, SWA, and GRA was greater. The interception of canopies and leaves during lessens rainwater infiltration and increases the surface soil moisture content. Th

Analysis of ET Differences for Different Land Cover Types
The factors affecting ET can be divided into three categories: water conditions, energy conditions, and dynamic conditions. The water conditions describe the amount of water output on the evaporation surface. The energy conditions indicate the speed of water conversion from liquid to gas, and the dynamic conditions indicate the speed of water vapor transport over the evaporation surface. As the main parameter to reflect pressure on water resources, the water stress index (WSI) is closely related to water resource quantity [73]. In order to further explore the reasons behind the obvious differences in ET, we used WSI to analyze the influence of water conditions on ET for different land cover types. As shown in Figure 17, the WSI for different land cover types varied significantly; notably, ET decreased with the increase of WSI. In PAD and CRO, artificial irrigation ensures sufficient water resources, i.e., a smaller WSI corresponds to greater ET. For other land cover types, PRCP and capillary action (which transports groundwater to the surface) were the main sources of surface water. Meanwhile, different species and densities of vegetation were key factors influencing the absorption of groundwater and the maintenance of the water content in surface soil. The ultimate depths after absorbing groundwater in FOR, SWA, and GRA were greater than those of SAN and BAR. The groundwater transported more water to the surface through capillary action. Moreover, the vegetation coverage of FOR, SWA, and GRA was greater. The interception of canopies and leaves during rainfall lessens rainwater infiltration and increases the surface soil moisture content. Therefore, the WSI of FOR, SWA, and GRA was smaller than that of SAN and BAR; as such, the former three types demonstrated greater ET. This indicated that the water conditions in high-latitude, arid areas significantly affects ET in various land cover types. The difference of ET among different land covers was not only due to the different water conditions, but also to the different mechanisms of action and degrees of the energy conditions and dynamic conditions required for evaporation in different land covers. TEMP, as the main energy condition, showed a different correlation with ET in different land covers. An increase of TEMP not only accelerated water evaporation in the soil, but also controlled the stomatal opening of plants, resulting in changes of vegetation evapotranspiration. This indicated that the influence of TEMP on ET was different in areas with different vegetation coverage. As the main dynamic conditions, the effects of RH and WDSP on ET were obviously different for different land covers. Only in SAN and BAR soil with serious water stress did WDSP have a negative effect. This may be because in areas with low water stress, WDSP accelerates the rising speed of water vapor above the evaporation surface, thus accelerating the rate of ET, while in SAN and BAR (without vegetation cover), the surface is gradually weathered, resulting in a negative impact. Therefore, studying the temporal and spatial changes of ET and its influencing factors from the perspective of the water cycle characteristics of different land cover types is necessary for understanding variability in the hydrologic cycle.
Remote Sens. 2022, 13, x energy condition, showed a different correlation with ET in different land cov crease of TEMP not only accelerated water evaporation in the soil, but also co stomatal opening of plants, resulting in changes of vegetation evapotranspi indicated that the influence of TEMP on ET was different in areas with differen coverage. As the main dynamic conditions, the effects of RH and WDSP on ET ously different for different land covers. Only in SAN and BAR soil with se stress did WDSP have a negative effect. This may be because in areas with stress, WDSP accelerates the rising speed of water vapor above the evaporat thus accelerating the rate of ET, while in SAN and BAR (without vegetation surface is gradually weathered, resulting in a negative impact. Therefore, st temporal and spatial changes of ET and its influencing factors from the perspe water cycle characteristics of different land cover types is necessary for und variability in the hydrologic cycle.
1 ET WSI PET = − Figure 17. ET and water stress index for different land covers.

Conclusions
The differences among water evaporation mechanisms and action condi ferent land cover types result in varying ET for different land cover types. T ences are reflected not only in the changing trends, but also in the main factors ET. We analyzed the spatiotemporal variation of ET for different land cover ty correlations with meteorological factors, NDVI, groundwater depth, elevation and applied the ridge regression model to quantify the contributions of these the past 20 years, due to vegetation greening and increasing precipitation, ET increasing trend at a rate of 6.72 mm/year. This trend was most obvious in CRO in additional pressure on agricultural irrigation. Therefore, we should pay mo to the soil moisture content in CRO and carry out agricultural irrigation in a t

Conclusions
The differences among water evaporation mechanisms and action conditions in different land cover types result in varying ET for different land cover types. These differences are reflected not only in the changing trends, but also in the main factors influencing ET. We analyzed the spatiotemporal variation of ET for different land cover types and its correlations with meteorological factors, NDVI, groundwater depth, elevation, and slope, and applied the ridge regression model to quantify the contributions of these factors. In the past 20 years, due to vegetation greening and increasing precipitation, ET showed an increasing trend at a rate of 6.72 mm/year. This trend was most obvious in CRO, resulting in additional pressure on agricultural irrigation. Therefore, we should pay more attention to the soil moisture content in CRO and carry out agricultural irrigation in a timely manner, in order to ensure that appropriate amounts of water are being delivered for crop growth. Although the correlations among ET and the aforementioned factors were quite different for different land cover types, NDVI, RH, PRCP, and WDSP were shown to have a cumulative contribution rate of more than 60% in all land cover types, and could serve as the main determinants of ET. In GRA, CRO, PAD, FOR and SWA, NDVI had the greatest influence, while PRCP was the dominant influencing factor for BAR and SAN. Therefore, policymakers should apply targeted measures for water resources management for different land cover types. Especially in SAN and BAR with severe water stress, afforestation should be increased to maintain soil water availability. Although there were some limitations in this study, i.e., the differences of ecological and climatic environments in different regions, our results emphasize the importance of comprehensively considering the differences among various hydrologic cycles according to land cover type to assess the contributions of influencing factors to ET. In the future, we will conduct further research on the hydrologic cycle in humid and semihumid regions.