Evaluating the Impacts of Climate Change and Vegetation Restoration on the Hydrological Cycle over the Loess Plateau, China

: In recent decades, both observation and simulation data have demonstrated an obvious decrease in runo ﬀ and soil moisture, with increasing evapotranspiration, over the Loess Plateau. In this study, we employed a Variable Inﬁltration Capacity model coupled with scenario simulation to explore the impact of change in climate and land cover on four hydrological variables (HVs) over the Loess Plateau, i.e., evapotranspiration (ET), runo ﬀ (Runo ﬀ ), shallow soil moisture (SM1), and deep soil moisture (SM2). Results showed precipitation, rather than temperature, had the closest relationship with the four HVs, with r ranging from 0.76 to 0.97 ( p < 0.01), and this was therefore presumed to be the dominant climate-based driving factor in the variation of hydrological regimes. Vegetation conversion, from cropland and grassland to woodland, signiﬁcantly reduced runo ﬀ and increased soil moisture consumption, to sustain an increased ET, and, assuming that the reduction of SM2 is entirely evaporated, we can attribute 71.28% ± 18.64%, 65.89% ± 24.14% of the ET increase to the water loss of SM2 in the two conversion modes, respectively. The variation in HVs, induced by land cover change, were higher than the expected climate change with respect to SM1, while di ﬀ erent factors were selected to determine HVs variation in six catchments, due to di ﬀ erences in the mode and intensity of vegetation conversion, and the degree of climate change. Our ﬁndings are critical for understanding and quantifying the impact of climate change and vegetation conversions, and provide a further basis for the design of water resources and land-use management strategies with respect to climate change, especially in the water-limited Loess Plateau.


Introduction
Climate and land cover are the two main factors regulating global and regional hydrological cycles [1,2]. Climate, which is primarily represented by the spatial-temporal distribution of net radiance, precipitation, temperature, and wind speed, controls the water flux and state of land [3]. Land cover, especially various vegetation classes, plays a key role in influencing the water and energy balance at the land surface, affecting a series of hydrologic processes, including precipitation, interception, runoff, transpiration, and evaporation [4]. Many studies have found that the hydrological cycle in many regions worldwide changes significantly in the context of climate change and land cover change [5,6]. Generally, changes in climate, especially precipitation, directly influence the amount of water input and output in a region, while change in land cover affects the timing and magnitude of evaporative

Study Area
The Loess Plateau (Figure 1), located in the upper and middle reaches of the Yellow River (100 • 54 -114 • 33 E and 33 • 43 -41 • 16 N), covers approximately 640,000 km 2 , and is characterized by its thick loess soils with a depth of 30 m to 80 m [40]. It is a typical semiarid and arid climate, with mean annual precipitation varying from 150 mm·yr −1 in the northwest to 800 mm·yr −1 in the southeast, and a mean value of about 450 mm·yr −1 . It is a typical water-limited landscape, with evaporation accounting for 85% of the precipitation, and this percentage increases along the precipitation gradient [26]. About 55%-78% of the annual precipitation falls during the main vegetation growing season between June and September, much of it in high-intensity rainstorms [41]. Precipitation is the primary water source for local vegetation, as groundwater depth in the Loess Plateau is typically deeper than 20 m below the land surface [42]. During the period 1961-2010, the trends of annual mean temperature and precipitation in the Loess Plateau were 0.29 ± 0.1 • C/decade and −11.03 ± 13 mm/decade, indicating an obvious warming and drying climate [17].
The Loess Plateau of China is well known for severe soil erosion and heavy sediment load because of its low vegetation cover, frequent high-rainfall storms in summer months, highly erodible loessal soil, steep topography, and long history of intensive cultivation and unsustainable land use [8]. Besides the Loess Plateau, another six catchments, with different levels of climate change and land cover change, were selected to analyze the effects of CC and LUCC on hydrologic processes in detail. The Loess Plateau of China is well known for severe soil erosion and heavy sediment load because of its low vegetation cover, frequent high-rainfall storms in summer months, highly erodible loessal soil, steep topography, and long history of intensive cultivation and unsustainable land use [8]. Besides the Loess Plateau, another six catchments, with different levels of climate change and land cover change, were selected to analyze the effects of CC and LUCC on hydrologic processes in detail.

Model Description
In this study, the Variable Infiltration Capacity (VIC) [43,44] model is employed, to simulate hydrologic processes in the Loess Plateau. VIC model is a physically-based macroscale hydrologic model, which has been applied in many basins all over the world [45][46][47][48][49]. As a semi-distributed hydrological model, it is able to resolve the complete water and energy balance of multiple land cover classes within each computational grid cell using a mosaic-type vegetation scheme. In this study, the Loess Plateau is divided into 9606 grid cells with a spatial resolution of 0.083° × 0.083°. The representation is relatively detailed, including multiple evaporation and transpiration pathways, snow accumulation and melt, and soil freeze/thaw in three soil layers [44], and 10 thermal nodes, with a thermal damping depth of 4 m, with the finite difference solution technique [50]. The variable infiltration curve is used to consider the spatial heterogeneity of the runoff generation [51]. ET and runoff are simulated for each type of land cover separately, and a weighted average, using the fractions of vegetation types, is calculated, to obtain grid-cell average ET and runoff. Baseflow, which only happens in the third layer, is described by ARNO model [52]. A conceptual surface runoff model with Philip infiltration formulation is used to generate runoff from the first and second layers [53,54]. A separate routing model (RVIC) is coupled with VIC model to simulate streamflow [55,56], where the runoff generated in each grid cell is routed to selected points through the channel network [57].

Model Description
In this study, the Variable Infiltration Capacity (VIC) [43,44] model is employed, to simulate hydrologic processes in the Loess Plateau. VIC model is a physically-based macroscale hydrologic model, which has been applied in many basins all over the world [45][46][47][48][49]. As a semi-distributed hydrological model, it is able to resolve the complete water and energy balance of multiple land cover classes within each computational grid cell using a mosaic-type vegetation scheme. In this study, the Loess Plateau is divided into 9606 grid cells with a spatial resolution of 0.083 • × 0.083 • . The representation is relatively detailed, including multiple evaporation and transpiration pathways, snow accumulation and melt, and soil freeze/thaw in three soil layers [44], and 10 thermal nodes, with a thermal damping depth of 4 m, with the finite difference solution technique [50]. The variable infiltration curve is used to consider the spatial heterogeneity of the runoff generation [51]. ET and runoff are simulated for each type of land cover separately, and a weighted average, using the fractions of vegetation types, is calculated, to obtain grid-cell average ET and runoff. Baseflow, which only happens in the third layer, is described by ARNO model [52]. A conceptual surface runoff model with Philip infiltration formulation is used to generate runoff from the first and second layers [53,54]. A separate routing model (RVIC) is coupled with VIC model to simulate streamflow [55,56], where the runoff generated in each grid cell is routed to selected points through the channel network [57].

Model Calibration and Validation
The six parameters of the VIC model are difficult to determine and thus require calibration, namely the variable infiltration curve parameter (b), the maximum velocity of base flow (Ds max ), the fraction of Ds max where non-linear base flow begins (Ds), the fraction of maximum soil moisture where non-linear baseflow occurs (Ws), and the thickness of two soil moisture layers (d i , i = 2, 3).
The VIC model was calibrated by adjusting the six parameters individually to optimize the objective function, which compared the simulated and observed streamflow [47]. We set 1960 to 1964 as the warm-up period, 1965 to 1974 as the calibration period, and 1975 to 1984 as the validation period for each catchment. The Nash-Sutcliffe coefficient (NSE) and relative bias (BIAS) were employed as the optimal functions for the model calibration and validation: where Q i,s and Q i,o are the simulated and observed streamflow series (m 3 /s), and Q s and Q o are the respective mean values.

Input Data
The data for driving VIC model contained two parts, i.e., forcing data and parameter database. Daily meteorological data of 112 stations within or around the Loess Plateau, obtained from China Meteorological Administration (CMA), were used in this study, which consisted of four variables: precipitation, maximum and minimum temperature, and wind speed for the period 1960-2015. These stations were selected based on data quality controls, to exclude stations with abnormal values. The meteorological data was interpolated to each 0.083 • × 0.083 • grid by a linear interpolation method, in which an inverse squared distance between the stations and the center of the target grid was used as a weight to transfer the meteorological variables from the station to the target grid [47]. Notably, to reflect the decrease in temperature with increasing elevation, the temperature data of each station were adjusted to the same elevation as the target grid, using a lapse rate of −6.5 • C·km −1 with respect to the elevation difference between the station and the target grid [58]. In addition, a trend-removed meteorological forcing data was generated, based on the interpolated dataset, following the method mentioned in previous studies [18,59].
The vegetation parameters were derived from land use/cover data (100 m), based on the Landsat TM data, acquired from the Resource and Environment Data Cloud Platform (RESDC), Chinese Academy of Sciences (http://www.resdc.cn) [60]. Two land use and land cover (LULC) maps for the two representative years of 1990 and 2010 were employed, to characterize the vegetation condition before and after 1999. Additionally, the primary characteristic of land cover that affected hydrologic fluxes in VIC model was LAI. To make the LAI more reasonable, to reflect the land cover change of each grid, two LAI datasets of period 1984-1999 and 2000-2015 were produced by GIMMS LAI, obtained from the Advanced Very High-Resolution Radiometer (AVHRR) satellites. Both LAI datasets contained 12-month LAI of each VIC grid and restored vegetation parameter files. Based on the two LULC maps and the LAI datasets, two sets of vegetation parameters (namely vege1990, vege2010) for different periods were produced for VIC simulation. The soil parameters were derived from Soil Map Based Harmonized World Soil Database (HWSD v1.1). For each 0.083 • × 0.083 • grid cell, the dominant soil type was used. The soil parameters that were not available from the soil data set were calculated according to methods by Saxton and Rawls [61]. Elevation data obtained from the Shuttle Radar Topography Mission (SRTM) high-resolution digital elevation dataset was used to extract each catchment and delineate river networks.

Other Data for Evaluation
Additionally, some other data were employed for model calibration, validation and further evaluations of model performance. The streamflow data from the period 1960-2015 of 16 hydrologic stations were obtained from the Annual Hydrological Report for P.R. China. Soil moisture observations for 28 stations across the Loess Plateau were obtained from the China Meteorology Administration (http://data.cma.cn/data/). The surface soil saturation degree (as a percentage) was presented as the 10-days average for 10, 20, 50, 70, and 100 cm depths from 1992 to 2007. Mean value of soil moisture at 20, 50, 70, and 100 cm depths was calculated as the average soil moisture of 10-100 cm depth. Close attention was paid, to ensure the observed soil moisture was expressed in gravimetric water content (kg/kg), which differed from the flux data of VIC with volumetric unit. Therefore, we focused on temporal variations rather than their absolute values [62] and employed the anomaly of soil moisture, calculated by the ratio of difference between observed and mean soil moisture to mean soil moisture, in the following evaluation. A three-layers (0-5 cm, 5-20 cm, 20-100 cm) of the soil moisture data with 0.25 • × 0.25 • resolution during 2002-2011 was provided by Yang [63], estimated through a Dual-Pass Microwave Land Data Assimilation System (ITPLDAS). For ET evaluation, we used three ET products derived from different methods, i.e., MOD16A2 [64], Zhang Ke ET datasets [65] and Jung global land ET product (http://www.bgc-jena.mpg.de/geodb/projects/Data.php) [66]. MOD16A2 and Zhang Ke ET datasets were estimated by diagnostic models, and used remote sensed vegetation information as an important input, using the Penman-Monteith (PM) algorithm (Zhang et al.; Mu et al.). Another ET dataset was Jung global land ET product, which was produced by upscaling FLUXNET observations from 198 global flux monitoring towers to the global scale, using a data-adaptive machine learning technique (i.e., model tree ensembles, MTE).

Scenario Simulation Design and Analysis
In this study, we focused on four key hydrological variables (HVs), i.e., evapotranspiration (ET), runoff (Runoff), shallow soil moisture (SM1, 0-50 cm), and deep soil moisture (SM2, 50-200 cm). Four scenarios, with different vegetation parameters and climate conditions, were adopted, to investigate the impact of CC and LUCC on hydrological processes (Table 1). The scenario0 (S0), using two vegetation parameters for two periods before and after 1999, and the real climate conditions of the period 1984-2015, simulated the real hydrological cycle, and was further used to detect the variations in hydrological fluxes and states during the period 1984-2015, as shown in Sections 3.1 and 3.2. The scenario1 (S1) also used two vegetation parameters as input, as well as detrended meteorological variables. Scenario2 (S2) and scenario3 (S3) both used one vegetation parameter (veg1990) for the whole period, without respect to LUCC, but used real and detrended meteorological variables, respectively. The simulation in S2 was a continuous time series dataset and correlation analysis, employed to detect the relationship between climate and hydrologic processes, and the variations in forcing data and simulation in S0 and S1 were used to explore the driving effects of CC on HVs, as shown in Section 3.3.1. The variations of simulation in S0 and S2 for the period 2000-2015 were induced by LUCC, as discussed in Section 3.3.2. In Section 3.4, six catchments, with different levels of CC and LUCC, were selected, to analyze the impacts of CC and LUCC on the hydrologic cycle in detail. Compared with S0, meteorological variables in S1 were detrended to identify the effect of CC. The vegetation parameter in S2 was static and the same as that of the period 1984-1999 (Veg1990) to identify the effect of LUCC, while S3 employed both static vegetation parameters and detrended meteorological variables, which simulated the hydrologic processes without CC and LUCC, to identify the co-effects of CC and LUCC. The effects were separated based on the scenario simulations, calculated as follows: where V s0 , V s1 , V s2 and V s3 are the mean value of four HVs under scenarios S0, S1, S2, and S3, respectively. As set in the vegetation parameters and meteorological forcing data, ∆V 0 C , ∆V 0 LC and ∆V 0 Co were the variations in HVs between S1, S2, S3, and S0 for the period 2000-2015, treated as contributions of CC, LUCC, and co-effects of CC and LUCC, respectively. A positive value represents an increase of each HV, while a negative value represents a decline.

Model Calibration and Validation
In this study, eight hydrological stations, with long-term monthly streamflow observations and minimal human disturbance, were selected to calibrate the model ( Figure 2). The NSE of simulated and observed streamflow was above 0.70 at six out of eight stations in the calibration and validation period (Table A1 in Appendix A). The BIAS between simulated and observed streamflow ranged from −15% to 16%, excepting Jiaokouhe station, in the validation period. The results indicated that the parameters of each catchment were reasonable, and the calibrated VIC model was generally acceptable at reproducing the observed monthly streamflow. Therefore, the model parameters of the uncalibrated catchments were set to be identical to the nearest calibrated catchments in the same major river basin for the simulation of the Loess Plateau [57]. Another eight hydrological stations, with larger catchment areas, covering the catchments used for model calibration, were chosen, to validate the parameter implementation (Table A2). The NSE were at or above 0.63 at seven of the eight stations, excepting Baijiachuan station, and the BIAS was in the range −18.94%-26.81%.

Comparisons of Simulated ET and SM with Other Datasets
Evapotranspiration and soil moisture were another two key hydrological variables investigated in this study. Evaluations of ET and SM were essential for the validation of model performance and further studies. Three ET products, with the common data period of 2001-2011, were employed to evaluate the simulate ET in this study (ETVIC). Results showed that ETVIC had a good correlation with ETMODIS (r = 0.72, p < 0.05) and ETZK (r = 0.59, p < 0.05), which indicated that ETVIC was consistent with ETMODIS and ETZK in capturing the inter-annual variation of ET. Compared with the ETVIC, ETZK and ETMTE underestimated ET by 18.07% and 10.89%, and ETMODIS overestimated ET by 8.54%. Though it was well accepted that ETMODIS underestimated actual ET, part of the desert with low ET in the northwest Loess Plateau was treated, as no data may account for the overestimation. The underestimation of ETZK may be because the model did not consider the evaporation from canopy interception. Generally, the spatial pattern of multi-year average ETVIC was similar to ETMODIS and ETZK.
Additionally, a soil moisture dataset from Yang (SMYK), and soil moisture observation from 28 stations (SMobs) were employed to evaluate the simulation by VIC (SMVIC). The comparison was conducted between SMYK and SMVIC within a depth of 0-100 cm, after conversion of units from m 3 /m 3 to mm. The multi-year average soil moisture of the two datasets was comparable, with a bias of -14.12%. The correlation coefficient was 0.38 (p > 0.05), which indicated an obvious difference in soil moisture simulation between VIC and ITPLDAS. The Pearson correlation coefficient (r) between simulated and observed soil moisture of 28 stations ranged from 0.44 to 0.81 for 0-10 cm, and from 0.40 to 0.72 for 10-100 cm ( Figure 3). In general, the correlation coefficients of 0-10 cm were higher than that of 10-100 cm, and showed a significant increase with increased precipitation, while a similar trend was not detected for 10-100 cm. Three stations with a long-time span of observation data and remarkable difference in precipitation were selected to represent the typical arid, semiarid and semihumid region, as shown in Figure 4. The variability of SMVIC was generally lower than that of SMobs, which may be partly attributed to the difference in time-space scale. SMVIC represented the average

Comparisons of Simulated ET and SM with Other Datasets
Evapotranspiration and soil moisture were another two key hydrological variables investigated in this study. Evaluations of ET and SM were essential for the validation of model performance and further studies. Three ET products, with the common data period of 2001-2011, were employed to evaluate the simulate ET in this study (ET VIC ). Results showed that ET VIC had a good correlation with ET MODIS (r = 0.72, p < 0.05) and ET ZK (r = 0.59, p < 0.05), which indicated that ET VIC was consistent with ET MODIS and ET ZK in capturing the inter-annual variation of ET. Compared with the ET VIC , ET ZK and ET MTE underestimated ET by 18.07% and 10.89%, and ET MODIS overestimated ET by 8.54%. Though it was well accepted that ET MODIS underestimated actual ET, part of the desert with low ET in the northwest Loess Plateau was treated, as no data may account for the overestimation. The underestimation of ET ZK may be because the model did not consider the evaporation from canopy interception. Generally, the spatial pattern of multi-year average ET VIC was similar to ET MODIS and ET ZK .
Additionally, a soil moisture dataset from Yang (SM YK ), and soil moisture observation from 28 stations (SM obs ) were employed to evaluate the simulation by VIC (SM VIC ). The comparison was conducted between SM YK and SM VIC within a depth of 0-100 cm, after conversion of units from m 3 /m 3 to mm. The multi-year average soil moisture of the two datasets was comparable, with a bias of -14.12%. The correlation coefficient was 0.38 (p > 0.05), which indicated an obvious difference in soil moisture simulation between VIC and ITPLDAS. The Pearson correlation coefficient (r) between simulated and observed soil moisture of 28 stations ranged from 0.44 to 0.81 for 0-10 cm, and from 0.40 to 0.72 for 10-100 cm ( Figure 3). In general, the correlation coefficients of 0-10 cm were higher than that of 10-100 cm, and showed a significant increase with increased precipitation, while a similar trend was not detected for 10-100 cm. Three stations with a long-time span of observation data and remarkable difference in precipitation were selected to represent the typical arid, semiarid and semi-humid region, Water 2019, 11, 2241 9 of 24 as shown in Figure 4. The variability of SM VIC was generally lower than that of SM obs , which may be partly attributed to the difference in time-space scale. SM VIC represented the average daily state of the entire grid cell, while SM obs was observed at a certain sampling point at a certain time, which was more susceptible to the local transient climate, and the vegetation and soil conditions of the sampling site. Generally, the results of evaluation confirmed a reasonable and acceptable performance from the calibrated VIC model.    daily state of the entire grid cell, while SMobs was observed at a certain sampling point at a certain time, which was more susceptible to the local transient climate, and the vegetation and soil conditions of the sampling site. Generally, the results of evaluation confirmed a reasonable and acceptable performance from the calibrated VIC model.

Spatial Patterns
The spatial distribution patterns of ET, Runoff, SM1, and SM2 throughout the period 1984-2015 in the Loess Plateau were shown in Figure 5. The multi-year average of ET, SM1 and SM2 showed a decreasing gradient trend from southeast to northwest, with an average of 321.95 mm, 136.02 mm and 319.24 mm, respectively, which had a high consistency with the gradient spread of precipitation. However, the runoff, ranging from 45.40 mm to 160.69 mm, showed an obviously different spatial pattern, with a higher runoff in the central Loess Plateau and a relatively smaller runoff in the southeastern forest and irrigated agricultural area. in the Loess Plateau were shown in Figure 5. The multi-year average of ET, SM1 and SM2 showed a decreasing gradient trend from southeast to northwest, with an average of 321.95 mm, 136.02 mm and 319.24 mm, respectively, which had a high consistency with the gradient spread of precipitation. However, the runoff, ranging from 45.40 mm to 160.69 mm, showed an obviously different spatial pattern, with a higher runoff in the central Loess Plateau and a relatively smaller runoff in the southeastern forest and irrigated agricultural area.
The variation in four HVs between the periods 1984-1999 and 2000-2015, with a mean value of −1.55 mm, 4.28 mm, 1.32 mm and −3.94 mm, respectively, were derived by CC and LUCC, as shown in Figure 6. The ET increased in the central and southern Loess Plateau, including the Ziwu Mountains, and part of the middle reaches of the Yellow River where vegetation restoration efforts by the GGP have taken place since 1999 [67], and an obvious decrease was found in the southeastern and southwestern Loess plateau. The variations of Runoff and SM2 showed a similar spatial pattern, with a larger fluctuation in the central Loess Plateau, and the SM2 experienced a sharp decrease, up to more than 60 mm, which demonstrates a great decline in terms of the regional multi-year average SM2 of a grid cell. The variation of SM1 ranged from -20 mm to 20 mm for most of the regions, but exhibited a larger spatial heterogeneity. The HVs variations of the six catchments were discussed in detail in Section 3.4.  The variation in four HVs between the periods 1984-1999 and 2000-2015, with a mean value of −1.55 mm, 4.28 mm, 1.32 mm and −3.94 mm, respectively, were derived by CC and LUCC, as shown in Figure 6. The ET increased in the central and southern Loess Plateau, including the Ziwu Mountains, and part of the middle reaches of the Yellow River where vegetation restoration efforts by the GGP have taken place since 1999 [67], and an obvious decrease was found in the southeastern and southwestern Loess plateau. The variations of Runoff and SM2 showed a similar spatial pattern, with a larger fluctuation in the central Loess Plateau, and the SM2 experienced a sharp decrease, up to more than 60 mm, which demonstrates a great decline in terms of the regional multi-year average SM2 of a grid cell. The variation of SM1 ranged from -20 mm to 20 mm for most of the regions, but exhibited a larger spatial heterogeneity. The HVs variations of the six catchments were discussed in detail in Section 3.4. The minimum values of the four HVs were found at 1997, with a minimum precipitation of 331.02 mm and maximum values were found at 2003, with a minimum precipitation of 575.14 mm. The fluctuation of the four HVs exhibited obvious uniformity in precipitation change, which implies that precipitation was the main driving factor for the change of HVs on the regional scale of the Loess Plateau. Among the four HVs, the runoff was most sensitive to precipitation changes, with a change rate of 67.23%, and the SM2 showed the maximum resistance to precipitation changes, with a change rate of 10.01%, which suggests that deep soil water can effectively resist the impact of precipitation changes. The minimum values of the four HVs were found at 1997, with a minimum precipitation of 331.02 mm and maximum values were found at 2003, with a minimum precipitation of 575.14 mm. The fluctuation of the four HVs exhibited obvious uniformity in precipitation change, which implies that precipitation was the main driving factor for the change of HVs on the regional scale of the Loess Plateau. Among the four HVs, the runoff was most sensitive to precipitation changes, with a change rate of 67.23%, and the SM2 showed the maximum resistance to precipitation changes, with a change rate of 10.01%, which suggests that deep soil water can effectively resist the impact of precipitation changes.

Climate Change (CC)
The relationship between HVs and precipitation on a regional scale in the Loess Plateau were shown in the previous section. In this section, which aims to explore the difference in the HV response to climate change of three vegetation types, the relationship between HVs and meteorological variables (precipitation and temperature) is discussed, based on the Scenario2 (S2) simulation and selected typical grid cells of cropland (uc), grassland (ug) and woodland (uw) (see Table S1 (Supplementary Materials) for more details). Wind speed was one of the key meteorological variables for VIC model, although, due to its large spatial and temporal variability, the effect of wind was not analyzed in this study. As shown in Table 2, a significant positive correlation between the four HVs and precipitation in the three vegetation types was confirmed, with r ranging from 0.76 to 0.97 at the significance level of p < 0.01. The Runoff had the maximum correlation, followed by ET, SM1, and SM2.

Climate Change (CC)
The relationship between HVs and precipitation on a regional scale in the Loess Plateau were shown in the previous section. In this section, which aims to explore the difference in the HV response to climate change of three vegetation types, the relationship between HVs and meteorological variables (precipitation and temperature) is discussed, based on the Scenario2 (S2) simulation and selected typical grid cells of cropland (uc), grassland (ug) and woodland (uw) (see Table S1 (Supplementary Materials) for more details). Wind speed was one of the key meteorological variables for VIC model, although, due to its large spatial and temporal variability, the effect of wind was not analyzed in this study. As shown in Table 2, a significant positive correlation between the four HVs and precipitation in the three vegetation types was confirmed, with r ranging from 0.76 to 0.97 at the significance level of p < 0.01. The Runoff had the maximum correlation, followed by ET, SM1, and SM2. The ET and Runoff were not related to temperature, while SM1 and SM2 showed a negative correlation, with r ranging from −0.29 to −0.55, but this was not significant except in the SM2 of woodland (uw), which indicated a lesser effect of temperature on hydrologic processes than in cropland and grassland. The variation in precipitation and temperature and corresponding variation in HVs between the periods 1984-1999 and 2000-2015 were analyzed to detect the driving role of climate change in three vegetation types (Table 3). The result further defined that the change in precipitation contributed the most to the variations in four HVs in cropland, with R 2 ranging from 0.58 to 0.74 (p < 0.01), which suggests that precipitation was the main and direct driving factor of cropland's hydrologic processes. The variations of Runoff and SM2 were also significantly correlated with variation in precipitation in grassland (R 2 = 0.76 and 0.64, p < 0.01) and woodland (R 2 = 0.92 and 0.72, p < 0.01), while a similarly significant role of precipitation was not found in ET and SM1, indicating a more complex response mechanism of ET and SM1 in grassland and woodland. Among the three vegetation types, precipitation had the strongest driving effect on Runoff and SM2 in woodland, followed by its effects on grassland and cropland. In comparison to precipitation, the effect of temperature changes on the hydrologic processes was very weak, with R 2 less than 0.12.
Since the HVs were significantly correlated with precipitation, the regression coefficients were calculated as the response level of HVs to precipitation changes, i.e., the variation in HVs induced by precipitation changing 1 mm ( Figure 8). The results showed that ET and Runoff were more sensitive to variation in precipitation than SM1 and SM2, indicating that most of the increased precipitation was converted into an increase in ET and Runoff. The ET of grassland exhibited the highest response, with a mean value of 0.51 mm/mm, followed by cropland (0.30 mm/mm) and woodland (0.21 mm/mm), which suggested that grassland could convert more precipitation increments into evapotranspiration, and SM1 was consistent with response of ET. However, the response levels of Runoff and SM2 to precipitation were ordered from high to low: woodland (0.36 mm/mm, 0.16 mm/mm) > cropland (0.29 mm/mm, 0.15 mm/mm) > grassland (0.22 mm/mm, 0.10 mm/mm).

Land Use and Land Cover Change (LUCC)
As shown in Table 1, the changes in hydrologic processes in Scenario0 (S0) and Scenario2 (S2) for the period 2000-2015 were considered to be induced by LUCC. Therefore, the differences in the annual mean HVs of three kinds of conversion grid cells (cg, cw, gw) (see Table S1 for more details) between S0 and S2 were employed, to analyze how the three typical vegetation conversion modes affected hydrologic processes. This was similar to the experiment designed for sensitivity analysis, as mentioned in Mao et al. [68], but was conducted under real climate, vegetation and soil conditions, which were more realistic and convincing, to reproduce the response of HVs to different vegetation

Land Use and Land Cover Change (LUCC)
As shown in Table 1, the changes in hydrologic processes in Scenario0 (S0) and Scenario2 (S2) for the period 2000-2015 were considered to be induced by LUCC. Therefore, the differences in the annual mean HVs of three kinds of conversion grid cells (cg, cw, gw) (see Table S1 for more details) between S0 and S2 were employed, to analyze how the three typical vegetation conversion modes affected hydrologic processes. This was similar to the experiment designed for sensitivity analysis, as mentioned in Mao et al. [68], but was conducted under real climate, vegetation and soil conditions, which were more realistic and convincing, to reproduce the response of HVs to different vegetation conversions. As shown in Figure 9a, the vegetation conversions had a greater impact on ET and SM2 than Runoff and SM1, especially for the conversion of cropland and grassland to woodland. Runoff experienced a slight decline after the three conversions, with a mean value of 9.14 mm, 14.31 mm, 10.18 mm for cg, cw and gw, respectively, while the variation in SM1 exhibited an obvious difference, fluctuating around 0, indicating a large spatial heterogeneity and complexity in shallow soil moisture response to vegetation conversion, which may closely relate to the degree of change in vegetation and regional climate, especially the temporal and spatial distribution of precipitation. The variations in ET and SM2, induced by the three vegetation conversions, were similar but opposite, which implied a close relationship in the water flux between ET and SM2. The ET increased most, by an average of 43.05 mm, and SM2 decreased by 33.32 mm in cw, followed by gw, with ET increasing by 28. followed by cg (2.08%) and gw (1.69%), and cg took more reduction in Runoff to supply SM1 than SM2, while this was comparable in cw and gw. Additionally, the effect of vegetation conversion on ET was relatively small and negligible during precipitation, which was likely, as evapotranspiration was close to saturated surface evapotranspiration, regardless of vegetation type. However, in nonrainfall days, ET experienced significant increases, in the order, from high to low, of cw (13.57%) > gw (11.12%) > cg (4.85%), which was mainly caused by the stronger evapotranspiration capacity of tree, compared to grass and crop. The results suggested that all three conversions had a reducing effect on runoff, and the conversions of cropland and grassland to woodland exacerbated the consumption of deep soil moisture to meet the increase in ET, causing an obvious decline in deep soil moisture. Assuming the reduction of SM2 was entirely evaporated, we can attribute 71.28% ± 18.64%, 65.89% ± 24.14% of the ET increase to the water loss of SM2 in cw and gw, respectively.

Contributions of Climate Change (CC) and Land Cover Change (LUCC)
In practice, hydrologic processes were concurrently affected by both climate change and vegetation changes in the Loess Plateau. The effects of different vegetation conversions and changes in meteorological variables, which occurred in a grid cell, may be strengthened or offset. The effects of vegetation change on hydrological processes were limited within the changing grid cells, while climate change had a wide-ranging impact on the entire basin and the region, though spatial heterogeneity may exist. This paper aimed to detect the separated effect of CC and LUCC, and their contribution to the variations in HVs.
Compared with detrended meteorological variables during the period 2000-2015, precipitation on the Loess Plateau showed a slight decline by 3.35 mm, with a smaller variability in annual mean precipitation, and the variation in the precipitation of each catchment was, in order from high to low, BL (11.25 mm) > JH (10.04 mm) > FH (4 mm) > WD (2.21 mm) > KY (−6.87 mm) > ZL (−13.5 mm) (Figure 10a). The land cover in the Loess Plateau has experienced obvious changes since 1999, and the change rates of cropland, grassland, and woodland between LULC2010 and LULC1990 were shown in Figure 10b. Benefiting from GGP, there was a notable increase in woodland in the Loess Plateau, in six catchments with varying degrees, whereas grassland and cropland significantly decreased, expecting in JH and ZL. Figure 10c-f showed the variations in four HVs, caused by different factors. In general, the ET and SM1 increased by an annual mean value of 2.4 mm and 0.18 mm, while Runoff and SM2 decreased by 1.91 mm and 2.77 mm in the Loess Plateau, influenced by the coeffects of CC and LUCC, and the variations in HVs induced by LUCC were higher than those Land cover changes influenced hydrological cycle by affecting the distribution of precipitation (mainly runoff and soil moisture) on rainfall days and the consumption of soil moisture (mainly ET) on non-rainfall days. In order to give a more detailed illustration of the effect of vegetation conversions on hydrological processes, we calculated the change rate of distribution of precipitation and consumption of soil moisture from June to September, noted that the change rates of ET, Runoff, SM1, and SM2 were calculated from the difference between S0 and S2 on rainfall days, and daily precipitation, while change rates of ET in non-rainfall days (ETnr) were calculated from the difference of ET in S0 and S2 to daily ET in S0 in non-rainfall days (Figure 9b). Results showed that all three vegetation conversion modes reduced runoff and caused more precipitation infiltration into soil, increasing soil moisture. Out of the three, cw exhibited the biggest drop, up to 2.62% in average, followed by cg (2.08%) and gw (1.69%), and cg took more reduction in Runoff to supply SM1 than SM2, while this was comparable in cw and gw. Additionally, the effect of vegetation conversion on ET was relatively small and negligible during precipitation, which was likely, as evapotranspiration was close to saturated surface evapotranspiration, regardless of vegetation type. However, in non-rainfall days, ET experienced significant increases, in the order, from high to low, of cw (13.57%) > gw (11.12%) > cg (4.85%), which was mainly caused by the stronger evapotranspiration capacity of tree, compared to grass and crop. The results suggested that all three conversions had a reducing effect on runoff, and the conversions of cropland and grassland to woodland exacerbated the consumption of deep soil moisture to meet the increase in ET, causing an obvious decline in deep soil moisture. Assuming the reduction of SM2 was entirely evaporated, we can attribute 71.28% ± 18.64%, 65.89% ± 24.14% of the ET increase to the water loss of SM2 in cw and gw, respectively.

Contributions of Climate Change (CC) and Land Cover Change (LUCC)
In practice, hydrologic processes were concurrently affected by both climate change and vegetation changes in the Loess Plateau. The effects of different vegetation conversions and changes in meteorological variables, which occurred in a grid cell, may be strengthened or offset. The effects of vegetation change on hydrological processes were limited within the changing grid cells, while climate change had a wide-ranging impact on the entire basin and the region, though spatial heterogeneity may exist. This paper aimed to detect the separated effect of CC and LUCC, and their contribution to the variations in HVs.
Compared with detrended meteorological variables during the period 2000-2015, precipitation on the Loess Plateau showed a slight decline by 3.35 mm, with a smaller variability in annual mean precipitation, and the variation in the precipitation of each catchment was, in order from high to low, BL (11.25 mm) > JH (10.04 mm) > FH (4 mm) > WD (2.21 mm) > KY (−6.87 mm) > ZL (−13.5 mm) (Figure 10a). The land cover in the Loess Plateau has experienced obvious changes since 1999, and the change rates of cropland, grassland, and woodland between LULC2010 and LULC1990 were shown in Figure 10b. Benefiting from GGP, there was a notable increase in woodland in the Loess Plateau, in six catchments with varying degrees, whereas grassland and cropland significantly decreased, expecting in JH and ZL. Figure 10c-f showed the variations in four HVs, caused by different factors. In general, the ET and SM1 increased by an annual mean value of 2.4 mm and 0.18 mm, while Runoff and SM2 decreased by 1.91 mm and 2.77 mm in the Loess Plateau, influenced by the coeffects of CC and LUCC, and the variations in HVs induced by LUCC were higher than those induced by CC, excepting SM1. The variations in HVs caused by CC were consistent with precipitation change, indicating that the increase in precipitation was a direct driving factor in the increase in HVs. The impact of vegetation changes on each hydrological variable was relatively uniform, with some variation in degrees, i.e., increasing ET and SM1 and decreasing Runoff and SM2, which mainly resulted from the conversion of cropland and grassland to woodland, expect Runoff in ZL.
Due to differences in climate change intensity and vegetation conversion modes, the impact of CC and LUCC on the HVs of Loess Plateau and six catchments exhibited obvious differences. The variation in HVs in the FH, where woodland increased most (2.77%), with a slight increase in precipitation, was dominated by vegetation change, followed by WD, whereas the hydrologic processes in ZL were more influenced by climate change, with an obvious and larger decrease in precipitation detected alongside the major vegetation conversion of grassland to cropland, which had a relatively small impact on hydrological processes. However, the variation in the four HVs in BL and JH was dominated by different factor; vegetation change contributed more to the increase in ET and decrease in Runoff, while climate change dominated variations in SM1 and SM2. Differences in the response level of the four HVs to changes in precipitation, as discussed in Section 3.3.2, likely to account for this. Among the six catchments, KY exhibited the greatest complexity, with comparable contributions of CC and LUCC in absolute values of variation in the four HVs, while the increase in other land covers may serve as a key uncertainty to influence hydrologic processes, which was not analyzed in this paper. The results suggest that the responses of the four HVs to CC and LUCC were significantly different in the six catchments, which may be related to the conversion mode of vegetation and the intensity of climate change.
induced by CC, excepting SM1. The variations in HVs caused by CC were consistent with precipitation change, indicating that the increase in precipitation was a direct driving factor in the increase in HVs. The impact of vegetation changes on each hydrological variable was relatively uniform, with some variation in degrees, i.e., increasing ET and SM1 and decreasing Runoff and SM2, which mainly resulted from the conversion of cropland and grassland to woodland, expect Runoff in ZL. Due to differences in climate change intensity and vegetation conversion modes, the impact of CC and LUCC on the HVs of Loess Plateau and six catchments exhibited obvious differences. The variation in HVs in the FH, where woodland increased most (2.77%), with a slight increase in precipitation, was dominated by vegetation change, followed by WD, whereas the hydrologic processes in ZL were more influenced by climate change, with an obvious and larger decrease in precipitation detected alongside the major vegetation conversion of grassland to cropland, which had a relatively small impact on hydrological processes. However, the variation in the four HVs in BL and JH was dominated by different factor; vegetation change contributed more to the increase in ET and decrease in Runoff, while climate change dominated variations in SM1 and SM2. Differences in the response level of the four HVs to changes in precipitation, as discussed in Section 3.3.2, likely to account for this. Among the six catchments, KY exhibited the greatest complexity, with comparable contributions of CC and LUCC in absolute values of variation in the four HVs, while the increase in other land covers may serve as a key uncertainty to influence hydrologic processes, which was not analyzed in this paper. The results suggest that the responses of the four HVs to CC and LUCC were significantly different in the six catchments, which may be related to the conversion mode of vegetation and the intensity of climate change.

Hydrological Effect of Vegetation Conversion
Vegetation is critical in the hydrological cycle [69,70]. Disturbance, both natural (e.g., wildfire, insect outbreaks, disease, windstorms, drought) and anthropogenic (e.g., timber harvesting, land conversion), can have a profound effect on hydrological processes, through impact on vegetation dynamics. The Loess Plateau is covered by nearly 100 m loess in thickness and a loose soil structure [71], and the groundwater level is relatively deep, up to 30-100 m below surface [72], making it difficult to effectively replenish soil moisture. Little of the groundwater at these depths can be used as a supply for soil evaporation and plant transpiration [28], and precipitation acts as the only natural source of water. Vegetation plays a redistributive role in precipitation, distributing precipitation into canopy interception, surface runoff, soil moisture, and subsurface runoff. Different vegetation classes have different precipitation allocation strategies, and vegetation change will lead to obvious changes in the precipitation distribution ratio between these hydrological variables [11,73,74]. Runoff is very sensitive to vegetation change. The study of Wang et al. [75] found that the conversion of grassland to forest decreased the annual water yield from slope by 50-100 mm, which was further confirmed in our study. Increased canopy interception by woodland accounted for a fraction of the decline in runoff, while a decrease in runoff infiltrates the soil and becomes soil moisture, especially during high-intensity rainfall, when water can replenish deeper soil. However, many studies have demonstrated that artificial trees can consume more water than natural native species [30,33,75]. Higher evapotranspiration of forests consumes more soil moisture than crop and grass and the increase in soil moisture was insufficient to support the increase in ET, which resulted in a decline in deep soil moisture [14,31,32,34]. A similar conclusion was found by Feng et al. [26], who detected a systematic, significant and negative trend in the soil moisture of Loess Plateau (p < 0.001) based on data estimated by the Community Land Model GLDAS, and confirmed that the ecosystem was mining soil water. These results were consistent with the reported soil desiccation in the Loess Plateau that had been attributed to the excessive depletion of deep soil water by planted vegetation and long-term insufficient precipitation input [76].
Our results suggested that the conversion of cropland and grassland to woodland, closely related to GGP, obviously changes hydrologic processes by increasing ET and decreasing runoff and deep soil moisture, which exacerbates soil water consumption and causes soil drying. The implementation of GGP aims to effectively reduce soil erosion, alleviate land degradation and improve ecosystem services, while soil drought, which was induced by more consumption of soil water to meet increased ET, and may directly cause ecosystem degradation, was not expected. By contrast, the conversion of cropland to grassland reduced runoff and increased soil moisture overall, which has a better effect on controlling soil erosion and water conservation.

Uncertainties and Limitations
The uncertainties of this study mainly came from two factors, i.e., meteorological forcing data and parameters. The meteorological forcing dataset was generated by employing the Synagraphic Mapping System (SYMAP) interpolation algorithm and based on observation data of 112 meteorological stations. Though we selected as many stations as possible, to ensure the accuracy of meteorological data, uncertainties were inevitably introduced in the interpolation from point to face, especially for wind speed, which has great and unpredictable variability in both time and space. Precipitation proved to be the most critical input for hydrological models and a tiny deviation is likely to result in great uncertainties. Parameters after calibration and validation were considered reasonable and applicable to simulate hydrologic processes in the catchments. Using a parameter dataset (i.e., six parameters with calibration) to characterize all grid cells of a catchment was a generalized description of the catchment features. Furthermore, parameter implementation, which was actually based on an idealized assumption for simulation in uncalibrated regions, was essential for hydrological simulation from catchment-scale up to region-scale [57,77]. In this study, though further validations were conducted in some catchments to ensure the reliability of parameter transplantation, and obtained a relatively satisfactory result, a potential uncertainty still existed, which may directly influence the model performance and call for improved parameterization.
In addition to these uncertainties, there were several limitations that need to be further examined for consideration and improvement in subsequent studies. The version of VIC model used in this study employed static vegetation parameters for the whole study period. In order to detect the impact of vegetation change on hydrologic processes, two vegetation parameter datasets were applied to characterize the land cover of Loess Plateau before and after the implementation of GGP. However, vegetation change was a continuous and serial process accompanied by GGP, and this application was considered a less than ideal alternative, with the assumption that land cover changes were completed in an instant at 1999. Thus, the changes of the four HVs over the period 1984-2015 were actually changes in HVs driven by climate change in the land covers of two periods before and after 1999, which were insufficient to reflect the actual trend of HVs under the influence of CC and LUCC. Therefore, to avoid this potential bias, we employed the variation in HVs between the two scenario simulations, rather than the trend, to detect the impact of CC and LUCC. An improvement in VIC model, updating the land cover information in simulations, which depends on yearly dynamic vegetation parameters (i.e., vegetation classes and LAI) that can be retrieved from remote sensing data, was requisite to solve this problem thoroughly. Besides, scenario simulation methods were widely used in the study of the effect of CC and LUCC on hydrological processes. As designed in this study, two vegetation parameter datasets and two forcing datasets were combined over different scenarios, to distinguish the impact of each factor. However, climate conditions and land cover were always closely connected, including feedback and response. For example, the effect of LUCC was separated by changing the land cover in a detrended climate condition, while the effect of LUCC was likely to be strengthened or weakened in the context of actual climate change. Several previous studies confirmed a close interaction between climate and land cover [78,79]. In this study, this interaction can be further confirmed by the fact that the sum of variations in each hydrological variable induced by CC and LUCC was not equal to that induced by the coeffect, indicating there was not a simple addition relationship between the effects of CC and LUCC. The scenario simulations unavoidably split this link and interaction [39], and inevitably introduce a certain bias in quantifying the variation in HVs, induced by CC and LUCC. Improvement of the VIC model by coupling dynamic vegetation processes in a changing climate background would be favorable to remedy this issue.

Conclusions
In this study, the VIC model was employed to simulate the water fluxes and states in the Loess Plateau from 1965 to 2015. With the evaluation of VIC simulations, using observed streamflow records and multi-source datasets of ET and soil moisture, VIC was confirmed to successfully capture the magnitude and inter-annual variability of streamflow, evapotranspiration, runoff, and soil moisture over the Loess Plateau. The following conclusions can be drawn:

1.
The ET and SM1 increased slightly during the period 1984-2015 over the entire Loess Plateau, while a decreasing trend was found in Runoff and SM2. The changes in hydrological variables across the Loess Plateau exhibits clear spatial heterogeneity, and the central loess plateau experienced an obvious increase in ET, and decrease in Runoff and SM2.

2.
The hydrological variables were significantly related to precipitation. The response levels of HVs to precipitation were different in the three vegetation classes.

3.
Vegetation conversions had a significant impact on hydrologic processes, especially conversion of cropland to woodland (cw), which causes a dramatic increase in ET and decrease in SM2. The three vegetation conversions decreased runoff and converted it to soil moisture, leading to an obvious increase in SM1 and SM2 on rainfall days, while ET, especially for cw and gw, was considerably higher than that before conversion on non-rainfall days.

4.
Land cover change contributed more to the variation in HVs in the Loess Plateau, while different factors were detected to control the change of hydrologic processes in six catchments. There existed a trade-off between mode and intensity of vegetation conversion versus degree of climate change.
Both climate change and vegetation conversions have been proved to have a notable impact on hydrologic processes, especially the mode and intensity of vegetation conversion. It is important to generate effective strategies to implement vegetation restoration scientifically in the context of climate change, including selecting a reasonable vegetation conversion mode and quantity, according to local conditions, which is critical to maintain a benign water cycle and achieve sustainable development of ecosystems.