Impact of revegetation of the Loess Plateau of China on the regional growing season water balance

To resolve a series of ecological and environmental problems over the Loess Plateau, the “Grain for Green Program” (GFGP) was initiated at the end of 1990s. Following the conversion of croplands and bare land on hillslopes to forests, the Loess Plateau has displayed a significant greening trend, which has resulted in soil erosion being reduced. However, the GFGP has also affected the hydrology of the Loess Plateau, which has raised questions regarding whether the GFGP should be continued in the future. We investigated the impact of revegetation on the hydrology of the Loess Plateau using relatively high-resolution simulations and multiple realizations with the Weather Research and Forecasting (WRF) model. Results suggest that revegetation since the launch of the GFGP has reduced runoff and soil moisture due to enhanced evapotranspiration. Further revegetation associated with the GFGP policy is likely to further increase evapotranspiration, and thereby reduce runoff and soil moisture. The increase in evapotranspiration is associated with biophysical changes, including deeper roots that deplete deep soil moisture stores. However, despite the increase in evapotranspiration, our results show no impact on rainfall. Our study cautions against further revegetation over the Loess Plateau given the reduction in water available for agriculture and human settlements and the lack of any significant compensation from rainfall.


Introduction
The Loess Plateau is a highland region of north central China, covering about 640 000 km 2 . The loess soils are well suited for agriculture, so natural forests have been progressively converted to farmland to support the growing population over the last 7000 years . However, the loess is also prone to wind and water erosion; thus, the area's long history of deforestation is associated with soil erosion, resulting in land degradation, low agricultural productivity, and significant local poverty in some farming communities (Bryan et al., 2018;Chen et al., 2015;Fu et al., 2017). The soil erosion aggravates the flux of sediment into the Yellow River Miao et al., 2010;Peng et al., 2010), increasing the risk of catastrophic flooding in some densely populated regions downstream (Bryan et al., 2018;Fu et al., 2017).
To minimize soil erosion, mitigate flood risk, store carbon, and improve livelihoods over the Loess Plateau, the "Grain for Green Program" (GFGP) was initiated by reforesting hillslopes in the late 1990s (Bryan et al., 2018;Fu et al., 2017;Liu et al., 2008). Consequently, the Loess Plateau has displayed a significant "greening" trend Fu et al., 2017;Li et al., 2017). The large-scale vegetation restoration program has also reduced soil erosion over the Published by Copernicus Publications on behalf of the European Geosciences Union.
As a consequence of the beneficial outcomes of the GFGP, further investment is planned with a commitment of around USD 33.9 billion from China through to 2050 . However, further revegetation over the Loess Plateau is controversial Chen et al., 2015;Fu et al., 2017) with evidence from field (Jia et al., 2017;Jin et al., 2011;Wang et al., 2012) and satellite (Feng et al., 2017;Lv et al., 2019a;Xiao, 2014) observations that revegetation has affected the hydrological balance of the region. Compared with croplands or barren surfaces, the planted forests enable higher evapotranspiration associated with a larger leaf area, higher aerodynamic roughness, and deeper roots (Anderson et al., 2011;Bonan, 2008;Bright et al., 2015). Consequently, revegetation tends to decrease soil moisture and runoff with the associated risk of limiting water availability for agriculture, human consumption, and industry Chen et al., 2015;Fu et al., 2017). Indeed, the present vegetation over the Loess Plateau, which to some extent reflects decades of reafforestation, may already exceed the limit that the local water supply can support; hence, further revegetation may not be sustainable S. L. Zhang et al., 2018).
Despite the increasing observational evidence demonstrating that revegetation tends to impair the hydrological balance of the Loess Plateau, the response of rainfall to revegetation over this region has commonly been overlooked. This is mainly due to the difficulty in detecting the impact of revegetation on rainfall from observations. As an important component of hydrological cycle of the Loess Plateau, rainfall not only controls the terrestrial water budget but also influences soil erosion and the discharge of sediment into the Yellow River (Liang et al., 2015;Miao et al., 2010;Peng et al., 2010;Wang et al., 2016). Therefore, information on how rainfall responds to revegetation is critical for a comprehensive assessment of the impact of revegetation on the hydrology of the region. Indeed, if rainfall responds to revegetation, this may influence national policies on whether to continue large-scale vegetation restoration programs. Afforestation or deforestation does have the potential to affect rainfall via changes in biogeophysical processes, but any impact of afforestation or deforestation on rainfall tends to be highly region specific (Findell et al., 2006;Lorenz et al., 2016;Winckler et al., 2017).
In contrast with observations, modeling can help disentangle the impact of revegetation on rainfall from the impact of other drivers. Cao et al. (2017) and Li et al. (2018) performed numerical experiments over the whole of China and demonstrated that the revegetation over the Loess Plateau can enhance rainfall locally. Very recently, Lv et al. (2019b) and Cao et al. (2019) performed simulations focused on the Loess Plateau in order to examine the impact of revegetation or afforestation on rainfall. Lv et al. (2019b) reported a significant increase in rainfall, while Cao et al. (2019) found spatially divergent changes in rainfall. We also note some earlier studies investigating the response of rainfall to land cover change across China (e.g., Chen et al., 2017;Ma et al., 2013;Wang et al., 2014). Unfortunately, these studies either focused less on the Loess Plateau (Ma et al., 2013) or applied land cover changes unable to reflect the revegetation of the Loess Plateau (Chen et al., 2017;Wang et al., 2014). Therefore, large uncertainties remain in the response of rainfall to the revegetation of the Loess Plateau owing to inconsistent conclusions derived from limited studies. We note that  reported that the increased rainfall due to revegetation over North China (covering but not limited to the Loess Plateau) was large enough to compensate for the increase in evapotranspiration and resulted in little impact on soil moisture. This simulated negligible soil moisture change associated with revegetation is contradicted by extensive studies based on observations (e.g., Feng et al., 2017;Jia et al., 2017;Wang et al., 2012). Here, we note that it might be unfair to directly compare the observational and modeling results because observational results commonly incorporate multiple factors and modeling results are subject to uncertainties in both land cover change and biophysical parametrization schemes implemented in models (de Noblet-Ducoudre et al., 2012;Pitman et al., 2009). These intrinsic differences between observational and modeling results cannot fully account for the disagreement on the runoff and soil moisture change due to revegetation over the Loess Plateau. Thus, the impact of revegetation on the hydrology of the Loess Plateau remains unclear and needs careful reevaluation.
In this study, we examine the impact of revegetation following the launch of the GFGP on the hydrology of the Loess Plateau using relatively high-resolution simulations with the Weather Research and Forecasting model. We also examine the impact of further revegetation on the hydrology of the Loess Plateau with the goal of providing helpful information to policymakers. As far as we know, there has been no study investigating how the regional hydrology would be affected by further revegetation over the Loess Plateau, which is an important factor for informing policymakers on the mitigation and adaptation of climate change for this region. Additionally, the vegetation over the Loess Plateau is fragile and highly dependent on water availability . How the hydrology would be impacted by further revegetation determines the water availability and, in turn, how much more revegetation can be sustained over the Loess Plateau. Neglecting this process risks errors in assessing the upper threshold for vegetation in the Loess Plateau S. L. Zhang et al., 2018). Given the importance of revegetation over the Loess Plateau now and in the future, we examine the impact of further revegetation on the hydrology of the Loess Plateau and pay particular attention to the response of rainfall to revegetation. J. Ge et al.: Impact of revegetation of the Loess Plateau on the regional growing season water balance 517 2 Methods

Model configuration
The Weather Research and Forecasting (WRF, version 3.9.1.1, Skamarock et al., 2008) model, a fully coupled land-atmosphere regional weather and climate model, was used in our study. WRF has been shown to perform well in the dynamic downscaling of regional climate over China (e.g., He et al., 2017;Sato and Xue, 2013;Yu et al., 2015). Additionally, WRF has been used to study the impact of land use and land cover change on the hydrological balance at regional scales (Deng et al., 2015;. Thus, while WRF is potentially suitable for evaluating the impact of revegetation on the hydrology of the Loess Plateau, we undertake an evaluation of WRF in simulating surface air temperature and rainfall for this region (see Sect. 3.1). To perform simulations at a high spatial resolution over the Loess Plateau region, we applied two-way nested runs with two domains at different grid resolutions running simultaneously. The ERA-Interim reanalysis data (Dee et al., 2011; Table 1) provided the boundary conditions for the larger and coarserresolution (30 km) domain, and the larger domain provided boundary conditions for the smaller and higher-resolution (10 km) domain. The ERA-Interim reanalysis data also provided the initial conditions for both domains. Using a Lambert projection, the larger domain was centered at 100 • E, 37 • N, with 180 grid points in the west-east direction and 155 grid points in the south-north direction, covering most of China and some of the surrounding regions (Fig. 1a). The inner domain covers the entire Loess Plateau with 166 grid points in the west-east direction and 151 grid points in the south-north direction (Fig. 1a, b). Both domains had 28 sigma levels in the vertical direction with the top level set at 70 hPa. Figure 1b shows the region analyzed in this paper.
The main physical parameterization schemes used in our study included the WRF single-moment 6-class microphysics scheme (Hong and Lim, 2006), the Dudhia shortwave radiation scheme (Dudhia, 1989), the Rapid Radiative Transfer Model (RRTM, Mlawer et al., 1997) for longwave radiation, a revised MM5 scheme (Jimenez et al., 2012) for the surface layer, the Noah land surface model (Ek, 2003), the Yonsei University scheme  for the planetary boundary layer, and the Kain-Fritsch scheme (Kain, 2004) for cumulus convection. The Noah land surface model used the Unified NCEP/NCAR/AFWA scheme with soil temperature and moisture in four layers (the first layer from 0 to 10 cm, the second layer from 10 to 40 cm, the third layer from 40 to 100 cm, and the fourth layer from 100 to 200 cm), fractional snow cover, and frozen soil physics. A sub-tiling option considering three land cover types within each grid cell was applied to help improve the simulations of the land surface fluxes and temperature (Li et al., 2013).

Satellite data
We used satellite-observed land cover type obtained from the Moderate Resolution Imaging Spectroradiometer (MODIS) Land Cover Type product (MCD12Q1, Version 6; Friedl and Sulla-Menashe, 2019; Table 1). This provides land cover types based on the "International Geosphere-Biosphere Programme" (IGBP) classification scheme (Table 2) globally at a spatial resolution of 500 m, and at yearly intervals from 2001 to 2017. The MCD12Q1 Version 6 is improved over previous versions via substantial improvements to algorithms, classification schemes, and spatial resolution (Sulla-Menashe et al., 2019). We changed the land cover type within the Loess Plateau while retaining the default land cover type for other regions in our experiments (see details in Sect. 2.3). Therefore, the MCD12Q1 data were reprojected to geographic grid data with a resolution of 30 s (approximately 0.9 km) by the MODIS Reprojection Tool to make them consistent with the default land cover map in WRF. Key land surface biogeophysical parameters include the green vegetation fraction (VEGFRA), snow free albedo (α), leaf area index (LAI), and the background roughness length (Z 0 ). The fraction of photosynthetically active radiation (FPAR) can be used as a proxy for VEGFRA (Kumar et al., 2014;Liu et al., 2006), enabling both VEGFRA and LAI data to be obtained from the MODIS Terra+Aqua LAI/FPAR product (MCD15A2H, Version 6;Myneni et al., 2015a;Table 1). This provides 8 d composite LAI and FPAR globally at a spatial resolution of 500 m from 4 July 2002. The MODIS Terra LAI/FPAR product (MOD15A2H, Version 6; Myneni et al., 2015b; Table 1) was also used to provide observations prior to 2002, as it started on 8 February 2000. Although MOD15A2H covers a longer time span, MCD15A2H is generally preferred. This is because only observations from the MODIS sensor on NASA's Terra satellite are used to generate MOD15A2H, but observations from sensors on both Terra and Aqua satellites are used for MCD15A2H. The MCD15A2H and MOD15A2H sinusoidal tile grid data were reprojected before use. The 8 d LAI and FPAR data were composited to monthly data to make them suitable for WRF.
As we only focus on the growing season (see Sect. 2.3.1), α can be assumed to be equivalent to satellite-observed snowfree albedo. The α data were derived from the blue sky albedo for shortwave provided by the Global Land Surface Satellite (GLASS) product (Liang and Liu, 2012; Table 1). This provides an 8 d composite albedo globally at a spatial resolution of 0.05 • from 1981 to present. Compared with the MODIS albedo product, the GLASS albedo product has a higher temporal resolution and more successfully captures the surface albedo variations . The 8 d α data were composited to monthly data.
The background roughness length (Z 0 ) was calculated following Eq. (1): where Z max and Z min are the land-cover-dependent maximum and minimum background roughness lengths, respectively, provided by lookup tables. VEGFRA, VEGFRA max , and VEGFRA min are the instantaneous, maximum, and minimum green vegetation fractions, respectively, which were calculated from satellite-observed VEGFRA (equal to FPAR) that would be implemented in WRF (see Sect. 2.3).

Observation data
To evaluate the WRF model performance with respect to simulating the surface air temperature and rainfall over the Loess Plateau, we used a gridded observation dataset developed by the National Meteorological Information Center of the China Meteorological Administration (Zhao et al., 2014; Table 1). The dataset provides monthly surface air temperature and rainfall at a spatial resolution of 0.5 • from 1961 to present and was produced by merging more than 2400 station observations across China using thin plate spline interpolation. The dataset has been widely used to analyze the surface air temperature and rainfall over the Loess Plateau Tang et al., 2018). To facilitate the comparison between simulations and observations, the observation data were bilinearly interpolated to the WRF inner domain grid. Savannas 9 Tree cover between 10 % and 30 % (canopy > 2 m).
Permanent Wetlands 11 Permanently inundated lands with between 30 % and 60 % water cover and greater than 10 % vegetated cover.
Cropland 12 At least 60 % of the area is cultivated cropland.
Urban and Built-up lands 13 At least 30 % impervious surface area, including building materials, asphalt, and vehicles.
Snow and Ice 15 At least 60 % of the area is covered by snow and ice for at least 10 months of the year.
Barren 16 At least 60 % of the area is non-vegetated and barren (sand, rock, soil) with less than 10 % vegetation.
Water Bodies 17 At least 60 % of the area is covered by permanent water bodies.

Experiment design 2.3.1 The impact of revegetation since the launch of the GFGP
To examine the impact of revegetation on the hydrology of the Loess Plateau since the launch of the GFGP, we conducted a control experiment (LC 2001 ) and a sensitivity experiment (LC 2015 ; Table 3). For LC 2001 , satellite-observed land cover type, VEGFRA, LAI, and α in 2001 were used to approximate the land cover type and land surface biogeophysical parameters before the launch of the GFGP. There is a 1-year gap between the launch of the GFGP (end of 1999) and 2001, but any bias introduced by this gap is small compared with the changes in land cover type and land surface biogeophysical parameters between 1999 and present. Satellite-observed land cover type, VEGFRA, LAI, and α in 2015, representing the current land cover type and land surface biogeophysical status, were used for the LC 2015 . Model configurations were identical for the LC 2001 and LC 2015 except for land cover type and land surface biogeophysical parameters. Thus, comparing the LC 2001 and LC 2015 isolates the impact of revegetation since the launch of the GFGP. We note that the difference between LC 2001 and LC 2015 should not be regarded as equivalent to the impact of GFGP for two reasons. First, actual changes in land cover type since the launch of the GFGP are highly spatially heterogeneous due to various anthropogenic activities, including GFGP, ir-  (Fig. 2a, c, e, g). In addition to the gain of forests (including evergreen needleleaf, evergreen broadleaf, deciduous needleleaf, deciduous broadleaf, and mixed forests) and savannas (including woody savannas and savannas), other changes in land cover type include the expansion of croplands (including croplands and cropland/natural vegetation mosaics) at the expense of grasslands and savannas (Fig. 2g). These increased croplands revealed by the MODIS land cover product, which seem unlikely, have been reported previously (Fan et al., 2015;Lv et al., 2019b), and are likely associated with expanded irrigation activities along the Yellow River (Fan et al., 2015;Zhai et al., 2015). Second, the observed VEGFRA, LAI, and α changes also incorporate other factors including improved agricultural management, climate variability, rising atmospheric CO 2 concentration, and nitrogen deposition Fan et al., 2015;Piao et al., 2015). As shown in Fig. 3a, c, e, and g, the biogeophysical changes are not strictly limited to the regions undergoing changes in land cover type. For example, the α decrease mostly occurs over grasslands in the northwest (Fig. 3e), where the land cover type rarely changes (Fig. 2c). This decreased α is attributed to increased precipitation as well as the restoration of grasslands benefiting from the Returning Rangeland to Grassland Program launched in 2003 over this region (Zhai et al., 2015). In contrast, the α change is negligible in the SLP and ELP owing to the combined effects of increased forests (Fig. 2a) and croplands (Fig. 2d). However, overall, the MCD12Q1 demonstrates a significant greening trend (increased VEGFRA, LAI, and Z 0 and decreased α) over the Loess Plateau since the launch of the GFGP (Fig. 3), which is spatially consistent with previous studies (e.g., Cao et al., 2019;Xiao, 2014;Zhai et al., 2015). Both LC 2001 and LC 2015 were run from 1 May to 30 September from 1996 to 2015 resulting in 20 realization members for LC 2001 and LC 2015 , respectively. We only run simulations for the growing season; any impact of revege-tation should be most apparent during the growing season given that over 70 % of the annual rainfall occurs over the Loess Plateau in this season Tang et al., 2018).

The impact of further revegetation on the Loess Plateau
If the GFGP is continued in the future, further revegetation could impact the hydrology of the Loess Plateau. Therefore, we conducted a third experiment (LC futr ) in which the coverage of forests was assumed to be at a maximum over the Loess Plateau following the policy of the GFGP (Table 3).
To maximize forests, we first assumed that all croplands and barren land on hillslopes were converted to forests. We then assumed that savannas or forests with low coverage (e.g., low VEGFRA) became dense forests. The land cover and land surface biogeophysical parameters for the LC futr were then constructed following two steps: 1. All cropland, barren, and savanna pixels on hillslopes (> 15 • ) were replaced by forest pixels over the Loess Plateau based on the land cover map from 2015. The slope was derived from the Shuttle Radar Topography Mission (SRTM version 2.0; Table 1) digital elevation model at a spatial resolution of 3 s (about 90 m). The pixel resolution of the land cover type was 30 s; thus, every land cover type pixel covered 100 (10 × 10) slope values. To maximize the revegetation, land cover type pixels with maximum slope values over 15 • were regarded as hillslopes. For a pixel to be changed, the forest class was determined by the class of neighboring forest pixels, considering the adaptation of planted trees to local climate. Using this strategy, forest pixels increased by 164 % and croplands pixels decreased by nearly half in the constructed land cover map compared with the land cover type in 2001, with most conversions occurring in SLP (Fig. 2b, h).  2. We constructed the VEGFRA, LAI, and α map in line with the land cover type constructed in step 1. For each forest class, we screened out the "dense forests" pixels with VEGFRA over the 95th percentile among the pixels labeled as the same forest class over the Loess Plateau. The monthly values of VEGFRA, LAI, and α of the "dense forest" pixels were calculated for each forest class. We then adjusted the monthly VEGFRA, LAI, and α of other "non-dense forests" pixels to the values of the "dense forests" pixels. Using this strategy, all forest pixels over the Loess Plateau were changed to more dense forest. Consequently, the Loess Plateau shows an amplified greening trend in LC futr , especially in SLP (Fig. 3b, d, f, h).
The LC futr was run from 1 May to 30 September from 1996 to 2015. Therefore, comparing LC 2001 and LC futr isolates the impact of further revegetation on the hydrology of the Loess Plateau.

Identification of the impact of revegetation
Model internal variability is defined as the difference between realization members where the only differences are the initial conditions. These differences result from nonlinearities in the model physics and dynamics (Giorgi and Bi, 2000;Christensen et al., 2001). This means that some differences between LC 2001 and LC 2015 (or LC futr ) will be caused by internal variability in addition to revegetation (Lorenz et al., 2016;Ge et al., 2019). To minimize the impact of internal model variability, we performed multiple simulations for the year 2001 by changing initial conditions. Specifi-cally, we carried out a pair of experiments named LCENS 2001 and LCENS 2015 (Table 3), which were the same as LC 2001 and LC 2015 except that LCENS 2001 and LCENS 2015 were only run for the year 2001 but initialized for each day between 21 and 30 April and ended on 30 September. This led to a total of 11 members (including the members with initial dates of 1 May in LC 2001 and LC 2015 ) for LCENS 2001 and LCENS 2015 , respectively. Comparing LCENS 2001 and LCENS 2015 , simulated changes were likely robust if the impact from revegetation was large and consistent relative to the differences caused by the change in the initial condition.
Results before 1 June was discarded as spin-up time in each simulation. Our analysis focuses on June, July, August, and September (JJAS) averages.

Local significance test
To test the statistical significance of the local impact of revegetation on the hydrology we calculate a grid point by grid point Student's t test. This tests the null hypothesis that the two groups of data are from independent random samples from normal distributions with equal means and equal but unknown variances. The local difference is regarded as statistically significant when the p value of the two-tailed t test passes the significance level of 95 %.

Evaluation of WRF's skill in simulating temperature and rainfall
We first evaluate WRF's simulation of surface 2 m air temperature (T 2) and rainfall (RAIN), the quantities with the most credible observations available over the Loess Plateau, by comparing the averaged value of the 11 members in LCENS 2001 with the observed values in 2001. After topographic correction (Zhao et al., 2008), WRF simulates T 2 over the Loess Plateau mostly within 2 • C of the observations (Fig. 4a, c, e), although there are small areas where WRF simulates warmer temperatures (by 4 • C) than the observations. The model also performs well with respect to simulating RAIN (Fig. 4b, d, f), including a region of higher observed rainfall from the southwest to the central Loess Plateau. The RAIN bias between the WRF simulations and the observations is below 0.5 mm d −1 for almost the entire Loess Plateau (Fig. 4f). Larger RAIN biases mostly occur around the eastern and southern borders of the Loess Plateau, most likely due to extremely complex topography in these locations. As we focus on the impact of land cover change on the hydrology of the region, the reasonable simulation of RAIN gives us confidence in the results from WRF, particularly in SLP.

Impacts on surface fluxes
We first examine the change in the land surface radiation budget, energy, and water fluxes, as these are directly impacted by changes in land cover type and the surface biogeophysical parameters. Comparing LC 2001 and LC 2015 (LC 2015 -LC 2001 ), land surface net radiation (R net ), latent heat flux (Q E ) and sensible heat flux (Q H ) changes mainly occur where the land cover type and land surface biogeophysical parameters are changed, suggesting a strong local effect on R net , Q E , and Q H . R net increases by around 5-20 W m −2 (Fig. 5a) over most of the region due to a reduction in α (Fig. 3e). While Q E increases by 10-30 W m −2 (Fig. 5c), Q H reduces by around 10 W m −2 (Fig. 5e), mostly in SLP and ELP as a result of increased VEGFRA, LAI, and Z 0 (Fig. 3a, c, g). Changes in R net and Q E are statistically significant at a 95 % confidence level over most of the region, but statistically significant changes in Q H are mostly limited to SLP and ELP (see the embedded subplots in each panel, Fig. 5a, c, e). As a consequence of further revegetation (LC futr -LC 2001 ), R net , Q E , and Q H changes are intensified (Fig. 5b, d, f), especially in SLP where large areas of croplands are converted to forest leading to large changes in land surface biogeophysical parameters in LC futr (Figs. 2, 3). Focusing on SLP, the increase in evapotranspiration (ET) is 0.49 mm d −1 between LC 2001 and LC 2015 (Fig. 6a). WRF simulates further water loss (0.85 mm d −1 ) through ET if the revegetation is continued in the future (Fig. 6c). For ELP, where relatively fewer croplands or barren areas can be further converted to forests in LC futr , the future ET increase is still considerable (0.72 mm d −1 ; Fig. 6b, d). The values of the regional mean ET change among the 20 members of LC 2015 -LC 2001 and LC futr -LC 2001 remain consistently positive over SLP and ELP. This indicates that the simulated higher ET is  Fig. 2. The first and second line members denote absolute and relative changes averaged by 20 members. The black asterisk denotes that the change is statistically significant at a 95 % confidence level using a two-tailed Student's t test.
a consistent result from WRF as a consequence of revegetation since the launch of the GFGP and is likely to be further strengthened by continued revegetation over the Loess Plateau.

Impacts on rainfall
Increased ET can contribute to the formation of clouds and rainfall; therefore, we examine whether this is the case for the Loess Plateau. The RAIN is composed of convective rainfall (RAINC), calculated by the cumulus convection scheme, and non-convective rainfall (RAINNC), calculated by microphysics scheme, in WRF. Thus we separate RAINC and RAINNC changes in addition to the RAIN change in Fig. 7. As for LC -LC 2001 , the change in RAIN is spatially heterogeneous, with an increase of up to 1.2 mm d −1 in small parts of the northeast and a decrease of around −1.0 mm d −1 along the southeast border of the Loess Plateau (Fig. 7a). The RAIN change is divided almost evenly between RAINC and RAINNC (Fig. 7c, e). However, most of the RAIN, RAINC, and RAINNC changes are not statistically significant. In terms of LC futr -LC 2001 , RAIN, RAINC, and RAINNC are not significantly changed by further revegetation (Fig. 7b, For both LC 2015 -LC 2001 and LC futr -LC 2001 , most RAIN changes seem to be randomly scattered around the Loess Plateau instead of being located coincident with SLP or ELP where land cover type, land surface biogeophysical parameters, and land surface fluxes are most strongly modified (Fig. 7a, b). In contrast, the RAIN change is negligible over SLP and ELP for both LC 2015 -LC 2001 and LC futr -LC 2001 (Figs. 6, 7). However, the RAIN change in individual realizations is not small, e.g., the RAIN change varies from −2.11 to 2.21 mm d −1 over the ELP for LC -LC 2001 (Fig. 6b). Thus, averaging the divergent RAIN changes among the 20 members causes a negligible RAIN change overall. This large variability in RAIN changes among the 20 members can be attributed to either different boundary conditions (background climate), which cause the impact of land cover change to diverge (Pitman et al., 2011), or model internal variability. This will be further analyzed in Sect. 3.6.

Impacts on runoff
As a consequence of the significant ET increase and negligible and statistically insignificant RAIN change, underground runoff (UDROFF) is reduced by up to 1.5 mm d −1 locally for LC -LC 2001 (Fig. 8c). Averaged over the SLP and ELP, the UDROFF decreases by 0.16 mm d −1 (−23 %) and 0.34 mm d −1 (−23 %) for SLP and ELP, respectively (Fig. 6a, b). These UDROFF changes are not statistically significant and vary strongly among the 20 members, suggesting a large uncertainty in the UDROFF change. WRF Figure 7. Same as in Fig. 5 but for (a, b) total rainfall (mm d −1 ), (c, d) convective rainfall (mm d −1 ), and (e, f) non-convective rainfall (mm d −1 ). The south Loess Plateau (SLP) and east Loess Plateau (ELP) regions are defined in Fig. 2. simulated a larger UDROFF decrease due to further revegetation (Fig. 8d), especially over SLP and ELP where the regional mean UDROFF decreases by 0.38 mm d −1 (−54 %) and 0.63 (−42 %), respectively (Fig. 6c, d). These UDROFF decreases are statistically significant at a 95 % confidence level for both SLP and ELP. Moreover, the upper quartile of UDROFF changes among the 20 members systematically shift below 0 mm d −1 for both the SLP and ELP. These results indicate a larger chance of a UDROFF decrease if revegetation is continued over the SLP and ELP. Moreover, the spatial change in UDROFF is consistent with that of the net budget of RAIN and ET (RAIN-ET) for both LC 2015 -LC 2001 and LC futr -LC 2001 (Fig. 8e, f), suggesting that the UDROFF change can be mostly explained by the change in RAIN-ET. We also note some UDROFF changes in adjacent regions of the Loess Plateau (Fig. 8c, d) associated with RAIN changes (Fig. 7a, b).
Compared with the UDROFF change, the surface runoff (SUROFF) change is mostly small for both LC -LC 2001 and LC futr -LC 2001 (Fig. 8a, b). However, the relative change of SUROFF is considerable, especially for LC futr -LC 2001 in which SUROFF decreased by 21 % for SLP and 14 % for ELP (Fig. 6c, d). We also find that the upper quartile of the SUROFF change systematically shifts below 0 mm d −1 , al- though the SUROFF change is not statistically significant for LC futr -LC 2001 .

Impacts on soil moisture
In addition to the decline in runoff, the soil moisture (SMOIS) of each layer is significantly reduced over the Loess Plateau for LC -LC 2001 (Fig. 9a, c, e, g) with larger decreases in the middle two layers. The regional mean SMOIS for the SLP decreases by 0.02 m m −3 (−8 %) and 0.03 m m −3 (−12 %) for the second and third layers, respectively (Fig. 6a). WRF simulated further falls in soil moisture following further revegetation, with a larger impact on deeper soil layer moisture (Fig. 9b, d, f, h). For example, the decrease in the regional mean soil moisture of the bottom layer for the SLP varies from −0.01 (or −5 %) in LC 2015 -LC 2001 (Fig. 6a) to −0.04 (or −17 %) in LC futr -LC 2001 (Fig. 6c). Similar to the UDROFF change, the spatial change in SMOIS for each layer is consistent with that of RAIN-ET for both LC 2015 -LC 2001 and LC futr -LC 2001 (Fig. 8e, f).

Robust identification of rainfall change
We found a large variability in changes in RAIN among the 20 members over the SLP and ELP for both LC 2015 -LC 2001 and LC futr -LC 2001 . We next examine whether these can be attributed to revegetation. We first show the RAIN change in individual members for LC -LC 2001 (Fig. 10). Large variability in RAIN changes among the 20 members occurs throughout the study region. Even the increase in RAIN over the northeast Loess Plateau (Fig. 7a), which is available by comparing multiyear mean RAIN between LC 2001 and LC 2015 , is not consistent for every year. As for the northeast Loess Plateau, the RAIN shows an increase in 8 years (1997, 2001, 2003, 2004, 2007, 2010, 2012, and 2015), a decrease in 5 years (1996, 1999, 2006, 2009, and 2014), and negligible change in the other 7 years. This results in a net increase in RAIN over the 20 years, but a different selection of years could show an overall decrease (the result is similar for LC futr -LC 2001 , not shown). Similarly, other statistically significant RAIN changes occur in the study region (e.g., decreased RAIN to the southwest Loess Plateau shown in Fig. 7a), but these are not consistent across the 20 years. As mentioned earlier, this large variability in RAIN changes among the 20 members is possibly attributed to different boundary conditions (background climate); we next examine whether this is true over the Loess Plateau.
We note that the pattern of RAIN change in 2001 is very similar to the multiyear averaged pattern, although with a larger magnitude (Figs. 7a, 10f). The RAIN increase of the northeast Loess Plateau in only 2001 explains about 30 % of the multiyear mean RAIN increase in the same region. Therefore, we show the RAIN change in each realization for LCENS 2015 -LCENS 2001 in Fig. 11. These 11 ensemble members share the same boundary conditions with small differences in initial conditions. In contrast with the increased RAIN obtained from setting the initial date to 1 May (Fig. 10f), the RAIN changes are modified by an advance of 1 to 10 d in initial conditions. For example, WRF cannot simulate the increased RAIN over the northeast Loess Plateau when using an initial date of 22, 25, 27, or 30 April, highlighting that the RAIN change is very sensitive to the initial conditions. Thus, the RAIN increase in 2001 with an initial date of 1 May is likely associated with internal variability rather than revegetation. In another words, the RAIN change due to revegetation is negligible relative to the RAIN change induced by internal variability. Thus, we conclude that the multiyear averaged RAIN increase over northeast Loess Plateau for LC -LC 2001 (Fig. 7a) cannot be robustly linked to revegetation.

How many members do we need to get a robust signal?
Model internal variability is inevitable when we use models to investigate the impact of land cover change on climate. The model internal variability can be minimized as the number of individual realizations is increased to form a larger sample to calculate any average. Therefore, we examine the relationship between the RAIN change and the number of realization members (Fig. 12)

Discussion
Following the launch of the GFGP by China in the late 1990s, the Loess Plateau has shown a significant greening trend, although with simultaneous concerns about water security for agriculture and other human activities. We investigated the impact of revegetation on the hydrology of the Loess Plateau since the launch of the GFGP using WRF. Simulations show that the revegetation of the plateau is associated with a decrease in runoff and soil moisture as a consequence of higher evapotranspiration and little feedback from rainfall. Our results on changes of evapotranspiration, soil moisture, and runoff are broadly consistent with both field (Jia et al., 2017;Jian et al., 2015;Jin et al., 2011) and satellite (Feng et al., 2017;Li et al., 2016;Xiao, 2014) observations. For example, the spatial pattern of our simulated soil moisture decline in the growing season is similar to observations from the Advanced Microwave Scanning Radiometer on the Earth Observing System by the Japanese Aerospace Exploration Agency (Feng et al., 2017). Although the increased evapotranspiration due to revegetation of the Loess Plateau has been examined before (e.g., Cao et al., 2017Cao et al., , 2019Li et al., 2018;Lv et al., 2019b), the reduction in runoff and soil moisture in response to revegetation of the Loess Plateau, which is consistent with observations, has rarely been reported in modeling results to date. Moreover, our simulated weak response of rainfall to the revegetation of the Loess Plateau, which is hard to determine from observations, is useful in assessing the hydrometeorology of this region. We also investigated the potential future impact on the hydrology of the Loess Plateau if revegetation was continued, which has not been assessed before but is important for both scientific communities and policymakers. WRF suggests that further revegetation would exacerbate soil moisture and runoff declines with particularly large effects on the underground runoff and soil moisture in deeper layers. Our simulations suggested that the potential revegetation that could still be achieved would have larger consequences than those simulated since the launch of the GFGP. Our results provide useful advances in our understanding of the impact of further revegetation on the Loess Plateau. For example, both Feng et al. (2016) and S. L.  estimated the current vegetation over the Loess Plateau is approaching or may have exceeded the threshold of ecological equilibrium. They omitted the potential response of rainfall to further revegetation over the Loess Plateau when predicting future thresholds S. L. Zhang et al., 2018). Our result demonstrate that there is almost no feedback of rainfall associated with further revegetation, supporting the approach of Feng et al. (2016) and S. L.  in this specific region. That said, our approach does not attempt to incorporate changes in climate over the Loess Plateau; thus, the viability of large-scale reforestation in this region is not something that we attempted to assess.
We focused on the response of rainfall to revegetation over the Loess Plateau, which is probably the most uncertain of the hydrological components. WRF shows little response of rainfall to revegetation since the launch of the GFGP, which contradicts earlier results (Cao et al., 2017,  2019; Li et al., 2018;Lv et al., 2019b). Moreover, the rainfall is weakly affected by further revegetation despite a large increase in evapotranspiration. We also demonstrate that the rainfall change is strongly affected by internal variability and a large number of realizations are required before any impact of revegetation on rainfall might be robustly identified. We suggest that some previous studies (Cao et al., 2017(Cao et al., , 2019Lv et al., 2019b) based on model simulations may have exaggerated the impact of revegetation on rainfall over the Loess Plateau due to the lack of sufficient realizations. For example, Cao et al. (2017Cao et al. ( , 2019 and Lv et al. (2019b) used the same WRF to perform only three-or five-member simulations, and concluded a significant change in rainfall caused by revegetation over the Loess Plateau. More interestingly, Cao et al. (2017Cao et al. ( , 2019 obtained different conclusions on the rainfall change over the Loess Plateau with same WRF model. They used a broadly similar experimental design but a different spatial resolution (30 and 10 km, respectively) and simulations from 2001 to 2002 with three ensembles and consecutive simulation from 2000 to 2004, respectively. We could also demonstrate large changes in rainfall over the plateau if we chose three to five members, but we could demonstrate either large increases or large decreases in three-to fivemember averages. Returning to Fig. 6, ET shows a highly consistent increase in the response to revegetation among the 20 years, suggesting that ET change is robustly linked to revegetation. Although changes in runoff and soil moisture also show large variability among the 20 years, the distribution of the runoff and soil moisture changes are negatively biased. More importantly, the distribution of the runoff and soil moisture changes systematically shift towards negative values. This suggests that runoff and soil moisture changes are very likely linked to revegetation. The large variability in runoff or soil moisture changes is induced by the large variability of rainfall. Given the tight linkage between rainfall and runoff or soil moisture, the changes in runoff or soil moisture tend to be mistakenly represented if the rainfall change is not robustly examined, and this requires internal model variability to be thoroughly addressed. Our studies are also subject to some caveats. First, observations of soil moisture declines associated with revegetation can be alleviated once trees mature (Jia et al., 2017;Jin et al., 2011). Our simulations only capture an initial decline in runoff and soil moisture linked to the higher evapotranspiration, and we note that the impact of revegetation on the long-time trend (25-50 years) would be valuable. Second, we used current boundary conditions (1996)(1997)(1998)(1999)(2000)(2001)(2002)(2003)(2004)(2005)(2006)(2007)(2008)(2009)(2010)(2011)(2012)(2013)(2014)(2015) for WRF to predict the impact of further revegetation on the hydrology, which means that the boundary conditions do not change in the future in response to climate change. This suggests that we might underestimate the impact of further revegetation in the future if the future climate of the Loess Plateau suffers from large changes in response to global warming. Third, uncertainties exist in the current land surface model used to represent the response of vegetation to climate change in future. While using satellite observations to construct the land surface biogeophysical parameters helps overcome some land surface parameter limitations, this approach is obviously lim-ited looking forward in terms of the status of future vegetation. Furthermore, we note that our results are likely model dependent, as we only used one model. Although we performed relatively high-resolution simulations (10 km for the nested domain), the cumulus convection scheme remains necessary and is a further potential source of uncertainty. These factors account for the discrepancy between our result and another model-based study .  found a positive rainfall feedback to greening and, consequently, small changes in runoff and soil moisture over north China using a global climate model. In contrast, we demonstrate the rainfall change is too small to compensate for the strongly enhanced evapotranspiration, causing a reduction of runoff and soil moisture in response to revegetation over the Loess Plateau. A large ensemble of models, each with a reasonable number of realizations, is needed to build a model-independent assessment of the impact of revegetation; however, this is clearly beyond the scope of this study. Last, we investigated the impact of revegetation  Fig. 2. For a given number of realizations, the rainfall is averaged over these members. The gray area denotes the range of rainfall changes from all possible combinations of a given number of members. The red dashed line denotes the 5th and 95th percentile of the rainfall changes from all possible combinations of a given number of members. or greening, rather than GFGP, on the hydrology of the Loess Plateau. Directly linking our results to the impact of GFGP on the hydrology of the Loess Plateau should be avoided.
Overall, our results highlight how revegetation of the Loess Plateau led to increased evapotranspiration and how, as a consequence, the runoff and soil moisture declined. This is consistent with the understanding of land surface processes and how they respond to land cover change (Bonan, 2008). Critical in this impact of revegetation on the hydrology is what happens to rainfall. If the higher evapotranspiration increases rainfall, revegetation has the potential to increase soil moisture and runoff. It is very likely this would be the consequences in some regions, such as Amazonia (Lawrence and Vandecar, 2015;Perugini et al., 2017;Spracklen et al., 2018) and Sahel (Kemena et al., 2018;Xue and Shukla, 1996;Yosef et al., 2018). However, over the Loess Plateau we find no such result; thus, the higher evapotranspiration simply leads to lower soil moisture and runoff. Additionally, Bargues Tobella et al. (2014) reported a positive impact of trees on soil hydraulic properties influencing groundwater recharge when termite mounds are taken into account in Africa. However, termite mounds are rare over the Loess Plateau suggesting that this positive impact of trees is unlikely to occur. An implication of this result is that further revegetation, which requires water to be sustained, may not be viable. We also recognize that afforestation can help to sequester carbon, mitigate warming, and alleviate soil erosion. Therefore, if and how to implement further revegetation should be cautiously determined with the pros and cons of afforestation being carefully weighted for the Loess Plateau.

Conclusions
We evaluated how the growing season hydrology of the Loess Plateau has been impacted by revegetation since the launch of the "Grain for Green Program" and by further revegetation in the future using the WRF model. We used satellite observations to describe key biophysical parameters including decreased albedo and increased leaf area index and the fraction of photosynthetically active radiation. The observed greening trend increased evapotranspiration, but because the impact on rainfall was negligible the underground runoff and soil moisture both decreased. Further future revegetation enhanced evapotranspiration but still had little impact on rainfall. Thus, overall, revegetation over the Loess Plateau leads to higher evapotranspiration and, as a consequence, lower water availability for agriculture or other human demands. Considering the negative impact of revegetation on runoff and soil moisture and the lack of benefits on rainfall, we caution that further revegetation may threaten local water security over the Loess Plateau.