Quantifying the Responses of Evapotranspiration and Its Components to Vegetation Restoration and Climate Change on the Loess Plateau of China

Quantitatively identifying the influences of vegetation restoration (VR) on water resources is crucial to ecological planning. Although vegetation coverage has improved on the Loess Plateau (LP) of China since the implementation of VR policy, the way vegetation dynamics influences regional evapotranspiration (ET) remains controversial. In this study, we first investigate long-term spatiotemporal trends of total ET (TET) components, including ground evaporation (GE) and canopy ET (CET, sum of canopy interception and canopy transpiration) based on the GLEAM-ET dataset. The ET changes are attributed to VR on the LP from 2000 to 2015 and these results are quantitatively evaluated here using the Community Land Model (CLM). Finally, the relative contributions of VR and climate change to ET are identified by combining climate scenarios and VR scenarios. The results show that the positive effect of VR on CET is offset by the negative effect of VR on GE, which results in a weak variation in TET at an annual scale and an increased TET is only shown in summer. Regardless of the representative concentration pathway (RCP4.5 or RCP8.5), differences resulted from the responses of TET to different vegetation conditions ranging from −3.7 to −1.2 mm, while climate change from RCP4.5 to RCP8.5 caused an increase in TET ranging from 0.1 to 65.3 mm. These findings imply that climate change might play a dominant role in ET variability on the LP, and this work emphasizes the importance of comprehensively considering the interactions among climate factors to assess the relative contributions of VR and climate change to ET.


Introduction
Evapotranspiration (ET) is an essential component of both hydrological processes and surface energy exchange, which regulates water allocation to a land surface and thus plays an important role in water resource management [1,2]. Over the past three decades, dynamic changes in ET have been investigated extensively, revealing that global annual ET generally exhibits an increasing trend, but the magnitude of the increase varies between different regions [3,4]. Although it has been demonstrated that the local climate is a key driver of ET, land cover changes also play an important role in ET changes as the ET process can be modulated by vegetation cover by altering canopy interception, canopy transpiration, and ground evaporation (GE) [5][6][7][8][9]. Quantitatively identifying the causes of ET change is crucial to better understand hydrological process change due to climatic and anthropogenic impacts and to support decision makers in water resource planning and management.
In the context of climate change, many studies have tried to explain ET change [10,11] and the responses of ET to various climatic factors, such as air temperature, solar radiation, precipitation, relative humidity, and wind speed, and these factors have been extensively investigated through observations and modeling; however, inconsistent conclusions have been drawn due to different climate zones, vegetation cover characteristics, and land management practices. For example, Guo et al. [12] assessed the sensitivity of ET to climate conditions based on 30 Australian locations in different climatic zones and concluded that wind and relative humidity are important variables for dry and humid catchments, respectively, while solar radiation has the greatest influence in warmer catchments. Wang et al. [13] found that wind speed was the most influential climatic variable for ET variability in China during 1961−2013, followed by the maximum daily temperature and sunlight duration. Pour et al. [14] reported that minimum temperature was the most influential factor for ET increase in peninsular Malaysia. Additionally, the rate of evaporation from open pans of water in the Northern Hemisphere has been steadily decreasing in past decades, which is in contrast with the expectation that a warming climate should cause an increase in the rate of evaporation from terrestrial open water bodies [15]. These findings reveal the complexity of the interaction between climate variables and ET and the uncertainty of existing knowledge.
The influences of land cover change on ET have also attracted much attention. There is consensus regarding the responses of ET to ecological restoration, i.e., afforestation is generally effective in improving land surface greenness and simultaneously decreases water yields and increases ET at the watershed scale [16][17][18][19]. In contrast, vegetation degradation exerts the opposite effects and leads to decreased ET [20][21][22][23]; however, the consequences of the restoration on ET at regional scales and in large river basins are still contentious [18,24,25], because different land cover types have different biophysical effects due to different leaf area index (LAI), surface albedo, and root depth characteristics [16]. Moreover, stomatal responses to increased atmospheric CO 2 concentrations vary in vegetation species, which correspondingly have different evaporation rates [26,27].
Identifying the contribution of climate change and vegetation coverage to ET change can help understand the causes of hydrological variations and explain the complex interactions between vegetation and climate. For this purpose, a number of attribution analyses have been conducted globally [1,21,[28][29][30], and several types of approaches have been used to investigate the linkage of ET and various driving variables. Field observations are commonly used to quantify the impacts of climate and vegetation on ET, such as eddy covariance techniques, porometry, and lysimeters, and scintillometry, but these observation methods only cover small areas at the site level with short time spans; thus, the response mechanisms depending on the multiple vegetation cover and long-term climate evolution at the regional scale cannot be completely reflected. In recent years, remote sensing has been used to explore ET changes over a large spatial scale. An increasing number of ET products have been developed using land surface models and data assimilation based on remote sensing. These products have been largely used to detect ET trends [4] and statistically regress and correlate ET and environmental factors, such as land cover types, surface temperature, precipitation, and wind speed [6,19,31,32]. The correlation values cannot quantitatively reflect the contribution of certain factors to ET, and assuming linear relationships between ET and climatic factors is inconsistent with reality because the relationship is nonlinear [33].
With the advancement of computer techniques and the understanding of land and atmosphere coupling, large-scale hydrological models [17,34], empirical models based on the Budyko framework [35], remote sensing-or energy balance-based ET models [36,37] and process-driven land surface models [38,39] have been created and applied to explore how environmental factors affect ET. In particularly, land surface models, as offline hydrologic models or land surface modules in tandem with Earth system models, have become powerful tools to study regional and global water cycles and climate change under heterogonous underlying surface conditions, because they have comprehensive biogeophysical parameterization mechanisms to capture the physical and biological interactions between land and atmosphere [40]. The Chinese Loess Plateau (LP), located in the transition zone between the southeast humid monsoon climate and northwest inland arid climate, covers an area of 6.4 × 10 5 km 2 [41]. The plateau has experienced severe vegetation degradation and soil erosion. To conserve soil and water resources in this region, a series of ecological restoration programs have been implemented since the 1960s, leading to increased vegetation coverage [41]. In particular, the largest program, called the Green for Grain Project, was launched in 1999, in which a large proportion of cropland was converted to forest and grassland. These actions largely modify the landscape and simultaneously caused decreases in both water discharge and sediment load [42]. Inevitably, large scale revegetation consumes more water from soil to maintain growth, accompanied by water loss through canopy transpiration [43]. Consequently, elucidating whether current vegetation restoration (VR) results in more water loss through ET in the LP has attracted the attention of many scholars [35,[44][45][46][47]; however, the quantitative analysis of the contribution of VR and climate change to ET remains limited. In this study, we quantitatively evaluate the contribution of VR and climate change to the ET on the LP at the regional scale based on remote sensing products and the Community Land Model (CLM). The specific objectives here are to (1) detect the spatiotemporal change in the ET component, including GE and canopy ET (CET) during 1980−2018, (2) evaluate the CLM performance in simulating ET change in the LP region, and (3) quantitatively differentiate the contributions of climate change and VR to ET based on sensitivity simulations. The results can not only help us gain a better understanding of the hydrological effect of VR, but also provide decision makers with accurate and timely information for restoration planning and water resource management.

Model Description
All simulations used version 4.0 of the CLM (CLM4.0, hereafter referred to as CLM), which is coupled in Community Earth System Model version 1.2 (CESM1.2) and has been widely used in different ecosystems to evaluate the effects of land cover change, atmospheric deposition, and climate change on soil and surface waters in terrestrial ecosystems. The model parameterizes solar radiation partitioning, water transfer between the atmosphere and land surface, and soil organic matter dynamics. The CLM is designed to represent complex vegetation structures by applying plant functional types (PFTs), and seventeen PFTs are defined to capture plant physiology and ecological function, including needleleaf evergreen temperate trees, needleleaf evergreen boreal trees, needleleaf deciduous boreal trees, broadleaf evergreen tropical trees, broadleaf evergreen temperate trees, broadleaf deciduous tropical trees, broadleaf deciduous temperate trees, broadleaf deciduous boreal trees, C3 arctic grasses, C3 non-arctic grasses, C4 grasses, broadleaf evergreen shrubs, broadleaf deciduous temperate shrubs, broadleaf deciduous boreal shrubs, corn, wheat, and bare ground. The CLM allows multiple PFTs to coexist in each grid cell and land cover heterogeneity is specified through the foliage projective cover and LAI of each PFT. The soil profiles are divided into 15 layers with different depths from 0.007 m to 35.178 m in each grid cell. More information about the CLM and its theory can be found in the theoretical document [48].

Data Sources
The ET dataset, based on the Global Land Evaporation Amsterdam Model (GLEAM), was used to detect the spatiotemporal change in ET in the LP region during the period of 1980−2018. The GLEAM encompasses a set of algorithms which use satellite forcing data to separately estimate ET components, including GE, canopy transpiration, canopy evaporation, snow sublimation, and open-water evaporation [49]. Actual ET in the GLEAM is expressed as a function of potential ET and a stress factor that accounts for water limitations, and the potential ET is calculated based on the Priestley and Taylor equation [50]. The newest version, GLEAM-ET 3.3a, is provided with a 0.25 • × 0.25 • latitude-longitude grid with a monthly temporal resolution (https://www.gleam.eu/; accessed on: 6 January Remote Sens. 2021, 13, 2358 4 of 18 2021). The GLEAM dataset has been widely used and verified to reflect the spatiotemporal characteristics of different terrestrial evaporation components. The GLEAM-ET dataset for the period from 1995-2004 was compared with ET components simulated by the CLM to evaluate the model performance in the LP region. The necessary atmospheric forcing data for driving the CLM were developed by Qian et al. [51], and these data have been widely tested and successfully used in different regions to explore the water cycle. The vegetation cover represented by PFT on the LP in 2000 and 2015 was acquired from satellite observation-based global land cover products at a 300 m resolution from the European Space Agency (ESA) as part of the Climate Change Initiative (CCI) land cover project (http://www.esa-landcover-cci.org/; accessed on: 16 September 2020). The LAI data, obtained from the eight-day composite Global Land Surface Satellite (GLASS) product with a spatial resolution of 1 km, were retrieved from the Moderate Resolution Imaging Spectroradiometer (MODIS) reflectance product (MOD09A1) [52] (downloaded from the National Earth System Science Data Center, http://www.geodata.cn; accessed on: 9 September 2020). The LAI datasets were combined with the vegetation cover dataset to obtain the monthly LAI of each PFT using spatial overlay analysis.

Experimental Designs
Our analyses of the impact of VR on ET dynamics are based on two modern simulations that are performed using an offline CLM with fixed atmospheric forcing. The first simulation was designed for the control experiment (LC2000), representing vegetation cover in the LP region at the 2000 level. The model was forced with Qian atmospheric forcing data during 1972-2004, and the greenhouse gases were set at the 2000 level, with a CO 2 concentration of 367 ppm. Then, the model was run for 99 years with three replicates of forcing data. The first two cycles were simulated for 66 years as spin-up periods to bring the carbon and nitrogen states of the land surface into dynamic equilibrium, and the third cycle was simulated for 33 years for analysis. The second simulation was a sensitivity simulation (LC2015) which considered the reality of VR in the LP region. The model set is similar to the control simulations, except that land surface conditions associated with vegetation distribution were substituted with the 2015 level. Finally, the differences were quantified by comparing two simulations to estimate the magnitude and direction of ET change due to VR.
To explore the impact of climate change on ET, high-frequency weather data under future representative concentration pathways (RCPs) of greenhouse gases were needed to use the CLM. Thus, the fully coupled CESM was run under a medium climate stabilization scenario (RCP4.5, stabilized radiative forcing at 4.5 W/m 2 in 2100) and high emission scenario (RCP8.5, increasing radiative forcing up to 8.5 W/m 2 in 2100) to generate 3-hourly weather data from 2005 to 2100 [53], including solar radiation, precipitation, humidity, wind, surface air pressure, and near-surface air temperature. We conducted four simulations that were designed by combining two future climate scenarios (RCP4.5 and RCP8.5) and two vegetation cover scenarios (LC2000 and LC2015).

Data Analysis and Model Evaluation
In this study, to focus on the contribution of vegetation to ET, we mainly analyzed the responses of GE, CET, and total ET (TET) to VR and climate change. CETs were defined as the sum of canopy evaporation and canopy transpiration, and TET was calculated as the sum of GE and CET. The spatiotemporal changes in GE, CET and TET during 1980-2018 were detected using linear regression analysis (least squares method), with a statistical significance analysis based on a t-test. To evaluate the model performance on the ET simulation, the CLM-simulated GE, CET and TET were compared with GLEAM-ET data sets from 1995 to 2004, and the statistical criteria, including the determination coefficient (R 2 ) and root-mean-square error (RMSE), were calculated to indicate how well the simulated ET fit the GLEAM data sets. R 2 mainly indicates the linear association, while the RMSE can reveal the deviations between different datasets [54]. To simplify the spatiotemporal information of ET under the combination of different vegetation coverages and future climates, empirical orthogonal function (EOF) analysis was used to extract the main spatiotemporal patterns of annual ET that were ranked based on their representations of data variance during the period of 1921-2100, and each pattern was associated with a series of time coefficients that described the time evolution of the particular spatial mode [55]. fined as the sum of canopy evaporation and canopy transpiration, and TET was calculated as the sum of GE and CET. The spatiotemporal changes in GE, CET and TET during 1980-2018 were detected using linear regression analysis (least squares method), with a statistical significance analysis based on a t-test. To evaluate the model performance on the ET simulation, the CLM-simulated GE, CET and TET were compared with GLEAM-ET data sets from 1995 to 2004, and the statistical criteria, including the determination coefficient (R 2 ) and root-mean-square error (RMSE), were calculated to indicate how well the simulated ET fit the GLEAM data sets. R 2 mainly indicates the linear association, while the RMSE can reveal the deviations between different datasets [54]. To simplify the spatiotemporal information of ET under the combination of different vegetation coverages and future climates, empirical orthogonal function (EOF) analysis was used to extract the main spatiotemporal patterns of annual ET that were ranked based on their representations of data variance during the period of 1921-2100, and each pattern was associated with a series of time coefficients that described the time evolution of the particular spatial mode [55].   increased (Figure 3f). Although the southern LP had a higher annual average CET, the increasing trend was not statistically significant, even showing a decreasing trend in some areas. Similarly, the area-averaged TET increased significantly with a gradient of 1.58 mm/year (p < 0.01) (Figure 3g), even though the increased CET was offset by a decreased GE. Overall, the spatiotemporal patterns of CET can reflect the corresponding pattern of TET (Figure 3h,i).

Model Validation
The CLM performance in the ET simulation was evaluated by comparing the simulated GE, CET, and TET values with GLEAM datasets from 1995 to 2004 ( Figure 4). In general, the CLM performed well in capturing the dominant spatial pattern of annual average GE, CET and TET, showing a high magnitude in the southeastern LP and low magnitude in the northwestern LP; however, there were overestimations for GE (Figure 4a), which were clearly observed in the area-averaged monthly time series (Figure 4g). The simulated CET and TET had consistent fluctuating trends with the GLEAM-based monthly time series (Figure 4h,i). To further evaluate the performances quantitatively, the R 2 and RMSE between simulated and GLEAM-based area-averaged monthly series were calculated. The obtained high R 2 (ranging from 0.89 to 0.92) and low RMSE (ranging from 4.72 to 7.73) indicated good performance for the CLM in reflecting the spatiotemporal pattern of ET on the LP. Remote Sens. 2021, 13, x FOR PEER REVIEW 7 of 18

Model Validation
The CLM performance in the ET simulation was evaluated by comparing the simulated GE, CET, and TET values with GLEAM datasets from 1995 to 2004 ( Figure 4). In general, the CLM performed well in capturing the dominant spatial pattern of annual average GE, CET and TET, showing a high magnitude in the southeastern LP and low magnitude in the northwestern LP; however, there were overestimations for GE (Figure 4a), which were clearly observed in the area-averaged monthly time series (Figure 4g). The simulated CET and TET had consistent fluctuating trends with the GLEAM-based monthly time series (Figure 4h,i). To further evaluate the performances quantitatively, the R 2 and RMSE between simulated and GLEAM-based area-averaged monthly series were calculated. The obtained high R 2 (ranging from 0.89 to 0.92) and low RMSE (ranging from 4.72 to 7.73) indicated good performance for the CLM in reflecting the spatiotemporal pattern of ET on the LP.  Figure 5 depicts a comparison of the area-averaged annual ET for the two simulations. The simulated GE showed clear differences in their annual time series, and the GE of LC2015 was higher than that of LC2000; however, the CET showed opposite responses to vegetation coverage change, with a larger CET under LC2015 in comparison with LC2000. Consequently, the simulated TET did not show evident differences between LC2000 and LC2015, because the positive effect of VR on CET was offset by the negative effect of VR on GE. The spatial variation in annual average ET associated with VR is illustrated in Figure 6. The whole LP region, except for the small area of the central part and northern edge, showed decreased GE due to vegetation coverage changes from LC2000 to LC2015, particularly in the eastern LP, where GE decreased by more than 15 mm, with a statistical significance level of 0.05 (Figure 6c). The CET showed opposite responses to vegetation coverage change compared with GE, and LC2015 had a higher CET in the western and southeastern LP, with magnitudes ranging from 2 to 6 mm (Figure 6f). Due to the offset between GE and CET in the process of VR, the TET variation at the annual scale showed a complex and inconsistent pattern in the LP region, with magnitudes ranging from −4 to 2 mm (Figure 6i).   Figure 5 depicts a comparison of the area-averaged annual ET for the two simulations. The simulated GE showed clear differences in their annual time series, and the GE of LC2015 was higher than that of LC2000; however, the CET showed opposite responses to vegetation coverage change, with a larger CET under LC2015 in comparison with LC2000. Consequently, the simulated TET did not show evident differences between LC2000 and LC2015, because the positive effect of VR on CET was offset by the negative effect of VR on GE. The spatial variation in annual average ET associated with VR is illustrated in Figure 6. The whole LP region, except for the small area of the central part and northern edge, showed decreased GE due to vegetation coverage changes from LC2000 to LC2015, particularly in the eastern LP, where GE decreased by more than 15 mm, with a statistical significance level of 0.05 (Figure 6c). The CET showed opposite responses to vegetation coverage change compared with GE, and LC2015 had a higher CET in the western and southeastern LP, with magnitudes ranging from 2 to 6 mm (Figure 6f). Due to the offset between GE and CET in the process of VR, the TET variation at the annual scale showed a complex and inconsistent pattern in the LP region, with magnitudes ranging from −4 to 2 mm (Figure 6i).    . Spatial pattern of the simulated annual average ET (unit: mm) and differences (LC2015 minus LC2000) between LC2000 and LC2015 on the LP. The first column (a,d,g) and second column (b,e,h) show simulated ground evaporation, canopy ET and total ET under LC2000 and LC2015, respectively, and the third column (c,f,i) shows differences in ground evaporation, canopy ET and total ET between LC2015 and LC2000. Areas with significant differences (p < 0.05) based on t-tests are stippled.

Spatiotemporal Change in ET Attributed to VR
Seasonally, the spatial patterns of the responses of GE to vegetation change were similar among the four seasons, showing a consistent decrease in the eastern LP, while the magnitudes varied. June, July, and August (JJA) had the largest decreases, ranging from −3.20 to −2.00 mm, followed by September, October, and November (SON) ranging from −2.00 to −0.80 mm, and December, January, and February (DJF) with approximately 0.40 mm, with a statistical significance level of 0.05 (Figure 7b-d); however, the responses of CET varied among the four seasons. Nearly the whole LP region showed a dominant increase in CET in JJA with vegetation cover change from 2000 to 2015. The largest increase was found in the eastern LP, with a value of approximately 0.8 mm, and the increased magnitudes decreased in the northwest direction (Figure 7f). During SON, mixed negative and positive responses were observed with magnitudes varying from −0.8 mm in the north to 0.8 mm in the southeast (Figure 7g). Notably, the area with increased ET changed in DJF, showing a significant increase in the northeastern LP and a nonsignificant increase in the western LP, where the grassland clearly increased (Figure 7h). The overall change in TET is inconsistent among the four seasons due to the different responses of GE and Figure 6. Spatial pattern of the simulated annual average ET (unit: mm) and differences (LC2015 minus LC2000) between LC2000 and LC2015 on the LP. The first column (a,d,g) and second column (b,e,h) show simulated ground evaporation, canopy ET and total ET under LC2000 and LC2015, respectively, and the third column (c,f,i) shows differences in ground evaporation, canopy ET and total ET between LC2015 and LC2000. Areas with significant differences (p < 0.05) based on t-tests are stippled.
Seasonally, the spatial patterns of the responses of GE to vegetation change were similar among the four seasons, showing a consistent decrease in the eastern LP, while the magnitudes varied. June, July, and August (JJA) had the largest decreases, ranging from −3.20 to −2.00 mm, followed by September, October, and November (SON) ranging from −2.00 to −0.80 mm, and December, January, and February (DJF) with approximately 0.40 mm, with a statistical significance level of 0.05 (Figure 7b-d); however, the responses of CET varied among the four seasons. Nearly the whole LP region showed a dominant increase in CET in JJA with vegetation cover change from 2000 to 2015. The largest increase was found in the eastern LP, with a value of approximately 0.8 mm, and the increased magnitudes decreased in the northwest direction (Figure 7f). During SON, mixed negative and positive responses were observed with magnitudes varying from −0.8 mm in the north to 0.8 mm in the southeast (Figure 7g). Notably, the area with increased ET changed in DJF, showing a significant increase in the northeastern LP and a nonsignificant increase in the western LP, where the grassland clearly increased (Figure 7h). The overall change in TET is inconsistent among the four seasons due to the different responses of GE and CET. Nearly the whole LP exhibited increased TET in JJA (Figure 7j), while the TET decreased in the whole LP during March, April, and May (MAM) and SON with vegetation coverage changes from 2000 to 2015, particularly in the eastern LP (Figure 7i,k). During DJF, only a small area in the southern LP had a decreased TET, and most LP regions exhibited an increased TET (Figure 7l).
CET. Nearly the whole LP exhibited increased TET in JJA (Figure 7j), while the TET decreased in the whole LP during March, April, and May (MAM) and SON with vegetation coverage changes from 2000 to 2015, particularly in the eastern LP (Figure 7i,k). During DJF, only a small area in the southern LP had a decreased TET, and most LP regions exhibited an increased TET (Figure 7l). Figure 7. Spatial distributions of the seasonal differences (LC2015 minus C2000) in ET (unit: mm) between LC2000 and LC2015 on the LP. The first column (a,e,i), second column (b,f,j), third column (c,g,k) and fourth column (d,h,l) show differences in the ground evaporation, canopy ET, and total ET for spring, summer, autumn, and winter, respectively. MAM, JJA, SOM, and DJF represent spring for the months from March to May, summer for the months from June to August, autumn for the months from September to November, and winter for the months from December to February of the following year, respectively.

Quantifying the Contributions of Climate Change and VR to ET Change on the LP
To evaluate the dynamic change in ET associated with VR and climate change over the LP region, the main leading modes of long-term variability were revealed using EOF analysis. We only presented the first EOF modes (EOF1) for analysis, because EOF1 explained more than 25% of the total ET variances, while the second EOF and third EOF modes explained less than 15% and 9% of the variance, respectively. The EOF1 modes of GE for LC2000RCP45 and LC2015RCP45 reflected consistent spatial patterns, showing opposite GE variations between the northern LP and zonal region extending from the central-western area to the east over the LP (Figure 8a,b). Combined with the EOF1 spatial pattern and temporal coefficients, GE exhibited an increasing trend before 2040 in the central-western and central-eastern area, while GE had a decreasing trend most of the time after 2040, in which LC2015RCP45 had a more pronounced change than LC2000RCP45. The EOF1 modes of CET for LC2000RCP45 and LC2015RCP45 showed consistent negative values in their respective spatial patterns, and the time coefficients featured downward trends (Figure 8c,d), indicating that CET increased during 2021-2100, with a more pronounced increase in LC2015RCP45. For TET, the spatial pattern of variability showed a positive value in the whole LP region (Figure 8e,f), which reflected an increasing trend of Figure 7. Spatial distributions of the seasonal differences (LC2015 minus C2000) in ET (unit: mm) between LC2000 and LC2015 on the LP. The first column (a,e,i), second column (b,f,j), third column (c,g,k) and fourth column (d,h,l) show differences in the ground evaporation, canopy ET, and total ET for spring, summer, autumn, and winter, respectively. MAM, JJA, SOM, and DJF represent spring for the months from March to May, summer for the months from June to August, autumn for the months from September to November, and winter for the months from December to February of the following year, respectively.

Quantifying the Contributions of Climate Change and VR to ET Change on the LP
To evaluate the dynamic change in ET associated with VR and climate change over the LP region, the main leading modes of long-term variability were revealed using EOF analysis. We only presented the first EOF modes (EOF1) for analysis, because EOF1 explained more than 25% of the total ET variances, while the second EOF and third EOF modes explained less than 15% and 9% of the variance, respectively. The EOF1 modes of GE for LC2000RCP45 and LC2015RCP45 reflected consistent spatial patterns, showing opposite GE variations between the northern LP and zonal region extending from the central-western area to the east over the LP (Figure 8a,b). Combined with the EOF1 spatial pattern and temporal coefficients, GE exhibited an increasing trend before 2040 in the central-western and central-eastern area, while GE had a decreasing trend most of the time after 2040, in which LC2015RCP45 had a more pronounced change than LC2000RCP45. The EOF1 modes of CET for LC2000RCP45 and LC2015RCP45 showed consistent negative values in their respective spatial patterns, and the time coefficients featured downward trends (Figure 8c,d), indicating that CET increased during 2021-2100, with a more pronounced increase in LC2015RCP45. For TET, the spatial pattern of variability showed a positive value in the whole LP region (Figure 8e,f), which reflected an increasing trend of TET during 2021-2100 by combining with the upward trends exhibited in the time coefficients. TET during 2021-2100 by combining with the upward trends exhibited in the time coefficients. The spatial patterns of ET variabilities under the RCP8.5 scenario were noticeably different from those under the RCP4.5 scenario (Figure 9). The EOF1 mode of GE reflected consistent variations over the northwestern LP, and opposite variations over the rest of the LP in LC2000RCP85 and LC2015RCP85 (Figure 9a,b). The time coefficients demonstrated the temporal behavior of this mode, characterized by an upward trend before 2060 and a downward trend after 2060, indicating that GE had an increasing trend before 2060 and a decreasing trend after 2060 in the southeastern LP, while the northwestern LP showed opposite trends. The spatial patterns of EOF1 for CET were different from those obtained from GE and showed consistent negative values in the whole LP (Figure 9c,d), and their time coefficients had downward trends during 2021-2100, reflecting general increasing trends of CET in LC2000RCP85 and LC2015RCP85. Notably, higher variabilities were detected in the southeastern LP under LC2000RCP85. For TET, both LC2000RCP85 and LC2015RCP85 had similar spatial patterns and time coefficients as CET, suggesting increased TET from 2021 to the end of the 21st century (Figure 9e,f).
Summarizing the results of the EOF analysis, we note that there were no clear differences in spatial patterns obtained from different land cover conditions under the same climate scenario, while the spatial pattern had noticeable differences under different climate scenarios with the same land cover conditions, indicating that the leading EOF modes of ET in LP were clearly climate-dependent. The spatial patterns of ET variabilities under the RCP8.5 scenario were noticeably different from those under the RCP4.5 scenario (Figure 9). The EOF1 mode of GE reflected consistent variations over the northwestern LP, and opposite variations over the rest of the LP in LC2000RCP85 and LC2015RCP85 (Figure 9a,b). The time coefficients demonstrated the temporal behavior of this mode, characterized by an upward trend before 2060 and a downward trend after 2060, indicating that GE had an increasing trend before 2060 and a decreasing trend after 2060 in the southeastern LP, while the northwestern LP showed opposite trends. The spatial patterns of EOF1 for CET were different from those obtained from GE and showed consistent negative values in the whole LP (Figure 9c,d), and their time coefficients had downward trends during 2021-2100, reflecting general increasing trends of CET in LC2000RCP85 and LC2015RCP85. Notably, higher variabilities were detected in the southeastern LP under LC2000RCP85. For TET, both LC2000RCP85 and LC2015RCP85 had similar spatial patterns and time coefficients as CET, suggesting increased TET from 2021 to the end of the 21st century (Figure 9e,f).
Summarizing the results of the EOF analysis, we note that there were no clear differences in spatial patterns obtained from different land cover conditions under the same climate scenario, while the spatial pattern had noticeable differences under different climate scenarios with the same land cover conditions, indicating that the leading EOF modes of ET in LP were clearly climate-dependent.
To confirm the findings from the EOF analysis, we further calculated the contributions of VR and climate change on ET variation through four groups of simulations. Table 1 shows the unique effect of VR on ET in different stages during 2021-2100 on the LP. VR had a negative effect on GE, leading to a decrease in the area-averaged annual GE ranging from −5.9 to −3.6 mm under RCP4.5, while VR had a positive effect on CET, causing an increase in the area-averaged annual CET ranging from 1.3 to 4.7 mm. Under RCP8.5, the negative effect on GE and positive effect on CET intensified during 2021-2080. Since the negative effects on GE were larger than the positive effects on CET under both RCP4.5 and RCP8.5, the area-averaged annual TET showed a weak decrease, with magnitudes ranging from −3.7 to −1.2 mm. Regardless of vegetation cover in 2000 or 2015, climate change from RCP4.5 to RCP8.5 enhanced the area-averaged annual GE during the 2021-2080 period, with increased magnitudes ranging from 3.6 to 10.0 mm ( Table 2). The area-averaged annual CET also had positive responses to climate change except for the period of 2041-2060, with increased magnitudes ranging from 5.0 to 68.4 mm. As a result, the area-averaged annual TET under the two different land cover conditions evidently increased, with magnitudes ranging from 0.1 to 65.3 mm. To confirm the findings from the EOF analysis, we further calculated the contributions of VR and climate change on ET variation through four groups of simulations. Table  1 shows the unique effect of VR on ET in different stages during 2021-2100 on the LP. VR had a negative effect on GE, leading to a decrease in the area-averaged annual GE ranging from −5.9 to −3.6 mm under RCP4.5, while VR had a positive effect on CET, causing an increase in the area-averaged annual CET ranging from 1.3 to 4.7 mm. Under RCP8.5, the negative effect on GE and positive effect on CET intensified during 2021-2080. Since the negative effects on GE were larger than the positive effects on CET under both RCP4.5 and RCP8.5, the area-averaged annual TET showed a weak decrease, with magnitudes ranging from −3.7 to −1.2 mm. Regardless of vegetation cover in 2000 or 2015, climate change from RCP4.5 to RCP8.5 enhanced the area-averaged annual GE during the 2021-2080 period, with increased magnitudes ranging from 3.6 to 10.0 mm ( Table 2). The areaaveraged annual CET also had positive responses to climate change except for the period of 2041-2060, with increased magnitudes ranging from 5.0 to 68.4 mm. As a result, the area-averaged annual TET under the two different land cover conditions evidently increased, with magnitudes ranging from 0.1 to 65.3 mm.   The quantitative contribution was calculated as differences in simulated area-averaged annual ET between vegetation cover conditions of 2000 and 2015 under the same climate scenario (unit: mm). The quantitative contribution was calculated as differences in simulated area-averaged annual ET between scenarios RCP8.5 and RCP4.5 under the same vegetation cover conditions (unit: mm).

Discussion
ET is difficult to monitor at a regional scale due to complex underlying surfaces and large observation costs; thus, we adopted GLEAM-ET products to describe the dynamic changes in ET on the LP, and the spatial characteristics matched with the existing reports well [56,57]. Moreover, a significant increasing trend of TET was detected on the LP during 1980-2018, and the gradients were in line with the results (ranging from 1.34 to 3.45 mm/year) from previous studies at different spatial scales [45,58].
To quantify the influences of VR on ET changes at the regional scale on the LP, sensitivity simulations were conducted using realistic land use data rather than plausible scenarios. We found that the CET was enhanced under vegetation coverage change from 2000 to 2015, which was particularly evident in the southwestern and eastern LP, where the grassland and forestland clearly increased. This finding is well supported by many previous studies, where forested ecosystems are generally observed to have higher canopy transpiration than crop ecosystems, because forests have deeper roots and longer transpiration periods [46,59]. Another widely recognized reason is increased canopy interception, which contributes more to evaporative flux [60,61].
A significant reduction in GE was found over regions with increased LAI of grass and forest. The reduction probably occurred because an increased canopy cover features more shadows and decreases the evaporative demand. There are two main driving factor variations involving decreased evaporative demand. On the one hand, increased canopy cover impedes radiation transmission to the ground surface, leading to a decrease in available heat for evaporation [62][63][64]. In this regard, Raz-Yaseef et al. [65] reported that GE fluxes measured in sun-exposed areas were on average double those in shaded areas. On the other hand, more precipitation was intercepted by increased canopy cover; thus, GE was decreased due to decreased water availability because precipitation is regarded as an important factor driving evaporation in water-limited areas [19,66]. Therefore, a larger reduction in GE relative to increased CET resulted in weak variation in TET at the annual scale due to the offset between the positive contribution and negative contribution of VR to TET. Notably, the TET increase of approximately 1.00 mm existed in most of the LP during the summer period. This widespread increase is mostly the result of improved water availability to satisfy the evaporative demand, because the precipitation on the LP is mainly concentrated in the summer [17,67]; thus, the soil moisture is recharged by precipitation in this period.
Based on four groups of simulations under future climate and land cover conditions, annual ET had increasing trends on the LP with varied magnitude under RCP4.5 and RCP8.5 scenarios. These findings were consistent with a previous study [68], in which the significantly increasing trends are predicted in southeastern LP under RCP4.5 and RCP8.5 scenarios using the 28 general circulation models. We found that climate change exerted a greater impact on TET than vegetation cover change, which matched previous studies well. For example, Peel et al. (2010) studied relationships between forestation and catchment evapotranspiration from over 200 paired catchments around the world and found that non-forested catchments generally have higher actual ET than forested catchments and suggested that the leading factor of ET varied in climate types and depended on the availability of energy and water [19]. Li et al. [24] quantitatively separated the effects of climate change and land cover change on ET in China during the period of 2001-2013 and concluded that the influences of climate change were greater than those of land cover change on ET. Gao et al. [58] used the Budyko Framework to simulate the actual ET in 161 subbasins from 1990 to 2014 on the LP and found that the increase in regional ET was mainly due to the increase in regional precipitation. The LP region has become wetter since the implementation of the Green for Grain Project, with a significant increasing precipitation trend at a rate of 4.46 mm/year during 1998-2014 [69], which is beneficial to the water availability for VR and further increases the ET on the LP.
Some inconsistent results can also be found in regions of the LP, showing that vegetation greening was the dominant driving factor of ET increase in most LP region [45,46]. These conclusions were obtained from a remote sensing-based ET calculation model and partial correlation analysis between ET and some driving factors; however, due to the limited consideration of the biogeophysical mechanism in ET partitioning, the model had uncertainties in estimating ET components [57]. Additionally, the controlled factor was identified based on correlations between independent climatic variables and ET. While the reality is that these variables are not totally independent, their interactions may also play important roles in the ET trends.
Although the impacts of VR and climate change on ET were identified by combining GLEAM-ET dataset with CLM, some uncertainties remain. As reference data for validating the ET simulation of CLM, GLEAM-ET was used with a set of algorithms that used satellite forcing data, and the ET components were calculated using the Priestley-Taylor method, which does not incorporate parameterization of stomatal and aerodynamic resistance. Although GLEAM-ET has been validated against eddy covariance towers worldwide, that datasets were closer to the flux tower measurements at grassland system followed by the forest and cropland dominated regions [70]; however, cropland covers approximately 38% of the land area on the LP [57], where the uncertainties existed in the process of CLM validation, but we did not have a better choice at the regional scale. The understanding of spatiotemporal pattern of ET were obtained using grid-based evaluation statistics, which is more or less subjective, because it precludes spatial analysis; thus, distribution-based field significance are more meaningful [71]. In the sensitivity simulations, the Qian ATM forcing data were used repeatedly to force the CLM, and the vegetation cover conditions of 2000 and 2015 were used for model inputs to assess the ET responses, while the magnitude of vegetation growth might not be large enough to cause significant changes in ET, leading to biased results. Additionally, the relative contributions of VR and climate change were identified through scenario simulation under predicted future climate conditions, while current land surface models commonly underestimate the ratio of plant transpiration to total ET [72], which might contribute to the uncertainty in ET prediction under future climate scenarios.

Conclusions
Along with the implementation of VR, vegetation coverage clearly increases on the LP, inevitably causing TET change. This change involves not only the magnitudes of TET, but also the directions of TET components. In this study, we have quantified the contributions of VR and climate change to TET by analyzing a GLEAM-ET dataset and conducting sensitivity simulations with the CLM. Although GE and CET had larger responses to VR in the southeastern LP, where the LAI of forest and grass clearly increased, the TET had weak change due to the balance between the negative response of GE and positive response of CET. We also found the TET, GE, and CET were especially sensitive to climate change relative to VR, implying that the relative contribution of VR to regional ET may be overestimated in the LP by previous studies, because considering climate factors separately to assess the relative contributions of climate change to ET may mask actual signals. Although there were some uncertainties regarding the data inputs and model parametrization, our results highlight the importance of evaluating ET response to VR under various climate conditions with the interactions among climate factors. Moreover, to ensure the sustainability of vegetation restoration on the LP, reducing evaporative water loss and maintaining soil water availability are required. For this purpose, policymakers should pay more attention to land cover management and take positive steps to slow climate warming.

Data Availability Statement:
The data presented in this study are available on request from the corresponding author for research purposes.