Contributions of Vegetation Greening and Climate Change to Evapotranspiration Trend after Large-Scale Vegetation Restoration on the Loess Plateau, China

: Since the early 2000s, the vegetation cover of the Loess Plateau (LP) has increased signiﬁcantly, which has been fully recorded. However, the effects on relevant eco-hydrological processes are still unclear. Here, we made an investigation on the changes of actual evapotranspiration (ET a ) during 2000–2018 and connected them with vegetation greening and climate change in the LP, based on the remote sensing data with correlation and attribution analysis. Results identiﬁed that the average annual ET a on the LP exhibited an obvious increasing trend with the value of 9.11 mm yr − 1 , and the annual ET a trend was dominated by the changes of ET a in the third quarter (July, August, and September). The future trend of ET a was predicted by the Hurst exponent. Partial correlation analysis indicated that annual ET a variations in 87.8% regions of the LP were controlled by vegetation greening. Multiple regression analysis suggested that the relative contributions of potential evapotranspiration (ET p ), precipitation, and normalized difference vegetation index (NDVI), to the trend of ET a were 5.7%, − 26.3%, and 61.4%, separately. Vegetation greening has a close relationship with the Grain for Green (GFG) project and acts as an essential driver for the long-term development trend of water consumption on the LP. In this research, the potential conﬂicts of water demanding between the natural ecosystem and social-economic system in the LP were highlighted, which were caused by the fast vegetation expansion.


Introduction
Evapotranspiration (ET), as a key link between surface energy balance and water balance, is an important factor connecting the ecological environment and hydrological cycle [1]. Terrestrial ET processes may significantly transform on multiple spatial and temporal scales with the changes of land use and climate [2], which impacts regional vegetation growth, water cycle, and the climate change feedback, particularly in arid and semiarid areas [3]. Additionally, the variation of the regional water cycle can be recognized by the long-term trend of ET [4]. Consequently, a comprehensive understanding of the interrelation between land-use changes, climate change, and the ET trend is critical for ecosystem conservation and water resource management.
The Loess Plateau (LP) is located at the junction of northwest and northern China and belongs to the arid continental monsoon climate zone, where about 85% precipitation is expended by ET [5]. Since the implementation of the Grain for Green (GFG) project, LP has come through considerable land-use changes under large-scale vegetation restoration, whose significant characteristic was the massive transformation from sparse vegetation and farmlands to grassland or shrubbery [6]. During 1999-2010, the GFG project had transformed rain-fed farmland about 16,000 km 2 into grassland and forest [7]. The increasing trend of the normalized difference vegetation index (NDVI) has been significant in the LP since 2000 [8]. As a result, the vegetation coverage of the LP has remarkably increased from 31.6% in 1999 to 59.6% in 2013 [9]. During this period, the change of climate is gradually occurring, with wetter and slightly cooling trends rather than the warming-drying trend before the late 1990s [10]. Therefore, the Plateau ecosystem is significantly influenced by human activities and climate change [11], whose indirect or direct impacts on the hydrologic cycle are still indistinct on multiple spatial-temporal scales.
Land-use changes can affect the ET process through changing surface roughness, interception, and albedo [12]. For the inter-annual trend of ET, the change of climate variable on the long-term scale is also an important driving factor [13]. The early researches on the responses of ET to vegetation change were normally based on the field experiment, which was by means of observing the distinction of site scale ET under different climate conditions or different vegetation cover conditions to confirm the impacts of climate changes and vegetation on ET evolution process [14]. A field experiment is an essential approach in revealing the interrelated mechanism of soil, vegetation, and atmosphere. However, the field scale brings uncertainty of the study universality. In different areas, the results of field experiments are usually inconsistent because of the limitation on specific vegetation types, soil types, and climatic conditions.
Consequently, hydrologic models based on the physical are usually used by many researchers in order to make clear the influences of land use and climate variations in ET, groundwater, soil water, and runoff under different scenes [15,16]. Such physically-based hydrological models are effective for the estimate of regionally hydrologic responses under different land surface conditions and climate. However, there are some noticeable disadvantages in model validation and calibration such as too many parameters, complicated model structures, uncertainty, and time-consuming. Additionally, in other studies, the impacts of underlying surface types and climate changes on the hydrological processes were usually assessed by the empirical models based on the Budyko framework or water balance by some scholars [17][18][19]. Using these empirical models, the uncertainties that existed in the models based on physics are avoided and therefore an effective and simple method for attribution analysis is provided [20,21]. However, due to their characteristic application based on the annual and basin-scale, it is normally not thought that they can obtain the spatial-temporal distributions of ET.
A lot of field experiments have shown that the response and corresponding mechanism of ET to the variation of vegetation conditions are different based on the vegetation composition and type [22]. In addition, the response may vary with various spatial scales [23]. Traditional empirical modeling methods or field experiments cannot fully reflect the response mechanism based on vegetation types and spatial scales. With the development of satellite technology, remote sensing provides a wide range of spatial coverage of land surface features and high spatial and temporal resolution [24]. Therefore, more and more research is utilizing remote sensing based on remote sensing retrievals or ET models to solve the effect of vegetation coverage changes on the ET process. For example, depending on the ET dataset based on satellite estimation, Li et al. [25] researched the impacts of climate variation and land surface pattern on actual ET in the Tarim River, China from 2002 to 2012.
Previous studies have indicated that with the increase of vegetation coverage in the LP, the regional ET has also obviously changed [26]. Additionally, the complexity of understanding and quantifying spatiotemporal patterns and drivers of ET is enhanced by the interaction between climate variation and vegetation greening. However, the impacts of climate variation on vegetation growth were usually ignored due to the fact that many researchers only studied the impacts of vegetation greening and climate change on ET. Consequently, the contribution of human activities to ET variations may be overestimated. In addition, many studies based on remote sensing focus on the response of ET to the vegetation restoration or other influencing variables, but few studies have investigated the process and mechanisms of historical ET, although understanding the evolution mecha-nism of ET is essential for estimating terrestrial ecosystems and land surface-atmosphere interactions and making clear its responses to land cover changes and climate change.
In the study, the precipitation data can be obtained from meteorological stations, and the calculation of vegetation index is relatively simple. However, the potential ET data is often based on the traditional method of evaporating dish measurement and is limited by natural and man-made conditions. It is difficult to obtain large-area ET observation and estimation. Remote sensing technology, because of its large-area, dynamic and rapid access to surface conditions, has had rapid development in recent years, so that large-scale non-uniform land surface ET research has made great progress. For example, Jung et al. [27] provided a data-driven, spatially explicit estimate of global terrestrial ET from 1982 to 2008; the Numerical Terra Dynamic Simulation Group of Montana University researched and developed the global MODIS ET data set, which has been widely used in the study of ET with high spatial and temporal resolution. Liu et al. [28,29] used hydrological and meteorological data to verify the applicability of MODIS ET products in Asia and China, and good results were obtained.
On this basis, the hydrological responses of vegetation greening and climate change in the LP were analyzed by using MOD16 products and NDVI satellite remote sensing data, combined with the meteorological data of 98 national meteorological stations inside and outside the LP from 2000 to 2018. This study is aimed at (1) analyzing the temporal and spatial changes and evolution mechanisms of ET; (2) identifying the leading driver that controls the inter-annual changes of ET; (3) quantitatively distinguishing the contribution of climate and vegetation changes to ET trends.

Study Area
The Loess Plateau is at the junction of northwest and northern China (33 •~4 2 • N, 100 •~1 14 • E), covering an area of 64,000 km 2 . The average altitude of LP is 1398 m. It is one of the four major plateaus in China and also the largest loess plateau in the world (Figure 1), spanning west from the Qinghai Province Sun Moon Mountain, east to Henan Province Taihang Mountain, south Qinling, and north to the Great Wall. The main coverage includes all regions of Shanxi Province, as well as some regions of Qinghai, Shaanxi, Inner Mongolia, Gansu, Ningxia, Henan, and other provinces, a total of 7 provinces (autonomous regions), 46 places (league, state, city), and 282 counties (flags, cities, districts). The area is in the semi-arid and arid transition zone, located on the edge of the warm temperate continental monsoonal climate area, having more prominent mainland and monsoon instability. Being cold and dry in spring and winter, hot and rainy in summer and autumn, it has an average annual precipitation of 150-750 mm and an average annual temperature of 3.6-14.3 • C over the whole Loess Plateau. From northwest to southeast, the vegetation types of the Loess Plateau are desert, desert steppe, typical steppe, and forest-steppe with a gradual increase of leaf area index.

Data
MOD16 global land surface evapotranspiration data set is one of the NASA/EOS projects, its accuracy of using satellite remote sensing data to retrieve global land evapotranspiration has been verified in many aspects. It has been widely used to analyze soil moisture, water, and energy balance in a watershed. Based on the long-term evapotranspira-

Data
MOD16 global land surface evapotranspiration data set is one of the NASA/EOS projects, its accuracy of using satellite remote sensing data to retrieve global land evapotranspiration has been verified in many aspects. It has been widely used to analyze soil moisture, water, and energy balance in a watershed. Based on the long-term evapotranspiration data obtained from the MOD16 product, the effects of meteorology, land cover type changes, and ecosystem disturbances on regional water resource changes can be quantified. In MOD16, a class 4 product of the MODIS series, after geometric correction and radiation correction, each pixel in the data is accompanied by geocoding, reflectivity, and radiance information, which can not only provide characteristic parameters related to evapotranspiration but also have the characteristics of high time resolution (up to 8-day) and low error (error less than 1 pixel), etc.
In the study, the surface evaporative monthly synthesis product set (MOD16A2), including surface evapotranspiration (ET a ) and potential evapotranspiration (ET p ) with 8-day temporal composite at 500 m spatial resolution from 2000 to 2018 in the LP, were obtained from the Level 1 and Atmosphere Archive and Distribution System (LAADS) DAAC of National Aeronautics and Space Administration (NASA). The algorithm used for the MOD16A2 product is based on the Penman-Monteith equation [30]. It is noted that for the spatial areas in the LP such as city, water, glacier, and bare land, data cannot be obtained in the study. The basic information of the data was shown in Table 1. Firstly, the original hierarchical data HDF formats were converted to GeoTiff format by using the MRT projection conversion tool, transforming SIN projection into WGS-1984/geographic longitude and latitude coordinate system and splicing the image. Then, based on the data provided from the NASA website, eliminate invalid values and restore true values. Finally, the actual evapotranspiration (ET a ) and potential evapotranspiration (ET p ) in each year and month of the study area were extracted by vector clipping through ArcGIS software. NDVI datasets were obtained from the Institute of Tibetan Plateau Research, Chinese Academy of Sciences. It is the SPOT VEGETATION NDVI satellite remote sensing dataset based on continuous time series, generated by the maximum value synthesis method since 1998. The distribution and changes of vegetation coverage in all regions of the country at spatial-temporal scale are reflected by this dataset effectively. In the study, NDVI grid files with 10-day temporal composite at 1000 m spatial resolution from 2000 to 2018 inside the Loess Plateau were obtained by splicing and recasting.
Meteorological data, including the temperature and precipitation data of the LP from 2000 to 2018, were obtained from the National Earth System Science Data Center of China. These data were strictly screened to eliminate the abnormal data, and then the univariate linear regression method was used to calculate the interpolation for the meteorological values of individual missing measurements. After a comparison with other interpolation methods such as ordinary kriging (OK) and empirical Bayesian kriging (EBK) [31][32][33], the precipitation data were gridded by inverse distance weighted (IDW) [34] with the spatial resolution of 1 km, which was more efficient and accurate in the cases of dealing a lot of data on the large region.

Methodology
The spatial patterns of the variables were determined by the annual mean of them from 2000 to 2018 in the LP. The slope of the univariate linear regression model [35] determined the trends of ET, NDVI, and climate variables: where k is the slope of the univariate linear regression model; n represents the span of the year; m represents the order of year, and X m represents the variable value for the mth year. The larger value of k is, the more significant the changing trend is.
Then we calculated the Hurst exponent [36] to analyze the self-similarity and longterm dependence of ET a in order to forecast its value trend in the existing circumstances: For {X(t)} (t = 1, 2, · · · , n), its mean of time sequence: Cumulative deviation: Range sequence: Standard deviation sequence: Hurst exponent calculation: where X presents ET a and c is a constant, the Hurst empirical formula can be obtained by taking the logarithm on both sides of Equation (6). The Hurst exponent (H) is obtained by using the Hurst empirical formula based on time series correction, and H represents the trend change with time series. When 0.5 < H < 1, it shows the future change trend is consistent with the past, and the closer to 1, the stronger the consistent tendency is; when H = 0.5, it indicates that the future change trend is random and cannot be predicted; when 0 < H < 0.5, it indicates that the future change trend is opposite to the past, and the closer to 0, the stronger the anti consistency is [37].
To determine the dominant driving forces controlling the inter-annual changes of ET a , the partial correlation analysis between ET a and NDVI, ET p , and P were conducted on the pixel scale. For each pixel, the variable with the highest partial correlation coefficient was judged as the dominant driving factor.
The influence of each factor on ET a change was quantitatively differentiated by multiple regression, and the relative contribution rate of each factor to ET a change was identified. The significant levels of independent variables entering and moving out of multivariate regression models were set to 0.05 and 0.1, respectively. In order to identify the relative contribution rate of each factor to ET a change, we firstly standardize the data by the z-score standardized method. Then, the regression equation of the standardized data sequence is obtained by regression analysis. According to the following formula, the contribution rate of each variable change to the dependent variable change can be calculated [38]: where Y s is the standard value of dependent variables; X 1s ,X 2s ,X 3s . . . is the standard values of independent variables; and a,b,c is the regression coefficient after sequence standardization; η 1 represents the relative contribution rate of X 1 change to Y change; η 2 represents the actual contribution rate of X 1 change to Y change; ∆Y s is the variation of Y s change; ∆X 1 is the variation of X 1s .
The specific technology roadmap is as Figure 2:

Spatial Patterns of ETa, ETₚ, NDVI, and Precipitation
With certain vegetation and climate conditions in the Loess Plateau, precipitatio available energy which is represented by ETₚ are the main driving factors of long annual ETa [39]. As shown in Figure 3a, there was a large spatial variability for annu on the LP, which gradually increased from northwest to southeast areas, ranging fro mm to 936 mm. Compared with the spatial distribution of annual precipitation as s in Figure 3d, the spatial pattern of ETa was consistent with it on the whole. Mul averaged ETa (370 mm) shared about 80% of the corresponding precipitation (466 This result made clear that evaporation was the way for most of the precipitation to to the atmosphere eventually and that precipitation on the multi-year scale mostly lated the changes of annual ETa. For the spatial pattern of annual ETa, the NDVI ( 3c) was more consistent with it (r = 0.66) by comparison with that of the precipitatio it indicated that the annual ETa was mainly regulated by different conditions of v tion. As shown in Figure 3b, the spatial pattern of annual ETp was obviously contr that of ETa (r = −0.48), and the value of annual ETₚ generally decreased from northw southeast areas.

Spatial Patterns of ET a , ET p , NDVI, and Precipitation
With certain vegetation and climate conditions in the Loess Plateau, precipitation and available energy which is represented by ET p are the main driving factors of long-term annual ET a [39]. As shown in Figure 3a, there was a large spatial variability for annual ET a on the LP, which gradually increased from northwest to southeast areas, ranging from 139 mm to 936 mm. Compared with the spatial distribution of annual precipitation as shown in Figure 3d, the spatial pattern of ET a was consistent with it on the whole. Multiyear averaged ET a (370 mm) shared about 80% of the corresponding precipitation (466 mm). This result made clear that evaporation was the way for most of the precipitation to return to the atmosphere eventually and that precipitation on the multi-year scale mostly regulated the changes of annual ET a . For the spatial pattern of annual ET a , the NDVI (Figure 3c) was more consistent with it (r = 0.66) by comparison with that of the precipitation and it indicated that the annual ET a was mainly regulated by different conditions of vegetation. As shown in Figure 3b, the spatial pattern of annual ET p was obviously contrary to that of ET a (r = −0.48), and the value of annual ET p generally decreased from northwest to southeast areas.
to the atmosphere eventually and that precipitation on the multi-year scale mostly regulated the changes of annual ETa. For the spatial pattern of annual ETa, the NDVI ( Figure  3c) was more consistent with it (r = 0.66) by comparison with that of the precipitation and it indicated that the annual ETa was mainly regulated by different conditions of vegetation. As shown in Figure 3b, the spatial pattern of annual ETp was obviously contrary to that of ETa (r = −0.48), and the value of annual ETₚ generally decreased from northwest to southeast areas.

Spatial and Temporal Changes of Climate, NDVI, and ET
For the processes of the regional eco-hydrological cycle, climate conditions and human activities are the important driving factors. Changes in climate variables may have a significant effect on the water consumption of vegetation. Human activities, which can indirectly and even directly impact the growth of natural vegetation by means of agricultural management and ecological restoration, also has a significant influence on the processes of eco-hydrological evolution.

Change in Climate
According to the previous studies, the climate of the Loess Plateau for the long term showed a warmer and dryer trend on the whole [40]. However, in the study, we found that the warming and drying trend had slowed down gradually since 2000 because precipitation showed an increasing trend (2.58 mm yr −1 ) while the rate of temperature rise tended to slow down (0.0128 °C yr −1 ) ( Figure 4). Spatially, there was obvious spatial heterogeneity for the change rate and the trend of precipitation and temperature from 2000 to 2018 in the LP. The annual precipitation change rate was decreasing in the southwest and southeast of the LP and a slowly growing trend in most other regions. In the eastern Shanxi and the central Shaanxi regions, an increasing trend of precipitation was the most obvious and also consistent with the changing trend of vegetation.

Spatial and Temporal Changes of Climate, NDVI, and ET
For the processes of the regional eco-hydrological cycle, climate conditions and human activities are the important driving factors. Changes in climate variables may have a significant effect on the water consumption of vegetation. Human activities, which can indirectly and even directly impact the growth of natural vegetation by means of agricultural management and ecological restoration, also has a significant influence on the processes of eco-hydrological evolution.

Change in Climate
According to the previous studies, the climate of the Loess Plateau for the long term showed a warmer and dryer trend on the whole [40]. However, in the study, we found that the warming and drying trend had slowed down gradually since 2000 because precipitation showed an increasing trend (2.58 mm yr −1 ) while the rate of temperature rise tended to slow down (0.0128 • C yr −1 ) (Figure 4). Spatially, there was obvious spatial heterogeneity for the change rate and the trend of precipitation and temperature from 2000 to 2018 in the LP. The annual precipitation change rate was decreasing in the southwest and southeast of the LP and a slowly growing trend in most other regions. In the eastern Shanxi and the central Shaanxi regions, an increasing trend of precipitation was the most obvious and also consistent with the changing trend of vegetation.
tended to slow down (0.0128 °C yr −1 ) (Figure 4). Spatially, there was obvious spatial heterogeneity for the change rate and the trend of precipitation and temperature from 2000 to 2018 in the LP. The annual precipitation change rate was decreasing in the southwest and southeast of the LP and a slowly growing trend in most other regions. In the eastern Shanxi and the central Shaanxi regions, an increasing trend of precipitation was the most obvious and also consistent with the changing trend of vegetation.

Change in Vegetation
The "Grain for Green" (GFG) project, which includes many measures of vegetation greening such as converting farmland to forest or grassland, planting trees on barren land, and promoting the natural restoration of the mountains, ensures the continuous improvement of vegetation conditions in the LP.
In the study, the NDVI was used to reveal the changes of historical vegetation from 2000 to 2018 in the LP. There was an obvious increasing trend for the NDVI from 2000 to 2018 as shown in Figure 5a (0.0078/year) (p < 0.01). On the quarterly scale (Figure 5b), the NDVI was also found to increase in all quarters. In the second quarter, the increase of NDVI was the most noticeable, being almost seven times that in the fourth quarter. According to the implementation area of the GFG project, the central and southeastern areas of the LP were the significantly increasing spatial extent for the NDVI. However, in the urban areas, the NDVI still significantly decreased because of urbanization.
Water 2021, 13, x FOR PEER REVIEW 9 of 17

Change in Vegetation
The "Grain for Green" (GFG) project, which includes many measures of vegetation greening such as converting farmland to forest or grassland, planting trees on barren land, and promoting the natural restoration of the mountains, ensures the continuous improvement of vegetation conditions in the LP.
In the study, the NDVI was used to reveal the changes of historical vegetation from 2000 to 2018 in the LP. There was an obvious increasing trend for the NDVI from 2000 to 2018 as shown in Figure 5a (0.0078/year) (p < 0.01). On the quarterly scale (Figure 5b), the NDVI was also found to increase in all quarters. In the second quarter, the increase of NDVI was the most noticeable, being almost seven times that in the fourth quarter. According to the implementation area of the GFG project, the central and southeastern areas of the LP were the significantly increasing spatial extent for the NDVI. However, in the urban areas, the NDVI still significantly decreased because of urbanization.

Change in ET
As shown in Figure 6, all the actual ET (ETa) was found to have a significant increasing trend on the annual and quarterly scale. On the annual scale (Figure 6a), the average annual ETa increased with the trend of 9.11 mm yr −1 from 2000 to 2018 in the LP. On the whole, the average annual ETa kept a relatively steady increasing trend and the trend of annual ETa was roughly consistent with that of the NDVI over the same period. Comparing the trend of average ETa and NDVI on the quarterly scale (Figure 6b), the consistency between them was still significantly existent, especially for the second quarter. Additionally, there were other factors influencing the annual ETa changes. For instance, in 2004, the rapidly declining annual ETa was controlled by the significant decline of precipitation, and the slight increase in annual ETa during 2011-2013 may have been the result of the increased precipitation and ETp. On the quarterly scale, it was noticed that the increasing trend of ETa in the third quarter (from July to September) was significantly faster than other quarters in one year. In these months, potential water stress may be caused.

Change in ET
As shown in Figure 6, all the actual ET (ET a ) was found to have a significant increasing trend on the annual and quarterly scale. On the annual scale (Figure 6a), the average annual ET a increased with the trend of 9.11 mm yr −1 from 2000 to 2018 in the LP. On the whole, the average annual ET a kept a relatively steady increasing trend and the trend of annual ET a was roughly consistent with that of the NDVI over the same period. Comparing the trend of average ET a and NDVI on the quarterly scale (Figure 6b), the consistency between them was still significantly existent, especially for the second quarter. Additionally, there were other factors influencing the annual ET a changes. For instance, in 2004, the rapidly declining annual ET a was controlled by the significant decline of precipitation, and the slight increase in annual ET a during 2011-2013 may have been the result of the increased precipitation and ET p . On the quarterly scale, it was noticed that the increasing trend of ET a in the third quarter (from July to September) was significantly faster than other quarters in one year. In these months, potential water stress may be caused. Spatially, during 2000-2018, the NDVI showed the trend of increase in most areas in the LP (93.3%) and presented a significant increasing trend in the 78.4% area of the LP (Figure 7a). The ETa showed the trend of increasing in most areas of the LP (99.4%) and presented a significant increasing trend in the 97.7% area of the LP (Figure 7b). Obviously, the spatial trends of annual ETa were consistent with that of the NDVI, which increased from northwest to southeast. On the contrary, as shown in Figure 7c, ETp in the LP showed a decreasing trend (−4.9 mm yr −1 ) and was found to have an obvious complementary relationship with ETa in the spatial trend pattern by comparing the spatial pattern of ETa trends. Spatially, during 2000-2018, the NDVI showed the trend of increase in most areas in the LP (93.3%) and presented a significant increasing trend in the 78.4% area of the LP (Figure 7a). The ET a showed the trend of increasing in most areas of the LP (99.4%) and presented a significant increasing trend in the 97.7% area of the LP (Figure 7b). Obviously, the spatial trends of annual ET a were consistent with that of the NDVI, which increased from northwest to southeast. On the contrary, as shown in Figure 7c, ET p in the LP showed a decreasing trend (−4.9 mm yr −1 ) and was found to have an obvious complementary relationship with ET a in the spatial trend pattern by comparing the spatial pattern of ET a trends. Spatially, during 2000-2018, the NDVI showed the trend of increase in most areas in the LP (93.3%) and presented a significant increasing trend in the 78.4% area of the LP (Figure 7a). The ETa showed the trend of increasing in most areas of the LP (99.4%) and presented a significant increasing trend in the 97.7% area of the LP (Figure 7b). Obviously, the spatial trends of annual ETa were consistent with that of the NDVI, which increased from northwest to southeast. On the contrary, as shown in Figure 7c, ETp in the LP showed a decreasing trend (−4.9 mm yr −1 ) and was found to have an obvious complementary relationship with ETa in the spatial trend pattern by comparing the spatial pattern of ETa trends.

Future Change Trend of ET a
The mean Hurst exponent of ET a in the LP was 0.46 and the area less than 0.5 accounted for 67.6%, indicating that the reverse characteristics of ET a change in the Loess Plateau were stronger than that in the same direction (Figure 8a). The high-value areas of the Hurst index (H > 0.65) were mainly distributed in the middle-south part of the LP, which exhibits that the ET a trend in these regions will be strongly consistent with the past trend in the future; the low-value areas of the Hurst exponent (H < 0.35) were mainly distributed in the northeast and northwest of LP, showing that the trend of ET a in these regions will be strongly contrary with the past trend in the future; and there were the characteristics of weak co-direction and weak reverse for the trend of ET a in other areas in the future. The trend of ET a and corresponding Hurst exponent were superposed and coupled to reveal the future trend of ET a (Figure 8b). On the whole, ET a with a continuously increasing trend was mainly distributed in the northeast of the Loess Plateau, which included the areas of both very low vegetation cover and very high vegetation cover, accounting for 32.1% area of the LP, while ET a with a trend from increase to decrease was mainly distributed in other areas with relatively abundant vegetation, which accounted for 67.2% area of LP.

Future Change Trend of ETa
The mean Hurst exponent of ETa in the LP was 0.46 and the area less than 0.5 accounted for 67.6%, indicating that the reverse characteristics of ETa change in the Loess Plateau were stronger than that in the same direction (Figure 8a). The high-value areas of the Hurst index (H > 0.65) were mainly distributed in the middle-south part of the LP, which exhibits that the ETa trend in these regions will be strongly consistent with the past trend in the future; the low-value areas of the Hurst exponent (H < 0.35) were mainly distributed in the northeast and northwest of LP, showing that the trend of ETa in these regions will be strongly contrary with the past trend in the future; and there were the characteristics of weak co-direction and weak reverse for the trend of ETa in other areas in the future. The trend of ETa and corresponding Hurst exponent were superposed and coupled to reveal the future trend of ETa (Figure 8b). On the whole, ETa with a continuously increasing trend was mainly distributed in the northeast of the Loess Plateau, which included the areas of both very low vegetation cover and very high vegetation cover, accounting for 32.1% area of the LP, while ETa with a trend from increase to decrease was mainly distributed in other areas with relatively abundant vegetation, which accounted for 67.2% area of LP.

Dominant Driving Forces in the Changes of Inter-Annual ETa
The fluctuations of annual ETa were high (Figure 6a) despite the increasing trend of annual ETa on the whole in the LP. Therefore, it was necessary to make clear the dominant forces of ETa for the causation of high fluctuations in annual ETa. The partial correlation analyses on the pixel scale between ETa and its driving forces (NDVI, precipitation, and ETp) in the LP were used to reveal the dominant driving forces that controlled the interannual change of ETa from 2000 to 2018. As shown in Figure 9a, the spatial consistency for partial correlation coefficients(r) between ETa and NDVI was high over the whole study area. The positive correlation areas were distributed on most areas of the LP (97.3%) and the average value of r was 0.65 on the pixel scale. High r value areas were mainly in the eastern and southern areas of the LP which were abundant in precipitation; low r value areas were mainly in the northwestern part of the LP where the vegetation was relatively sparse. For the ETa and precipitation, their correlations were positive on the whole ( Figure  9c). However, the average value of r on the pixel scale was only 0.12, mainly located in the mid-north part of the LP. The correlation between ETa and P was very weak and even negative in some areas such as the southwest of the LP which is rich in farmland where water pressure was significantly alleviated by artificial irrigation. The partial correlation

Dominant Driving Forces in the Changes of Inter-Annual ET a
The fluctuations of annual ET a were high (Figure 6a) despite the increasing trend of annual ET a on the whole in the LP. Therefore, it was necessary to make clear the dominant forces of ET a for the causation of high fluctuations in annual ET a . The partial correlation analyses on the pixel scale between ET a and its driving forces (NDVI, precipitation, and ET p ) in the LP were used to reveal the dominant driving forces that controlled the interannual change of ET a from 2000 to 2018. As shown in Figure 9a, the spatial consistency for partial correlation coefficients(r) between ET a and NDVI was high over the whole study area. The positive correlation areas were distributed on most areas of the LP (97.3%) and the average value of r was 0.65 on the pixel scale. High r value areas were mainly in the eastern and southern areas of the LP which were abundant in precipitation; low r value areas were mainly in the northwestern part of the LP where the vegetation was relatively sparse. For the ET a and precipitation, their correlations were positive on the whole (Figure 9c). However, the average value of r on the pixel scale was only 0.12, mainly located in the mid-north part of the LP. The correlation between ET a and P was very weak and even negative in some areas such as the southwest of the LP which is rich in farmland where water pressure was significantly alleviated by artificial irrigation. The partial correlation coefficient values between ET a and ET p presented conspicuous heterogeneity on the spatial pattern (Figure 9b), where positive correlations were mainly in the southeast part of the LP and negative correlations were mainly in the northwest part of the LP.
According to the spatial patterns of partial correlation coefficients, the dominant driving factor of ETa changes for each grid cell is shown in Figure 9d. The inter-annual changes of ETa were controlled by the NDVI in most areas of the LP (87.9%), while ETp controlled that only in 0.1% of LP areas. For the precipitation, it played a leading role in 6.7% of the whole study area which was mainly distributed in the northwestern part of the LP.

Contributions of Vegetation Greening and Climate Variation to the Trend of ETa
As shown in Figure 10, the spatial heterogeneities in the contribution of the NDVI, ETp, and precipitation to ETa were significant. In general, the ETa trend was affected positively by the vegetation greening (Figure 10a). The contribution of the NDVI to the ETa trend was relatively high in the southeastern areas of the LP, which was rich in farmlands and woodlands with contribution values of more than 10 mm yr −2 and even more than 23 mm yr −2 in some regions. The relatively low contribution values of the NDVI to the ETa trend were mainly distributed in the northwestern part of the LP, where the trend of greening was weak and woodlands were dominant with the contribution values of about 3 mm yr −2 . For the precipitation (Figure 10b), climate variation also had an effect on the According to the spatial patterns of partial correlation coefficients, the dominant driving factor of ET a changes for each grid cell is shown in Figure 9d. The inter-annual changes of ET a were controlled by the NDVI in most areas of the LP (87.9%), while ET p controlled that only in 0.1% of LP areas. For the precipitation, it played a leading role in 6.7% of the whole study area which was mainly distributed in the northwestern part of the LP.

Contributions of Vegetation Greening and Climate Variation to the Trend of ET a
As shown in Figure 10, the spatial heterogeneities in the contribution of the NDVI, ET p, and precipitation to ET a were significant. In general, the ET a trend was affected positively by the vegetation greening (Figure 10a). The contribution of the NDVI to the ET a trend was relatively high in the southeastern areas of the LP, which was rich in farmlands and woodlands with contribution values of more than 10 mm yr −2 and even more than 23 mm yr −2 in some regions. The relatively low contribution values of the NDVI to the ET a trend were mainly distributed in the northwestern part of the LP, where the trend of greening was weak and woodlands were dominant with the contribution values of about 3 mm yr −2 . For the precipitation (Figure 10b), climate variation also had an effect on the increase of ET a . But its contributions were generally low generally, of only 0.52 mm yr −2 . In contrast, ET p change had a negative impact on the trend of ET a in most areas of the LP with the average contribution values of −2.4 mm yr −2 .
Compared with the value of spatially averaged ETa trend (9.11 mm yr −2 ), it was clear to find that the relative contributions of ETp, P, and NDVI to ETa change were −26.3% (−2.4 mm yr −2 ), 5.7% (0.52 mm yr −2 ), and 61.4% (5.59 mm yr −2 ) respectively over the whole Loess Plateau. The fact indicated that the vegetation greening represented by the NDVI was the dominant driving force to the long-term trends of ETa during 2000-2018 on the Loess Plateau (Figure 10c).

Discussion
It is generally recognized that ETp, vegetation condition, and soil water stress (mainly dominated by precipitation) are three main factors that control ETa [41]. In the study, we assessed the spatial patterns of ET, NDVI, and precipitation, and analyzed their temporalspatial trends, respectively. From 2000 to 2018, ETa in the Loess Plateau showed a significant increasing trend of 9.11 mm yr −1 with the pattern of low value in the northwest and high value in the southeast. At the same time, NDVI, P and ETp respectively exhibited the trend of 0.0078 yr −1 , 2.58 mm yr −1 , and −4.9 mm yr −1 , and the trend of inter-annual ETₚ was Compared with the value of spatially averaged ET a trend (9.11 mm yr −2 ), it was clear to find that the relative contributions of ET p , P, and NDVI to ET a change were −26.3% (−2.4 mm yr −2 ), 5.7% (0.52 mm yr −2 ), and 61.4% (5.59 mm yr −2 ) respectively over the whole Loess Plateau. The fact indicated that the vegetation greening represented by the NDVI was the dominant driving force to the long-term trends of ET a during 2000-2018 on the Loess Plateau (Figure 10c).

Discussion
It is generally recognized that ET p , vegetation condition, and soil water stress (mainly dominated by precipitation) are three main factors that control ET a [41]. In the study, we assessed the spatial patterns of ET, NDVI, and precipitation, and analyzed their temporalspatial trends, respectively. From 2000 to 2018, ET a in the Loess Plateau showed a significant increasing trend of 9.11 mm yr −1 with the pattern of low value in the northwest and high value in the southeast. At the same time, NDVI, P and ET p respectively exhibited the trend of 0.0078 yr −1 , 2.58 mm yr −1 , and −4.9 mm yr −1 , and the trend of inter-annual ET p was found to have an obvious complementary relationship with ET a in the LP which was coincident with the theory of evapotranspiration complementarity [42]. Considering the trend of NDVI and P in the LP, they are all consistent with that of ET a on the whole, and the similar results can also be found in some other studies [26,43]. There is no doubt that the improvement of climatic conditions (an increasing trend in precipitation and tending to slow down for the rate of temperature rise) on the LP in the past two decades is a matter worthy of celebration. However, according to the calculation of the Hurst exponent, under the existing change trends of climate and vegetation greening, the ET a in a third of the LP will continue to increase in the future and it may aggravate the problem of regional water resources shortage. More attention should be paid to the significant increasing tendency of ET a in the LP.
Then, we further researched the causes for the temporal and spatial variation of ET a . From 2000 to 2018, the contribution of vegetation greening to ET a increase reached 61.4%, and vegetation greening was identified as the main cause of the ET a trend. At the same time, there were −26.3% and 5.7% of the ET a trend being contributed by ET p and climate change respectively. These results showed that the reduction of surface albedo and the increase of net radiation absorbed by the surface on the LP with the increase of green leaf area and much deep soil water absorbed by plants during growth can bring a positive effect to ET a . At the same time, increasing precipitation and temperature also promoted the increase of ET a . However, our prediction showed that the spatial patterns of annual ET a trends on the LP will change a lot and its tendency of decrease in the future will be stronger than that of increase. It indicated that the increase of the amount of actual moisture infiltration into the soil after large-scale vegetation restoration will be higher than the increase of vegetation leaf evaporation caused by the increase of vegetation cover in most of the LP at some point in the future. According to the comparison of spatial patterns between the NDVI and the future trend of ET a in the LP, it was found that there are probably three turning points for the trend of ET a with the increase of vegetation coverage.
We roughly divide the current vegetation coverage on the LP into three categories: low, medium, and high, and then divide it into four stages to illustrate the corresponding tendency of ET a in the future. When the vegetation coverage is lower than the low level, the future ET a will decrease; when it is lower than medium grade, the ET a will continue to increase; when it is lower than high grade, the future ET a will exhibit a decreasing trend; when it is higher than high grade, the future ET a will exhibit an increasing trend again. Obviously, in order to play a good role in improving the ecological environment and avoid the conflict between supply and demand of regional water resources, it is the most reasonable afforestation activity to cause the third stage of vegetation coverage which is higher than medium grade and lower than high grade. It exhibits the importance of reasonable scales of vegetation construction. These phenomenons showed the possible existence of the turning point for the ET a trend, and the quantitative analysis based on this will be carried out in the future study.
Large-scale vegetation restoration can improve the diversity of the ecosystem and the capacity of water-holding, but the increase of water requirement caused by early vegetation seedling may lead to the change of regional water budget and cause the conflict of potential water demand between the socio-economic system and natural ecosystem. When the actual ET was generally affected by the change of vegetation pattern, the distributive proportion of runoff and ET a in precipitation also changed with it. Multi-year averaged ET a accounted for nearly 80% of the corresponding precipitation during 2000-2018 on the LP, and this proportion may continue to increase for a long time according to the significant increasing trend of ET a in recent years. As a result, the situation of runoff was worrying. For example, the study of Li et al. [44] identified that when the vegetation coverage had an increase of 1% since the 1970s, the runoff coefficient would respectively decrease by 0.0051 and 0.0009 in the semi-dry and semi-humid regions of the LP. For soil water, its responses to afforestation are relatively complex. On the one hand, the increase of vegetation cover will lead to the increase of the soil water storage capacity and surface roughness, and the infiltration of rainstorm events will also be increased accordingly. On the other hand, the increase of vegetation cover can also mean the increase of rainfall interception, which leads to the decrease of net precipitation. In summer, the rainfall in the LP mainly occurs in the form of a rainstorm, so it is speculated that the soil moisture in the corresponding period will increase. However, according to the demonstrations of many studies, the soil desiccation on the LP is widespread [45,46], indicating that the additional water consumption of higher vegetation coverage has exceeded the amounts of increased infiltration. In order to maintain both the ecological sustainability and socio-economic benefits, in future reforestation, the threshold of vegetation coverage in the LP should be paid more attention, because the soil water availability actually controlled the degree of vegetation restoration in arid and semiarid regions [47]. Therefore, more reasonable types, densities, and scale of vegetation construction based on the balance of water resources consumption and supply and more adequate utilization of rainfall-runoff in vegetation restoration should be identified to guide the future implementation of ecological restoration strategies.

Conclusions
This study revealed the spatial-temporal pattern of ET, vegetation greening, and precipitation in the LP from 2000 to 2018 respectively, and identified their spatial-temporal variation characteristics using the univariate linear regression model. Then according to the Hurst exponent, the future trend of ET a was predicted. Through partial correlation analysis, the dominant driving factor was ascertained, which controlled the inter-annual change of ET a on each pixel. Additionally, contributions of vegetation greening and climate variation to ET a trend were quantitatively distinguished on the basis of multiple regression analysis.
During 2000-2018, the average annual ET a on the LP exhibited an obvious increasing trend with the value of 9.11 mm yr −1 , and the annual ET a trend was dominated by the changes of the third quarterly ET a . It is noticed that the reverse characteristics of ET a change in the Loess Plateau are stronger than that in the same direction, but ET a in a third of the LP may continue to increase in the future. Partial correlation analysis indicated that the inter-annual changes of ET a were controlled by the NDVI in most areas of the LP (87.9%), while ET p controlled that only in 0.1% of areas in the LP. For the precipitation, it controlled that of ET a in 6.7% of the whole study area, which was mainly in the northwestern part of the LP with relatively low precipitation. Over the whole LP, the relative contributions of ET p , precipitation, and the NDVI to ET a trends were −26.3% (−2.4 mm yr −2 ), 5.7% (0.52 mm yr −2 ), and 61.4% (5.59 mm yr −2 ), respectively. This fact indicated that the vegetation greening was the dominant driving force to the long-term trend of ET a on the LP. In order to keep the sustainable eco-hydrological environment in the LP, we suggest adopting more reasonable types, densities, and scale of vegetation construction based on the balance of water resources consumption and supply, and making more adequate utilization of rainfall-runoff in the future implementation of ecological restoration strategies.