InVEST Model-Based Estimation of Water Yield in North China and Its Sensitivities to Climate Variables

Abstract

A revegetation program in North China could potentially increase carbon sequestration and mitigate climate change. However, the responses of water yield ecosystem services to climate factors are still unclear among different vegetation types, which is critically important to select appropriate species for revegetation. Based on the Integrated Valuation of Ecosystem Services and Tradeoffs (InVEST) model, we estimated the temporal variations and associated factors in water yield ecosystem services in North China. The result showed that the InVEST model performed well in water yield estimation (R² = 0.93), and thus can be successfully applied across the study area. The total water yield across North China is 6.19 × 10¹⁰ m³/year, with a mean water yield (MWY) of 47.15 mm/year. A large spatial difference in the MWY was found, which is strongly related to temperature, precipitation, and land use types. The responses of the MWY to mean annual precipitation (MAP) are closely tied to temperature conditions in forests and grasslands. The sensitivities of the MWY to climate variables indicated that temperature fluctuation had a positive influence on the forest MWY in humid regions, and the influence of precipitation on grassland water yield was enhanced in warmer regions. We suggest shrub and grass would be more suitable revegetation programs to improve water yield capacity, and that climate warming might increase the water yield of forests and grasslands in humid regions in North China.

1. Introduction

Revegetation has taken place in degraded ecosystems in North China during previous decades, and subsequently resulted in positive ecosystem service sequences, such as increasing carbon sequestration [1], protecting soil conservation and biodiversity, and influencing hydrological processes [2]. With the problem of water scarcity, hydrological ecosystem services are becoming a hot issue in ecosystem management in North China. Water yield, as a key component of the hydrological ecosystem service, represents a key role in ecosystem management and hydrological balance [3]. However, a large uncertainty exists in water yield estimation, especially as large temporal–spatial variations in water yield can be expected under climate change conditions [4]. Hence, the accurate estimation and quantification of influential factors on water yield are urgently needed to make proper ecosystem service decisions, like revegetation strategies, and to secure the water demand of the socio-economic system [5].

Water yield is mainly controlled by precipitation and evapotranspiration (ET), and human-induced land use changes also indirectly affect the water yield [6]. There are many cases of the estimation of the impacts on water yield in regions where both climatic factors and land use changes have occurred, and the combined impacts of climate and land use change exhibit positive or negative feedbacks. For instance, many studies have reported an accordant positive correlation between precipitation and water yield [7], while land use changes (e.g., revegetation or afforestation) showed different impacts on water yield [8,9]. It should be noted that the sensitivities of water yield to climate variation among different land use types needs to be considered.

Previous studies of the hydrological consequences of revegetation reported that both canopy interception and ET increased, while the runoff decreased, leading to declining stream flow and soil moisture [10]. Revegetation reduced the water yield from 5% to 8% during the last three decades in North China’s semi-arid Loess Plateau [11]. By contract, a recent study found that afforestation might have positive influences on the hydrological cycle by facilitating the large-scale transport of water vapor and increasing precipitation at regional to global scales [9]. This is may be because the natural forests, with high transpiration fluxes, enhanced the ascending air motion over the forest and “sucked in” moist air that was transported from the remote ocean, so that the runoff losses from optimally moistened soil were fully compensated by increased precipitation [12]. Positive hydrological feedback associated with revegetation was also found based on land–atmosphere global climate models and forest inventory data in temperate humid climate zones [13]. Hence, considerable uncertainty remains in water yield estimation and its response to climate change among different vegetation types [14]. Therefore, it is important to estimate the sensitivity of water yield to climate change across different vegetation types.

Distributed physical hydrological models, based on coupling hydrologic processes with biological and geochemical processes, have been applied to forecast the water resource dynamics at the basin scale [15], such as the Soil and Water Assessment Tool (SWAT), Variable Infiltration Capacity (VIC), and Integrated Value Environmental Services Tradeoff (InVEST). Despite providing extensive outputs, the SWAT model was developed to estimate the variables at a watershed scale rather than a pixel scale. The VIC model is suitable for research at a large scale and the outputs vary by time step (form hourly to daily) and spatial scale (including gridded and watershed results). Compared with other models, the InVEST model is based on the Budyko curve [16] and annual precipitation and runs with relatively low data requirements, providing effective means to estimate an ecosystem service, like water yield, at a high spatial resolution and at different scales [17,18]. At the basin scale, the InVEST model had been applied to several studies to explore general patterns and variations in ecosystem services caused by land use changes or climate change impacts, and performed well with various geographic and climatic characteristics [19].

There have been many studies focused on the assessment of water yield sensitivity to changes in precipitation and land use types. Some studies assessed the sensitivity to climate change by modifying temperature and precipitation to certain magnitudes as model inputs [5,20,21]. For example, Redhead and Stratford validated the InVEST water yield model at several independent basins across the UK, and investigated model sensitivity by subjectively adjusting the amplitudes of the input data [21]. However, the vegetation changes and vegetation–climate feedback were omitted in these studies [9]. Other studies analyzed the variations of climate and its contributions to water yield based on Budyko’s framework approach using a theoretical model—the first order Taylor expansion of the Choudhury–Yang equation [22,23]—to quantify the direct contributions of precipitation, potential evapotranspiration(PET), and land surface characteristics to runoff changes. Earlier studies, based on theoretical models, estimated the change in annual runoff based on drought metrics characterized only by precipitation amount and PET [24,25]. Recent studies focused on the effect of vegetation change on regional water balance [26] and indicated that the sensitivity of runoff to change in the landscape parameter n (e.g., land use types) is higher in dryland regions [27], and the landscape parameter n is more sensitive to vegetation change in relatively dry catchments [28]. However, this approach has its limitations, including being incapable of quantifying past and future drivers of change or providing uncertainty estimations of water yield. Compared to these theoretical models that mainly estimated the water yield based on climate change and land cover change scenarios, the statistical models might be simple and effective enough to explore the sensitivity of the phenomenon to a few influential factors. Besides that, the spatial sensitivity of water yield capacity also needs to be explored, especially in the arid and semi-arid regions of North China in the grip of dramatic climate change. Therefore, the performance of the InVEST model at a large spatial scale with different climate types still needs to be explored, especially in arid and semi-arid regions.

In this study, we therefore used the InVEST water yield ecosystem model to estimate the spatial variation of water yield capacity in North China, which has different climate and ecosystem types and is highly sensitive to climate change [4]. The main objectives of this study are: (1) to validate the reliability of the InVEST model at a large scale in arid and semi-arid regions of North China, (2) to quantify the spatial patterns of the mean water yield (MWY), and (3) to explore the sensitivity of water yield capacity to climate variables among different land use types.

2. Materials and Methods 2.1. Study Region

The study region covers about one-sixth of the total land area of China, with a total area of 1.53 × 10⁶ km², which extends from 34.58° N to 53.36° N and 97.01° E to 125.93° E. The study region includes five large basins: Songhua River basin (SR), Liao River basin (LR), Yellow River basin (YR), Northwest River basin (NWR), and Hai River basin (HR) (Figure 1). The mean annual temperature and precipitation exhibit significant spatial heterogeneity, ranging from 2 °C to 1 °C and 0 to 750 mm, respectively. The study regions are characterized by water shortages, while the climate type in the southeast is warm, temperate, and humid with monsoons. The land cover includes desert, semi-arid steppe, irrigation agriculture, shrub, and forest. Since the 1990s, North China became a key region for China’s revegetation programs, which were designed to hold back the expansion of deserts, such as the Hobq desert, Tengger desert, and Ulan Buh desert. The influences of revegetation on hydrological cycles in this region are very controversial.

2.2. Water Yield Model

The water yield module of the InVEST model runs on a gridded map and estimates the quantity of water in each pixel based on the principle of water balance. The annual water yield of each pixel is determined as follows:

_Yxj=(1−_AETxjPx)×_Px

where_Yxjand_AETxjare the annual water yield and actual evapotranspiration for pixel x on the landscape j, respectively._Pxis the annual precipitation on pixel x.

The relationship between_AETxjand_Px is based on methodology developed by Budyko [16], and later adapted by Fu [29] and Zhang et al. [30]._AETxj/_PETxj, is estimated in a spatially explicit way on pixel x, that is

_AETxjPx=1+_PETxjPx−^{[1+(_PETxjPx)ω]1/ω}

where_PETxjis defined as the potential evapotranspiration for pixel x on landscape j, and ω is a nonphysical parameter that characterizes the natural climate–soil properties._PETxjcan be calculated as follows:

_PETx=_Kc(_ℓx)×E_T0(x)

where_Kc(_ℓx)is the vegetation evapotranspiration coefficient associated with the land use type (_ℓx) in each pixel,_Kc is extracted from the InVEST user’s guide [31], the detail is shown as follows (Table 1):

_ET0(x) is the reference evapotranspiration for pixel x, which was derived using the Hargreaves equation [32]. The Hargreaves equation is given as follows:

_ET0(x)=0.0023×Ra×[(_Tmax+_Tmin)/2+17.8]×^{(_Tmax+_Tmin)1/2}

where Ra is extraterrestrial radiation (MJ m⁻² d⁻¹);_Tmaxand_Tminare the mean maximum and minimum temperature in °C, respectively.

ω(x) characterizes the natural climatic–soil properties, where AWC is the volumetric plant-available water content and N is the number of rainfall events per year. ω(x) is calculated as follows:

ω(x)=Z(AWC(x)P(x))+1.25

where AWC(x) is estimated as the product of plant-available water capacity (PAWC) on pixel x, the minimum of root restricting layer depth and vegetation rooting depth [31]._Pxis the annual precipitation on pixel x. Z is an empirical factor that captures the local precipitation and additional hydrogeological characteristics, which are positively correlated with the number of rainfall events (N) per year (Z = 0.2 × N).

2.3. Data Sources and Processing

The data input to run the water yield module include meteorological data, land use, soil depth, plant-available water fraction, watersheds, and the biophysical table. For this study, the meteorological data and land use data (Figure S1) encompassed the years 2000, 2005, 2010, and 2015. For each year, the distribution of the water yield was estimated based on annual meteorological data and land use data, and the mean water yield was obtained as the mean of the water yield of the studied years (2000, 2005, 2010, and 2015). All of the data were resampled at a spatial resolution of 1 km and projected coordinate system of the World Geodetic System 84 (WGS84).

The meteorological data, including temperature, precipitation, and radiation data, were acquired from the China National Meteorological Information Center (http://data.cma.cn/), and spatially interpolated into 1 km × 1 km grid pixels according to the geographically weighted regression Kriging algorithm (Figures S2–S4). The annual potential evapotranspiration was calculated using the Hargreaves equation for each station, then spatially interpolated into 1 km × 1 km grid pixels (Figure S5).

The land use layers were obtained from the Resource and Environment Data Cloud Platform (RESDC) (http://www.resdc.cn), which are classified into 25 kinds of land use categories. For this study, the land use data of each year was aggregated into seven land use types, including forest, shrub, grassland, cropland, water, built-up area, and bare land.

The plant-available water capacity (PAWC) is defined as the difference between the wilting point and field capacity, which can be estimated based on the physical and chemical properties of soil [33]. The equation is given as follows:

PAWC=54.509−0.132×SAND−0.003×^SAND2−0.055×SILT−0.006×^SILT2−0.738×CLAY+0.007×^CLAY2−2.688×SOC+0.501×^SOC2

where the SAND, SILT, CLAY, and SOC are the proportion of sand, silt, clay, and organic matters in the soil, respectively. The data were converted to a fraction from 0 to 1 in accordance with the input requirements (Figure S6). The layers of soil depth and soil texture were obtained from the Harmonized World Soil Database [34].

Parameter Z, as a constant, represents the precipitation characteristics with a value from 1 to 30 and it is larger when the rainfall events are more frequent. The constant can be estimated as 0.2 × N, where N is the number of rainfall events (N > 1 mm) per year [21,35]. In this study, N was estimated based on daily precipitation from observation stations in study area.

2.4. Sensitivity Analysis To investigate the relationship between water yield and climate factors, we extracted the pixels of five land use types, including cropland (CL), forest (FO), shrub (SH), grassland (GL), and bare land (BL). All the pixels of each type were classified into 20 mm mean annual precipitation (MAP) and 0.5 °C mean annual temperature (MAT) bins, and the average water yield of each bin was calculated in a MAT–MAP space. To further investigate how water yield varies with precipitation and temperature, respectively, linear regressions were performed to calculate the sensitivity (slope) of the water yield to precipitation (Sp) in each 1 °C temperature bin, and the sensitivity of the water yield to temperature (St) in each 100 mm/year precipitation bin. We compared our result against the observed data using linear regression analysis. The p-value, R² value, and the root mean square error (RMSE) were estimated in the validation. All statistical analyses were conducted using MATLAB R2016a (MathWorks, Natick, MA, USA). 3. Results 3.1. Model Validation

To evaluate the performance of the InVEST model, we collected the measured actual evapotranspiration (AET) from a dataset (actual evapotranspiration and water use efficiency of typical terrestrial ecosystems in China) [36], and observed the yearly data in five independent basins from which the major rivers originate in the same basin (Figure 2A). The six sites of the dataset include one forest site, two dry cropland sites, and three grassland sites. The basins include the Inner Mongolia plateau inward flowing river basin (IMPIR), the Yellow River inward flowing river basin (YIR), the north part of the Hai River basin (NHR), the Luan River basin (LR) and the western Liao River basin (WLR). The yearly data of water resources, including the period of 2000, 2005, 2010, and 2015, were obtained from the statistical data of local water resources bulletin (http://slt.nmg.gov.cn/index.html, http://www.hwcc.gov.cn/). The data were converted to the water yield with the unit of mm/year (referred to the observed MWY), based on the size of each basin.

As shown in Figure 2B, the result shows that strong linear relationships exist between the estimated MWY and the observed MWY (R² = 0.93, p < 0.01, RMSE = 1.08). The estimated AET also exhibited strong linear relationships with the observed AET (R² = 0.95, p < 0.01, RMSE = 34). The result indicated that the InVEST model can be used to estimate the annual MWY at the regional scale.

3.2. The Water Yield in North China

Based on the InVEST model, we estimated the MWY and AET of the studied years (2000, 2005, 2010, and 2015) (Figures S7 and S8). We found that the total water yield of North China in 2015 was 6.19 × 10¹⁰ m³/year, with a mean water yield (MWY) of 47.15 mm. The water yield of the Hai River basin is higher than others, which accounts for 44% of the total water yield in North China. The water yield of the Songhua River is second and accounts for 28% of the total water yield. The high water yield of the Hai River basin can be mainly attributed to its higher MWY, which is 101.28 mm. Despite occupying a third of the total area, the water yield of the Northwest River basin is lower than the others, with an average MWY of 8.22 mm, significantly lower than the Songhua River basin, Liao River basin, and Hai River basin (p < 0.01) (Table 2).

As shown in Figure 3A, the spatial patterns of the MWY showed a strong heterogeneity across North China. The MWY of the northeast and mid-south areas of the study region was relatively higher than others. The highest MWY was found in construction land, with values higher than 300 mm/year, because of the impermeable underlying surface and strong evaporation [14]. In contrast, for more than seventy percent of the pixels, the MWY is lower than 50 mm/year (Figure 3B), and are mainly distributed across the mid-north and west areas of the study region, west of the Songhua River basin and the middle of the Liao River basin.

3.3. Difference in Water Yield and Associated Climate Variables among Land Use Types

To analyze the difference in water yield among vegetation types, we extracted the pixels including cropland (CL), shrub (SH), forest (FO), grassland (GL), and bare land (BA) and the distribution of land use types is showed in Figure 4A. Cropland was mainly distributed in the central south of the study region, while grassland was widespread across the central north of the study region. Forest was mainly distributed in the northeast of the region, and bare land was the main land use type in the west of the study region. Figure 4B illustrates that cropland, forest, and shrub were mainly distributed in humid regions where the MAP was larger than 430 mm/year, while the MAP of bare land was generally lower than 200 mm. Grassland was widely distributed in humid, arid, and semi-arid regions. Across all land use types, the MWY of cropland, forest, and shrub was higher than for grassland and bare land (Figure 4D). In detail, cropland produced the largest MWY (58.64 mm/year) and the MWY of shrub (48.58 mm/year) and forest (46.31 mm/year) was relative small. The MWY of grassland (28.11 mm/year) was lower than for the former, but higher than for BA (15.84 mm/year).

Similar variations were also found in the comparison of the mean actual evapotranspiration (AET) (Figure 4C), The AET of forest, shrub, and cropland was higher than for grassland and bare land. In detail, the AET of most forests was higher than 450 mm/year, and the AET of cropland and shrub was mainly higher than 400 mm/year. Meanwhile, the AET of grassland was also higher than for BA, which was higher than 220 mm/year. The AET of bare land was the lowest, which was normally lower than 200 mm/year.

3.4. Spatial Relationships between the MWY and Climatic Factors

The distribution of the MWY among different climatic variables is shown in Figure 5. Besides bare land, each of the other four land use types were distributed in regions with MATs ranging from −2.5 °C to 15 °C and MAP ranging from 200 mm/year to 680 mm/year. Bare land was mainly distributed in the west of the study region, with a MAP ranging from 50 mm/year to 550 mm/year. The higher MWY levels (> 100 mm/year) mainly occurred in regions with a MAP higher than 500 mm/year and a MAT higher than 5 °C. The relatively low MWY levels (< 25 mm/year) were mainly found in arid and semi-arid regions with a MAP less than 400 mm/year.

To further investigate the sensitivity of the MWY to climatic factors, we calculated the temperature sensitivity (St) of the MWY along the MAT gradient and precipitation sensitivity (Sp) along the MAP gradient for each land use type. As shown in Figure 5, Sp was mainly positive along with the MAT gradient of all the land use types, while there was no similar trend in St along the MAP gradient. Notably, the St of forests was close to 0 with a MAP lower than 500 mm/year, but it increased with a MAP higher than 500 mm/year. Meanwhile, an increase in Sp with a MAT higher than −2.5 °C was found in grassland.

4. Discussion 4.1. Evaluation of the Water Yield Model in the InVEST Model

The accurate assessment of water yield capacity on a large spatial scale is fundamental in ecosystem service research and water resource management, especially in the context of global climate change. Previous studies in China mainly assessed the performance of the InVEST water yield model in watersheds or sub-watersheds within a single river basin [14,37]. In this study, we applied the InVEST model to North China across five basins with substantial differences in climate and land type, compared our results against the annual water yield for five river basins in 2000, 2005, 2010, and 2015 (which was obtained from the local water resources bulletin), and measured the actual evapotranspiration (AET) from a dataset of typical terrestrial ecosystems. We also found that the estimated MWY was comparable to previous studies. Based on a process-based ecosystem model, Liu et al. found that the MWY in North China fell in a range of 0–200 mm/year during the period 2000–2014 [38]. The spatial pattern of the MWY was also consistent with ours, MWYs larger than 100 mm/year were mainly located in the central south of the HR basin and the east of the SR basin.

Despite the good performance of the InVEST model, its acknowledged simplification and limitation might affect the accuracy of the modeled water yield values, including its inability to account for inter- or intra-annual variation. Negative influences by issues like irrigation, mining [39,40], and hydropower supplement on water resource variation were found in many studies, which was not considered in the water yield model of the InVEST model [31]. A relatively poor performance of the models was found in the WLR and LR basins with a large amount of irrigated farmland, which indicated that water consumption and management might increase the uncertainty of the water yield estimation [41]. Meanwhile, the InVEST model is incapable of accounting for groundwater and water resource infrastructures that redistribute water resources, which might be important to describe the water resources status of North China [42]. Compared to the InVEST model, VIC and SWAT are more suitable for providing multiple and reliable estimations depending on available input data and the objective of the modeling exercise. As an important input parameter of the model, the reference evapotranspiration was calculated using the Hargreaves method because of the sparse climate data. Nevertheless, the Penman–Monteith or Thornthwaite methods should be applied where accurate weather data collection is available when conducting higher temporal solution estimations [43].

4.2. Water Yield and Vegetation Types

For fhe comparison of precipitation, the AET and MWY showed that the characteristics of water yield among different levels of vegetation cover were closely related to precipitation (Figure 4). Bare land was mainly distributed in arid and semi-arid regions with a MAP lower than 200 mm, so the water yield is significantly lower than for the others. However, attention should be paid to understanding the difference of the MWY among vegetation types distributed in similar rainfall conditions, which is critically important to select appropriate species for revegetation. It is well known that the MWY reflects the combined effects of the MAP and MAET, hence the regions with moderate to high precipitation do not necessarily have a higher MWY. Previous studies indicated that land use change had the largest absolute impact on runoff in humid areas [44,45]. In line with this argument, we found that, compared to forest, the MWY of shrub is slightly higher but not statistically significant, with a higher MAP but a lower AET. This may be because of a lower capacity for water loss that is associated with a lower leaf area [45,46] and lower canopy interception [47]. We also found that the MWY of grassland is not significantly lower than that of cropland, forest, and shrub, even though the MAP and AET are significantly lower. This may be because grassland is mainly distributed in arid and semi-arid regions, where the evapotranspiration is generally lower than in forest and shrub [47].

Recent revegetation in North China explains the spatial pattern of vegetation greening [48], including reforesting and regrassing in cropland and bare land. Many studies indicated that afforestation might increase canopy interception and the loss of water by transpiration, and hence reduce stream runoff and water yield [11,49], especially in arid and semi-arid regions [50]. Our result showed that choosing the appropriate vegetation type to revegetate might improve the water yield capacity of an ecosystem. Compared with afforestation, planting shrubs and grasslands might be more suitable to increase regional water yields in arid and semi-arid regions.

4.3. Sensitivity of the MWY to Temperature and Precipitation

Many studies have assessed the response of water yield to climate change and found a higher sensitivity of water yield to precipitation comparing to temperature [5,7], but the interactive effect between precipitation and temperature on water yield is still unclear. Many studies assessed the sensitivity of the MWY to climate factors by modifying the temperature and precipitation to certain magnitudes as model inputs [20,21,27,28,51], but the interactive effect between precipitation and temperature on water yield is still unclear. The former approach to the sensitivity analysis might be not sufficient due to a lack of practice guidance because the water yield capacity could respond to the adaptation of the ecosystem to the local climate [9]. The latter provided a quantitative assessment of the impact of drivers on water resource changes in the Budyko’s framework and it also had the inherent limitations of theoretical estimations, including a lack of uncertainty estimates, and more inputs were needed when quantifying past and future drivers of change [27].

Using a uniform method to explore the differences in hydrologic changes across environmental gradients provides an alternative approach [5,52]. Meanwhile, the water yield model is based on a gridded map, which allows us to take the heterogeneity of physiological characteristics into consideration, such as evapotranspiration characteristics and morphological characteristics, which are the key driving factors of the MWY and AET variation among different vegetation types. The validation of the water resource and actual evapotranspiration also provides valuable feedback and helps to improve the accuracy of the estimation of the annual datasets, so the sensitivity analysis could be conducted based on the estimation in the absence of water yield observation data.

In line with previous studies, the gradient distributions of the MWY are consistent with that of the MAP [53], i.e., the MWY increased with the MAP in almost all MAT ranges (Figure 5). Interestingly, we found that the sensitivities of the MWY to climate variables were different. The sensitivities of the MWY to the MAT (St) in forests was close to 0 with a MAP less than 500 mm, while it increased with a MAP larger than 500 mm. The variation of St indicated that the temperature fluctuation had less influence on the forest MWY in relatively arid regions but a positive influence in humid regions, which might be because of the higher leaf cover of forests in warmer zones, which decreased the transpiration and evaporation resulting from subcanopy radiation [54]. We also found the Sp in grassland increased with the MAT, which showed that the influence of precipitation on water yield was enhanced in warmer zones. That may be because the evapotranspiration of grassland is limited to the transpiration capacity of grass. We therefore suggest that climate warming might increase the water yield of forests and grassland in North China, especially in humid regions.

5. Conclusions In this study, the water yield ecosystem service was estimated using the InVEST model across North China. The results suggest that a relatively simple InVEST water yield model can provide reliable results across a large climatic gradient, as long as the parameters are representative of the spatial and temporal scale concerned. The total water yield of North China was 6.19 × 10¹⁰ m³/year in 2015, with a mean water yield (MWY) of 47.15 mm/year. The total water yield of each basin is influenced by area and the MWY, and the distribution of the MWY is closely related to precipitation. The results also indicate that uncertainties might exist in estimating the water yield of cropland, such as irrigation and cultivation mode. Estimation and sensitivity exploration based on a pixel scale are conducted in this research, offering new insights for understanding the water yield response to ongoing climate change. The spatial analysis of the MWY showed that shrub and grass were widespread in arid and semi-arid areas of North China. In these areas, a higher water yield capacity was observed in many MAP and MAT bins, which suggested that shrub and grass would be more suitable in North China’s revegetation programs. Sensitivity analysis indicated that the water yield capacity of forests and grass land would be enhanced in warmer and humid regions, suggested that climate warming might increase the water yield of forests and grassland in humid regions in North China.

Figure 1. Map of China's physiography regionalization (A) and the studied region (B). China's physiography regionalization includes North China (NC), North East China (NEC), Eastern China (EC), Central South China (CSC), North South China (NSC), and North West China (NWC). The distribution of five river basins is showed in (B), including the Songhua River basin (SR), Liao River basin (LR), Yellow River basin (YR), Northwest River basin (NWR), and Hai River basin (HR).

Figure 2. Model evaluation using the estimated mean water yield (MWY) vs. observed MWY from a water resources bulletin, and the estimated actual evapotranspiration (AET) vs. observed AET from the sites' records. (A) Spatial distribution of the selected basins and sites of the dataset, the areas of the basins vary from 44,095 to 311,853 km2. (B) Evaluation of the estimated MWY versus observed MWY, including data from 2000, 2005, 2010, and 2015. (C) Evaluation of the estimated AET versus observed AET, including data from 2005 and 2010 available in the dataset. Red dashed line indicates a relationship with intercept = zero and slope = 1.

Figure 3. Spatial distribution (A) and frequency histograms (B) of the MWY in North China in 2015. The dotted line in subfigure B refers to the average MWY across the whole region (B).

Figure 4. The spatial patterns of five land use types (A), and the comparison of mean annual precipitation (B), mean actual evapotranspiration (C), and mean water yield (D), for cropland (CL), shrub (SH), forest (FO), grassland (GL), and bare land (BA). In each box, the central line marks the median, the edges of the box correspond to the 25th and 75th percentiles, and the whiskers extend to the range of the data.

Figure 5. Distribution of the MWY in a two-dimensional space with (MAT) and (MAP) binned into intervals of 0.5 °C for the MAT and 20 mm for the MAP. For each figure, the sensitivity of the MWY to temperature (St) along the precipitation gradient is shown under the figure and the sensitivity of the MWY to precipitation (Sp) along the temperature gradient is shown on the right of the figure. The shaped area indicates 95% significance intervals of St and Sp, while the bins have more than 100 grid pixels.

Figure 5. Distribution of the MWY in a two-dimensional space with (MAT) and (MAP) binned into intervals of 0.5 °C for the MAT and 20 mm for the MAP. For each figure, the sensitivity of the MWY to temperature (St) along the precipitation gradient is shown under the figure and the sensitivity of the MWY to precipitation (Sp) along the temperature gradient is shown on the right of the figure. The shaped area indicates 95% significance intervals of St and Sp, while the bins have more than 100 grid pixels.

Land Cover	Kc	Land Cover	Kc
Irrigated cropland	1.2	Stream	1
Field cropland	1	Permanent lentic water	1
Forest closed	1	Flooded/marsh	1
Forest	0.9	Urban non-vegetated	0.1
Natural shrub	0.8	Residential	0.1

Basins	Area (104 km2)	Total Water Yield (109 m3/year)	MWY (mm)
SR	30.93	17.02	55.02
LR	14.00	6.19	44.30
YR	24.78	6.72	26.99
NWR	54.50	4.55	8.22
HR	27.19	27.46	101.28

Supplementary Materials

The following are available online at https://www.mdpi.com/2073-4441/12/6/1692/s1, Figure S1: Spatial distribution of land use types in North China, Figure S2: Spatial distribution of meteorological stations in North China, Figure S3: Spatial distribution of annual precipitation (AP) in North China, Figure S4: Spatial distribution of mean annual temperature (MAT) in North China, Figure S5: Spatial distribution of potential evapotranspiration in North China, Figure S6: Spatial distribution of plant available water capacity in North China, Figure S7: Spatial distribution of mean annual actual evapotranspiration (MAET) in North China, Figure S8: Spatial distribution of mean water yield (MWY) in North China.

Author Contributions

Conceptualization, X.W. and G.Y.; methodology, G.Y. and X.W.; software, G.Y.; validation, G.Y. and X.Z.; formal analysis, G.Y. and X.W.; resources, X.Z. and F.H.; data curation, G.Y. and Y.F.; writing-original draft preparation, G.Y.; writing-review and editing, X.Z., Q.H. and Y.F.; visualization, G.Y.; supervision, X.W.; project administration, X.W. and X.Z.; funding acquisition, X.Z. and F.H. All authors have read and agreed to the published version of the manuscript.

Funding

This research was funded by the National Key Research and Development Program of China (No. 2019YFC0409201), the Chinese National Special Science and Technology Program of Water Pollution Control and Treatment (Grant No. 2017ZX07302004) and Higher Education Discipline Innovation Project (No. B18006).

Conflicts of Interest

The authors declare no conflict of interest.

References

1. Fang, J.; Chen, A.; Peng, C.; Zhao, S.; Ci, L. Changes in forest biomass carbon storage in China between 1949 and 1998. Science 2001, 292, 2320-2322.

2. Chen, Y.; Wang, K.; Lin, Y.; Shi, W.; Song, Y.; He, X. Balancing green and grain trade. Nat. Geosci. 2015, 8, 739-741.

3. Brauman, K.A. Hydrologic ecosystem services: Linking ecohydrologic processes to human well-being in water research and watershed management. Water 2015, 2, 345-358.

4. Piao, S.; Ciais, P.; Huang, Y.; Shen, Z.; Peng, S.; Li, J.; Zhou, L.; Liu, H.; Ma, Y.; Ding, Y.; et al. The impacts of climate change on water resources and agriculture in China. Nature 2010, 467, 43-51.

5. Lu, N.; Sun, G.; Feng, X.M.; Fu, B.J. Water yield responses to climate change and variability across the North-South Transect of Eastern China (NSTEC). J. Hydrol. 2013, 481, 96-105.

6. Bonan, G. Ecological Climatology: Concepts and Applications; Cambridge University Press: Cambridge, UK, 2015.

7. Pessacg, N.; Flaherty, S.; Brandizi, L.; Solman, S.; Pascual, M. Getting water right: A case study in water yield modelling based on precipitation data. Sci. Total Enviorn. 2015, 537, 225-234.

8. Zhang, L.; Cheng, L.; Chiew, F.; Fu, B. Understanding the impacts of climate and landuse change on water yield. Curr. Opin. Environ. Sustain. 2018, 33, 167-174.

9. Li, Y.; Piao, S.; Li, L.Z.X.; Chen, A.; Wang, X.; Ciais, P.; Huang, L.; Lian, X.; Peng, S.; Zeng, Z.; et al. Divergent hydrological response to large-scale afforestation and vegetation greening in China. Sci. Adv. 2018, 4, eaar4182.

10. Brown, A.E.; Zhang, L.; McMahon, T.A.; Western, A.W.; Vertessy, R.A. A review of paired catchment studies for determining changes in water yield resulting from alterations in vegetation. J. Hydrol. 2005, 310, 28-61.

11. Feng, X.; Fu, B.; Piao, S.; Wang, S.; Ciais, P.; Zeng, Z.; Lü, Y.; Zeng, Y.; Li, Y.; Jiang, X.; et al. Revegetation in China's Loess Plateau is approaching sustainable water resource limits. Nat. Clim. Chang. 2016, 6, 1019-1022.

12. Makarieva, A.M.; Gorshkov, V.G. Biotic pump of atmospheric moisture as driver of the hydrological cycle on land. Hydrol. Earth Syst. Sci. 2007, 11, 1013-1033.

13. Liu, Y.; Stanturf, J.; Lu, H. Modeling the Potential of the Northern China Forest Shelterbelt in Improving Hydroclimate Conditions1. J. Am. Water Resour. Assoc. 2008, 44, 1176-1192.

14. Lang, Y.; Song, W.; Deng, X. Projected land use changes impacts on water yields in the karst mountain areas of China. Phys. Chem. Earth Parts A/B/C 2018, 104, 66-75.

15. Legesse, D.; Vallet-Coulomb, C.; Gasse, F. Hydrological response of a catchment to climate and land use changes in Tropical Africa: Case study South Central Ethiopia. J. Hydrol. 2003, 275, 67-85.

16. Budyko, M.I. Climate and Life; Academic: San Diego, CA, USA, 1974.

17. Vigerstol, K.L.; Aukema, J.E. A comparison of tools for modeling freshwater ecosystem services. J. Environ. Manag. 2011, 92, 2403-2409.

18. Dennedy-Frank, P.J.; Muenich, R.L.; Chaubey, I.; Ziv, G. Comparing two tools for ecosystem service assessments regarding water resources decisions. J. Environ. Manag. 2016, 177, 331-340.

19. Lang, Y.; Song, W.; Zhang, Y. Responses of the water-yield ecosystem service to climate and land use change in Sancha River Basin, China. Phys. Chem. Earth Parts A/B/C 2017, 101, 102-111.

20. Kaczmarek, Z. On the Sensitivity of Runoff to Climate Change. IIASA Working Paper. WP-90-05. 1990. Available online: https://core.ac.uk/download/pdf/33894868.pdf (accessed on 11 June 2020).

21. Redhead, J.W.; Stratford, C.; Sharps, K.; Jones, L.; Ziv, G.; Clarke, D.; Oliver, T.H.; Bullock, J.M. Empirical validation of the InVEST water yield ecosystem service model at a national scale. Sci. Total Environ. 2016, 569-570, 1418-1426.

22. Choudhury, B. Evaluation of an empirical equation for annual evaporation using field observations and results from a biophysical model. J. Hydrol. 1999, 216, 99-110.

23. Yang, H.B.; Yang, D.W.; Lei, Z.D.; Sun, F.B. New analytical derivation of the mean annual water-energy balance equation. Water Resour. Res. 2008, 44, 131-139.

24. Teng, J.; Chiew, F.H.S.; Vaze, J.; Marvanek, S.; Kirono, D.G.C. Estimation of Climate Change Impact on Mean Annual Runoff across Continental Australia Using Budyko and Fu Equations and Hydrological Models. J. Hydrometeorol. 2012, 13, 1094-1106.

25. Feng, S.; Fu, Q. Expansion of global drylands under a warming climate. Atmos. Chem. Phys. 2013, 13, 10081-10094.

26. Yang, H.; Piao, S.; Huntingford, C.; Ciais, P.; Li, Y.; Wang, T.; Peng, S.; Yang, Y.; Yang, D.; Chang, J. Changing the retention properties of catchments and their influence on runoff under climate change. Environ. Res. Lett. 2018, 13, 094019.

27. Berghuijs, W.R.; Larsen, J.R.; van Emmerik, T.H.M.; Woods, R.A. A Global Assessment of Runoff Sensitivity to Changes in Precipitation, Potential Evaporation, and Other Factors. Water Resour. Res. 2017, 53, 8475-8486.

28. Zhang, S.; Yang, H.; Yang, D.; Jayawardena, A.W. Quantifying the effect of vegetation change on the regional water balance within the Budyko framework. Geophys. Res. Lett. 2016, 43, 1140-1148.

29. Fu, B.P. On the calculation of the evaporation from land surface. Chin. J. Atmos. Sci. 1981, 5, 23-31.

30. Zhang, L.; Dawes, W.R.; Walker, G.R. Response of mean annual evapotranspiration to vegetation changes at catchment scale. Water Resour. Res. 2001, 37, 701-708.

31. Sharp, R.; Tallis, H.T.; Ricketts, T.; Guerry, A.D.; Wood, S.A.; Chaplin-Kramer, R.; Nelson, E.; Ennaanay, D.; Wolny, S.; Olwero, N.; et al. VEST 3.2.0 User's Guide; The Natural Capital Project: Stanford, CA, USA, 2015.

32. Zhang, L.; Hickel, K.; Dawes, W.R.; Chiew, F.H.S.; Western, A.W.; Briggs, P.R. A rational function approach for estimating mean annual evapotranspiration. Water Resour. Res. 2004, 40.

33. Zhou, W.; Liu, G.; Pan, J.; Feng, X. Distribution of available soil water capacity in China. J. Geogr. Sci. 2005, 15, 3-12.

34. Fischer, G.; Nachtergaele, F.; Prieler, S.; Van Velthuizen, H.T.; Verelst, L.; Wiberg, D. Global Agro-Ecological Zones Assessment for Agriculture (GAEZ 2008); IIASA: Laxenburg, Austria; FAO: Rome, Italy, 2008.

35. Donohue, R.J.; Roderick, M.L.; McVicar, T.R. Roots, storms and soil pores: Incorporating key ecohydrological processes into Budyko's hydrological model. J. Hydrol. 2012, 436-437, 35-50.

36. Han, Z.; Guirui, Y.; Xianjin, Z.; Qiufeng, W.; Leiming, Z.; Zhi, C.; Xiaomin, S.; Honglin, H.; Wen, S.; Yanfen, W.; et al. A Dataset of Actual Evapotranspiration and Water Use Efficiency of Typical Terrestrial Ecosystems in China (2000-2010); China Scientific Data: Beijing, China, 2019.

37. Yu, J.; Yuan, Y.; Nie, Y.; Ma, E.; Li, H.; Geng, X. The Temporal and Spatial Evolution of Water Yield in Dali County. Sustainability 2015, 7, 6069-6085.

38. Liu, Y.; Xiao, J.; Ju, W.; Xu, K.; Zhou, Y.; Zhao, Y. Recent trends in vegetation greenness in China significantly altered annual evapotranspiration and water yield. Environ. Res. Lett. 2016, 11, 094010.

39. Tao, S.; Fang, J.; Ma, S.; Cai, Q.; Xiong, X.; Tian, D.; Zhao, X.; Fang, L.; Zhang, H.; Zhu, J.; et al. Changes in China's lakes: Climate and human impacts. Natl. Sci. Rev. 2020, 7, 132-140.

40. Tao, S.; Fang, J.; Zhao, X.; Zhao, S.; Shen, H.; Hu, H.; Tang, Z.; Wang, Z.; Guo, Q. Rapid loss of lakes on the Mongolian Plateau. Proc. Natl. Acad. Sci. USA 2015, 112, 2281-2286.

41. Lopez-Moreno, J.I.; Vicente-Serrano, S.M.; Moran-Tejeda, E.; Zabalza, J.; Lorenzo-Lacruz, J.; Garcia-Ruiz, J.M. Impact of climate evolution and land use changes on water yield in the ebro basin. Hydrol. Earth Syst. Sci. 2011, 15, 311-322.

42. Tao, S.; Zhang, H.; Feng, Y.; Zhu, J.; Cai, Q.; Xiong, X.; Ma, S.; Fang, L.; Fang, W.; Tian, D.; et al. Changes in China's water resources in the early 21st century. Front. Ecol. Environ. 2020, 18, 188-193.

43. Droogers, P.; Allen, R.G. Estimating Reference Evapotranspiration under inaccurate data conditions. Irrig. Drain. Syst. 2002, 16, 33-45.

44. Bosch, J.M.; Hewlett, J.D. A review of catchment experiments to determine the effect of vegetation changes on water yield and evapotranspiration. J. Hydrol. 1982, 55, 3-23.

45. Farley, K.A.; Jobbagy, E.G.; Jackson, R.B. Effects of afforestation on water yield: A global synthesis with implications for policy. Glob. Chang. Biol. 2005, 11, 1565-1576.

46. Calder, I.R. Water use of eucalypts-A review with special reference to South India. Agric. Water Manag. 1986, 11, 333-342.

47. Wang, Y.H.; Yu, P.T.; Xiong, W.; Shen, Z.X.; Guo, M.C.; Shi, Z.J.; Du, A.; Wang, L.M. Water-yield reduction after afforestation and related processes in the semiarid Liupan Mountains, Northwest China. JAWRA J. Am. Water Resour. Assoc. 2008, 44, 1086-1097.

48. Piao, S.; Yin, G.; Tan, J.; Cheng, L.; Huang, M.; Li, Y.; Liu, R.; Mao, J.; Myneni, R.B.; Peng, S.; et al. Detection and attribution of vegetation greening trend in China over the last 30 years. Glob. Chang. Biol. 2015, 21, 1601-1609.

49. Yu, P.; Wang, Y.; Wu, X.; Dong, X.; Xiong, W.; Bu, G.; Wang, S.; Wang, J.; Liu, X.; Xu, L. Water yield reduction due to forestation in arid mountainous regions, northwest China. Int. J. Sediment Res. 2010, 25, 423-430.

50. Yang, L.; Wei, W.; Chen, L.; Chen, W.; Wang, J. Response of temporal variation of soil moisture to vegetation restoration in semi-arid Loess Plateau, China. Catena 2014, 115, 123-133.

51. Ning, T.; Li, Z.; Feng, Q.; Liu, W.; Li, Z. Comparison of the effectiveness of four Budyko-based methods in attributing long-term changes in actual evapotranspiration. Sci. Rep. 2018, 8, 12665.

52. Serrat-Capdevila, A.; Scott, R.L.; Shuttleworth, W.J.; Valdes, J.B. Estimating evapotranspiration under warmer climates: Insights from a semi-arid riparian system. J. Hydrol. 2011, 399, 1-11.

53. Liu, Y.; Zhou, Y.; Ju, W.; Chen, J.; Wang, S.; He, H.; Wang, H.; Guan, D.; Zhao, F.; Li, Y.; et al. Evapotranspiration and water yield over China's landmass from 2000 to 2010. Hydrol. Earth Syst. Sci. 2013, 17, 4957-4980.

54. Tarboton, D.G.; Goeking, S.A. Forests and Water Yield: A Synthesis of Disturbance Effects on Streamflow and Snowpack in Western Coniferous Forests. J. For. 2020, 118, 172-192.