*2.2. Hydrological Model*

The VIC model is a semi-distributed physically-based hydrological model. It can simulate water and energy balance. In the model, a study area is divided into grids according to latitudes and longitudes, and all the calculations are done in each grid separately. In each grid, all elements are computed based on different land cover types, before being averaged according to their correspondent fractions [25,26]. Since the model has been proved to perform well over a range of scales [30,31], we employed it in our study. Here, the spatial resolution of modeling is at 0.5◦ × 0.5◦ to ensure consistency with the climate forcing data.

## *2.3. Data Availability*

#### 2.3.1. Vegetation and Soil Parameters

To produce the vegetation parameters for VIC, land use maps and the average leaf area index (LAI) over 12 months are necessary. The two land use maps obtained were created by merging Landsat Thematic Mapper (TM) images for the periods ranging from 1983–1986 and 2010–2015. In these studies, land cover conditions were divided into 12 types as in Figure 1. A human–computer interactive interpretation method of remotely sensed land use cover information was used to interpret these maps, and the accuracy was over 90% [29,32]. The LAI datasets were available from the Global Land Surface Satellite (GLASS) products (http://www.bnu-datacenter.com/) that were retrieved from the Moderate Resolution Imaging Spectroradiometer reflectance data (MOD-09A1) using a general regression neural network algorithm [33]. The resulting combinations of land use and LAI datasets for the two collection periods (1983–1986 and 2010–2015) were referred to as "LC1986" and "LC2015," respectively. The other parameters, except for the LAI, in the VIC vegetation library data were set as in the study of a previous study [24]. We assumed that the differences between the two combination datasets could accurately represent the progress accomplished by the ecological programs and other factors from 1986 to 2015. The soil parameters were set according to those of a previous study [24], in which the parameters were derived based on the Food and Agriculture Organization of the United Nations (FAO) digital soil map of the world, and they were previously evaluated at a global scale [34,35].

#### 2.3.2. Bias-Corrected Climate Datasets

The climate data sets used in this study are from the Inter-Sectoral Impact Model Inter-comparison Project (ISI-MIP), and are available at: https://esg.pik-potsdam.de/. This project provides a framework to compare climate impact projections in different sectors and at different scales, which enable quantitative synthesis of climate change impacts at different level of global warming [27]. A new bias-correction method was developed within the first stage of this project, which can reduce the bias between daily or monthly simulated and observed climate data, while maintaining the absolute or relevant long-term trend much better than previous methods [36]. In this study, bias-corrected climate data from five GCMs (HadGEM2-ES, GFDL-ESM2M, IPSL-CM5A-LR, MIROC-ESM-CHEM, and NorESM1-M) with a daily time step were selected in order to lessen the uncertainty of a single model. Considering that the RCP8.5 scenario, in which the anthropogenic radiative forcing will approach 8.5W m<sup>−</sup><sup>2</sup> by 2100 and the concentration of CO2 will be 3–4 times the present value, estimates extremely severe circumstances in future, evaluation of this scenario was assumed to be more important than for others scenarios [37]. Five meteorological datasets were provided under this scenario: precipitation, maximum temperature, minimum temperature, mean wind speed, and downward shortwave radiation at a 0.5◦ × 0.5◦ spatial resolution, daily time steps from 2006 to 2099 were collected to generate the climate forcing data. These data are accessible at https: //www.isimip.org/outputdata/caveats-fast-track/. Among these data, precipitation, temperature, and wind speed are necessary to force the VIC model, and the downward shortwave radiation plays an important role in detecting the effects of land cover on the hydrological cycle, as it will be redistributed differently according to the various vegetation types.

## *2.4. Experimental Design*

In this study, the VIC model was run based on climate forcing datasets from the five GCMs, as described in Section 2.3.2, and the simulation results were averaged to reduce the uncertainties from the forcing data. Moreover, the model parameters described in Section 2.3 have been extensively validated by a previous study [24], with 15 stations of streamflow data and 10 stations of evapotranspiration (ET) data. The model has shown a favorable performance, with an average Nash–Sutcliffe efficiency for streamflow of 0.55 and a Pearson correlation coefficient for the ET of greater than 0.5 [24]. Therefore, we did not further validate the model parameters.

To identify the impact of historical LCC on the future hydrological regime in the TNR, the present land cover state (corresponding to the LC2015 data) was assumed to remain constant throughout the study period. Thus, two simulation experiments were conducted with different land cover and vegetation input data. The first experiment involved running the model with the LC2015 data to explore changes in the hydrological cycle under the influence of climate change alone. We paid the most attention to the ET, runoff (R), and soil moisture (SM), as these three hydrologic features play important roles in the water cycle, and their values can reflect the ecological conditions of the area to some degree. The time series of the climatic (i.e., precipitation and temperature) and the three hydrological variables were divided into four periods: 2020–2039, 2040–2059, 2060–2079, and 2080–2099 to lessen the uncertainties of a single year in the analysis. The spatial distribution for the four periods, as well as the annual mean values over the entire region, were analyzed to determine the hydrological response to climate change across the TNR.

The second simulation experiment ran the model with the LC1986 data for the same period as in the first simulation. Therefore, the two experiments differ only with respect to their land cover and vegetation data, respectively regarded as the land cover conditions after afforestation (LC2015) and before afforestation (LC1986). Comparison of the two experiments can be used to identify the effects of past LCC on the future hydrological regime under RCP8.5 scenario. The comparative analyses between the two experiments were conducted over different basins and periods. The annual differences and the spatial distribution of the annual mean differences between the two experiments were evaluated to measure the effects of afforestation on the hydrology across the TNR. We also analyzed the differences in four seasons, so as to explore the different effects of LCC in different seasons.
