The Impact of Land Use Change on Water-Related Ecosystem Services in the Bashang Area of Hebei Province, China

Land use change is an important scientific issue recognized for its potential to alter ecosystem services (ESs), especially water-related ecosystem services (WRESs). Using the integrated valuation of ecosystem services and trade-offs (InVEST) model, this study quantified and mapped spatiotemporal variations in land use and corresponding WRESs in the Bashang area of Hebei Province, China (BAHP) to investigate how land use change impacted WRESs by means of scenario analysis, especially, in which a new evaluation indicator, average ecology effect (AEE) was proposed and well applied. The results indicated that woodland expansion (+602.61 km2) and grassland shrinkage (−500.57 km2) dominated the land use change in the BAHP in 2000–2018, which altered local WRESs, including the moderate declines in water purification and water yield, as well as a significant enhancement in soil conservation. In scenario analysis, compared to baseline levels, riparian woodland buffer and planting trees scenarios slightly decreased water yield but strengthened water purification and soil conservation; reclaiming wasteland and integrated development scenarios significantly enhanced soil conservation but lowered water yield and water purification; fertilizer reduction scenario effectively mitigated water deterioration. According to AEE, the riparian woodland buffer (RWB) scenario performed greater than the planting trees (PT) scenario on variations of WRESs per unit area, which differed completely from the results based on total variations. Overall, a multiple-scale indicator for a comprehensive evaluation of ESs should receive more attention.


Introduction
Ecosystem services (ESs) represent the varieties of benefits that ecosystem offered to humanity directly or indirectly to maintain living and welfare, in the form of culture, provisioning, regulating, and supporting services [1,2]. As the Intergovernmental Science-Policy Platform on Biodiversity and Ecosystem Services (IPBES) reported, under the influence of human activities, global ecosystems are suffering from ESs losses at an unprecedented speed [3], thus arousing growing attention in ESs worldwide [4,5]. The reduction of ESs, especially the water-related ecosystem services (WRESs), e.g., water yield [6], soil conservation [7], and water purification [8], would severely threaten human health and welfare. Hence, utilizing and maintaining ESs sustainably is a hot but challenging issue [9,10]. Land use change was an essential factor to recognize its potential impacts on ESs through altering land use type, pattern, and intensity [11,12]. For instance, the conversions of land use types of Hebei Province launched a new plan for land use in 2019 [32], aiming to facilitate the construction of water conservation areas and ecological environment support areas and ensure the local ecological security, as well as sustainable water resources. Therefore, studies on better understanding the impact of land use change on local WRESs is necessary to help policymakers formulate optimal schemes in advance and bridge the knowledge gap between WRESs and land use change.
In this study, three WRESs, including water purification, water yield, and soil conservation in the BAHP from 2000 to 2018, were quantified and mapped by the InVEST model. Besides, five alternative scenarios and a base landscape (baseline) based on land use in 2018 were developed into scenario analyses. The objectives of this study were (1) to assess the spatiotemporal dynamics of land use and WRESs; (2) to analyze the impacts of future scenarios on WRESs; and (3) to reveal the potential mechanism of how land use change impacted WRESs in arid and semi-arid area and provided a reference for related studies in other similar areas.

Study Area
The BAHP (40 • 70 -42 • 60 N, 113 • 70 -117 • 90 E) is located in the northwestern part of Hebei Province, China and across six counties (e.g., Guyuan, Zhangbei, Kangbao, Shangyi, Weichang, and Fengning), with a total area of 1.9 × 10 4 km 2 and an elevation range of 831-2215 m ( Figure 1). The BAHP is featured by numerous seasonal rivers and more lakes. Nowadays, most lakes and rivers are in danger of drying up because of a fast decline in groundwater levels led by human activities [31]. The ever largest plateau inland lake in northern China-the Anguli Lake, had totally degenerated into saline-alkali land in 2004. The BAHP belongs to continental monsoon plateau climate; the annual average temperature is 1-2 • C, and the annual average precipitation is approximately 400 mm, with nearly 70% of the rainfall appearing from June to September [30]. In the past several decades, the BAHP has experienced a long-term and high-intensity vegetable planting owing to agricultural transformation under the stimulation of the market economy, which causes a huge consumption of local water resources and seriously threatens water ecological environment. The BAHP is confronting a series of water-related ecological problems, including soil erosion, water shortage, water quality deterioration, and so on.

Data Sources and Pre-Processing
The main data used in this study and their sources are introduced as follows: (1) 30-m resolution land use datasets in 2000 and 2018 were from the Resource and Environment Science Data Center of Chinese Academy of Sciences; They are divided into

Data Sources and Pre-Processing
The main data used in this study and their sources are introduced as follows: (1) 30-m resolution land use datasets in 2000 and 2018 were from the Resource and Environment Science Data Center of Chinese Academy of Sciences; They are divided into six land use types, including cropland, woodland, grassland, water, construction land and unused land; (2) A digital elevation model (DEM) data, ASTER GDEM 30M, was obtained from the Geospatial Data Cloud (http://www.gscloud.cn); (3) 1-km resolution soil data, namely, Harmonized World Soil Database (HWSD), was from the National Cryosphere Desert center (http://www.crensed.ac.cn); (4) meteorological data (e.g., daily precipitation and evaporation data) were from China Meteorological Data Service Center (http://data.cma. cn), including the data of 11 meteorological stations in and around the BAHP; raster maps of meteorological data were derived from the interpolation utilizing the Kriging tool in ArcMap10.2 software; (5) the vector map of sub-watershed was extracted by Hydrological Analysis tool in ArcMap10.2 software.

Water Yield
Water yield is based on water balance principles, obtained by subtracting the actual annual evapotranspiration from the annual precipitation in each grid [33]. The calculation was based on the Budyko curve [34], and the water yield module algorithm in the InVEST model is as follows: where Y x is the water yield of grid x, AET x is the annual actual evapotranspiration on grid x, and P x is the annual precipitation in grid x.

Soil Conservation
Soil export, an opposite representation to soil conservation, is calculated based on the Universal Soil Loss Equation (USLE) [35]. In terms of sediment retention, the InVEST model takes the important hydrological process into account, in which sediment in the upper slope can be intercepted, making the final results more scientific and reasonable [33]. The average soil export rate on pixel x is quantified by the following: where, USLE x is the annual soil loss of grid x; R x is the rainfall-runoff erosivity of grid x, K x is the soil erodibility factor of grid x, LS x is the slope-length factor of grid x, C x is the crop management factor of grid x, and P x is the support practice factor of grid x.

Water Purification
Water purification, a critical regulating service, measured by the annual nitrogen or phosphorus exports based on the corresponding export coefficient of each land use type; it represents the capacity of the ecosystem in purifying water quality by retaining non-point source pollutants through the effect of soil and vegetation [33]. The average nitrogen export at grid-scale was calculated in the InVEST model by the following equations: where ALV x is the adjusted loading value of grid x; HSS x is the hydrologic sensitivity score of grid x, and pol x is the output coefficient of grid x; λ x is runoff index of grid x; λ x is mean runoff index in the study area; Σ U Y U is the sum of the water yield of the grid along the flow path above grid x.

New Indicator-Average Ecology Effect
In scenario analyses, a common indicator-total value variations of ecosystem services-was often utilized to evaluate how much one scenario affects ESs. Consequently, a consensus was that the smaller the total values of ESs changed, the more slight this scenario performed. In some cases, the area of some regions where land use was altered by one scenario rule is limited, thus causing a small fluctuation in total values of ESs, but actually, these local regions generated intense impacts on ESs under this scenario. To quantify this part of impacts, we proposed a new indicator, Average Ecology Effect (AEE). AEE aimed to quantify per unit area of ESs variations only for the part of land use change and to provide a new evaluation indicator to compare different scenarios' impacts on ESs. However, the limitation of AEE was that it was only applicable to the scenario where only one land use type came into being as a result of land use change. The formula of AEE is as follows: where AEE i is the unit area variation of ecosystem service i under one scenario, and its unit depends on the unit of ecosystem service i; ∆ES i is the total variation of ecosystem service i under one scenario; ∆S is the area variation of desired land use type under one scenario.

Scenario Setting
In order to simulate the future impacts of land use change on WRESs, this study designed five alternative scenarios relevant to local land use planning to estimate WRESs in each scenario; besides, land use in 2018 was assumed as a baseline for further comparison among different scenarios. The detailed process was displayed in Table 1.  (1) cropland = −0.76%; (2) woodland = +4.1%; (3) grassland = −1.01%; (4) unused land = −3.08%

Planting trees (PT) Scenario
Cropland and unused land that have a slope gradient of >15 • were converted into woodland.

Reclaiming wasteland (RW) Scenario
Unused land that has a slope gradient of <6 • were converted into cropland for agricultural development, but woodland and grassland with a slope gradient of <6 • were not converted due to eco-environmental protection.

Land Use Change in the BAHP from 2000 to 2018
Maps of land use in 2000 and 2018 were displayed in Figure 2. On the whole, six land use types kept decreasing except for woodland and construction land. Among them, cropland was the primary land use type, accounting for the highest proportion (49.27%) of the entire BAHP in 2018, and showed only a 1.07% decrease from 2000. Accordingly, although waters and construction land occupied the least proportion of the total area in 2018 (1.05% and 2.82%, respectively), their change rates were the highest compared with the other four land use types, reaching −31.75% and 43.16%, respectively. Over the past 18 years, the area of woodland expansion and that of grassland shrinkage were 602.61 km 2 (+16.60%) and 500.57 km 2 (−11.30%), respectively, dominating the most main land use change for the BAHP. Regarding unused land, it covered more than 929 km 2 in 2018 and decreased by 5.08% compared with 2000. Generally, land use change in the BAHP from 2000 to 2018 was characterized by a high change rate in both construction land and waters, and significant area change in woodland and grassland. From 2000 to 2018, the distribution pattern of land use and its spatial change varied among the six counties in the BAHP. In 2018, croplands were mainly distributed in Zhangbei (30.41%), Kangbao (24.86%), Guyuan (21.25%) and Shangyi (13.58%), and decreased by 31.55% and 48.26% in Shangyi and Zhangbei, with new 32.85% and 30.28% increases from 2000, respectively. Nearly 83.03% of Woodland and 53.63% of grassland were mainly observed in Fengning, Weichang, and Shangyi in 2018, respectively. Among them, grassland shrunk appeared primarily in western BAHP, and woodlands expansion was located in Shangyi (71.31%). Waters were mainly found in Guyuan (25.56%), Shangyi (24.90%), and Zhangbei (23.76%) whose most distinct change was the disappearance of Anguli Lake. Construction lands concentrated mainly in the urban core region of every county and extended outwards over the past two decades. Unused lands occurred mainly in Guyuan (40.71%) and Zhangbei (20.07%), and its spatial change in Zhangbei could be described as "a decrease in the east, but an increase in the west".
From 2000 to 2018, the conversion from cropland to woodland and grassland reached 230.97 km 2 and 191.19 km 2 , respectively ( Table 2). In addition, nearly 201.59 km 2 of construction land came mainly from cropland. In the same period, approximately 295.63 km 2 , 534.96 km 2 , and 130.25 km 2 of grassland were transformed to cropland, From 2000 to 2018, the distribution pattern of land use and its spatial change varied among the six counties in the BAHP. In 2018, croplands were mainly distributed in Zhangbei (30.41%), Kangbao (24.86%), Guyuan (21.25%) and Shangyi (13.58%), and decreased by 31.55% and 48.26% in Shangyi and Zhangbei, with new 32.85% and 30.28% increases from 2000, respectively. Nearly 83.03% of Woodland and 53.63% of grassland were mainly observed in Fengning, Weichang, and Shangyi in 2018, respectively. Among them, grassland shrunk appeared primarily in western BAHP, and woodlands expansion was located in Shangyi (71.31%). Waters were mainly found in Guyuan (25.56%), Shangyi (24.90%), and Zhangbei (23.76%) whose most distinct change was the disappearance of Anguli Lake. Construction lands concentrated mainly in the urban core region of every county and extended outwards over the past two decades. Unused lands occurred mainly in Guyuan (40.71%) and Zhangbei (20.07%), and its spatial change in Zhangbei could be described as "a decrease in the east, but an increase in the west".
From 2000 to 2018, the conversion from cropland to woodland and grassland reached 230.97 km 2 and 191.19 km 2 , respectively ( Table 2). In addition, nearly 201.59 km 2 of construction land came mainly from cropland. In the same period, approximately 295.63 km 2 , 534.96 km 2 , and 130.25 km 2 of grassland were transformed to cropland, woodland, and unused land, respectively. Moreover, the area of waters degenerating into unused land was 68.48 km 2 , which was caused primarily by the drying up of the Anguli Lake ( Figure 2c). Besides, newly grown grassland (140.56 km 2 ) and cropland (76.05 km 2 ) were both converted from unused land. It could be concluded that the local government had made some progress in revegetation practices. However, notably, such success was at the expense of a dramatic decline in grassland.

Evaluation of WRESs Based on Past Land Use Changes
Water yield was a critical provision service for the BAHP, a semi-arid area, and its supply shortage would seriously threaten the local water resources sustainability and ecology security. Driven by land use change, the water yield increased from 29.64 × 10 8 m 3 in 2000 to 29.76 × 10 8 m 3 in 2018, an increase of only 0.40% (Figure 3). This was mainly related to precipitation and land use change. On the one hand, water yield was highly sensitive to variations in precipitation [36], whereas this study used annual average precipitation to eliminate meteorological influences on water yield. On the other hand, there was a trade-off between afforestation and urbanization for their impacts on water yield, finally causing a slight net increase.

Evaluation of WRESs Based on Past Land Use Changes
Water yield was a critical provision service for the BAHP, a semi-arid area, and its supply shortage would seriously threaten the local water resources sustainability and ecology security. Driven by land use change, the water yield increased from 29.64 × 10 8 m 3 in 2000 to 29.76 × 10 8 m 3 in 2018, an increase of only 0.40% (Figure 3). This was mainly related to precipitation and land use change. On the one hand, water yield was highly sensitive to variations in precipitation [36], whereas this study used annual average precipitation to eliminate meteorological influences on water yield. On the other hand, there was a trade-off between afforestation and urbanization for their impacts on water yield, finally causing a slight net increase. Water purification was also modestly improved, as nitrogen export reached 6.37×10 6 kg in 2018, showing a 0.93% decrease from 2000. Although ever-increasing construction land produced high nitrogen exports, newly generated woodlands in practice intercepted more nutrients, resulting in a positive influence on water quality.
Soil conservation in the BAHP was significantly strengthened due to a 14.05% decrease in soil exports during the past 18 years. By 2018, total soil export was reduced to 5.14 × 10 6 tons. Generally, woodland possessed the strongest soil retention capacity among all land use types, however, it was the only land use type to maintain the growth Water purification was also modestly improved, as nitrogen export reached 6.37 × 10 6 kg in 2018, showing a 0.93% decrease from 2000. Although ever-increasing construction land produced high nitrogen exports, newly generated woodlands in practice intercepted more nutrients, resulting in a positive influence on water quality.
Soil conservation in the BAHP was significantly strengthened due to a 14.05% decrease in soil exports during the past 18 years. By 2018, total soil export was reduced to 5.14 × 10 6 tons. Generally, woodland possessed the strongest soil retention capacity among all land use types, however, it was the only land use type to maintain the growth trend except for construction land in BAHP during 2000-2018. Hence, woodland expansion was reasonably identified as the main factor for the enhancement of soil retention.

The Comparison of WRESs in Five Future Scenarios
In this study, we selected the land use in 2018 as a baseline. Figure 3 displayed the relative variations of WRESs among different alternative scenarios. The riparian woodland buffer (RWB) scenario performed well in soil conservation and water purification because of the slight declines in nitrogen export (−2.83%) and soil export (−0.58%), respectively, but it performed poorly on water yield because water supply decreased by 0.47% compared with the baseline level. The planting trees (PT) scenario generated similar effects on WRESs with the RWB scenario, and as compared with the baseline, soil export, nitrogen export, and water yield further decreased by 0.67%, 3.61%, and 1.75%, respectively. These were mainly because under the PT and RWB scenarios, the expansion of woodland led to the overall increases in evapotranspiration, nutrient interception capacity, and sediment retention capacity. Besides, for woodland, the litter layer on its soil surface also effectively slowed down the runoff accumulation and direct erosion from rainwater, thereby avoiding more sediment and nutrient losses.
As unused land with a slope below 6 • was transformed to cropland, the reclaiming wasteland (RW) scenario performed worse on water yield due to a 4.81% decrease compared with the baseline, which mainly resulted from higher evapotranspiration of vegetation coverage in cropland. Similarly, the RW scenario also showed a worse performance on water purification because nitrogen export was 7.70% higher than that in the baseline, at a total value of 6.86 × 10 6 kg/a, indicating that a high proportion of cropland generated more nitrogen sources and nutrient exports. However, different from the previous two, soil conservation in the RW scenario exhibited a significantly increasing trend by a dramatic decline of 15.95% (0.82 × 10 6 tons) in soil loss compared with the baseline level, which was mainly due to the large gap between cropland and unused land on soil retention capacity.
Regarding the integrated development (ID) scenario, it performed best in soil conservation due to a minimum value of 4.22 × 10 6 kg/a in sediment export but did worse in water yield, which was also because of a minimum value of 28.04 m 3 /a in the water supply. Compared with baseline, the ID scenario slightly degraded water purification as a result of a 1.41% increase in nitrogen export.
Finally, we added an extra fertilizer reduction (FR) scenario to investigate whether it would effectively mitigate water quality deterioration. The result was interesting because the FR scenario exhibited the largest decline in nitrogen export (−18.68%), far more than that in other scenarios. We attributed this to the high proportion of cropland in the BAHP (more than 49%) and its high nitrogen loading (20.02 kg/ha) before the implementation of fertilizer reduction. This result released a message that limiting the use of chemical fertilizers in cropland was an effective means for controlling pollutant emissions. Figure 4 showed the spatial distribution of per unit area of WRESs in different scenarios based on a subwatershed scale (Figure 1). For water yield, its reductions occurred in subwatershed 5 under the RWB scenario and subwatershed 4 under the PT scenario, respectively, because high evapotranspiration of woodland resulted in reductions of the surface water in these regions. In addition, water yield in the RW scenario followed the spatial changes similar to that in the ID scenario, and under these two scenarios, subwatershed 5 and 11 jointly experienced declines in water yield, as a higher proportion of unused land was replaced by cropland and contributed to more evapotranspiration compared to other subwatersheds. For water purification, there was no significant spatial difference observed in the distribution pattern of nitrogen export between the PT scenario and the baseline. Under the RWB scenario, the only decrease occurred in subwatershed 7, as the large area of woodland buffer was applied near waters in this region and strengthened the nutrient export control. In the RW scenario, the nitrogen exports increased mainly in For water purification, there was no significant spatial difference observed in the distribution pattern of nitrogen export between the PT scenario and the baseline. Under the RWB scenario, the only decrease occurred in subwatershed 7, as the large area of woodland buffer was applied near waters in this region and strengthened the nutrient export control. In the RW scenario, the nitrogen exports increased mainly in subwatershed 5, 8, 9, and 11 because unused land located here was almost converted to cropland, thus contributing to a larger amount of nutrient export. The ID scenario performed similarly to the RW scenario in the spatial change of nitrogen export, except in subwatershed 8 and 11, which indicated that vegetation restoration mitigated water quality deterioration in these two regions because of the RW scenario. Finally, the FR scenario showed a declining trend in the subwatershed 4, 5, 7, 8, and 11 compared to the baseline level, as these five subwatersheds covered almost the entire cropland of the BAHP, decreasing nitrogen loading in cropland and further limited nitrogen exports from nutrient sources.

The Spatial Distribution Pattern of WRESs in Baseline and Five Future Alternative Scenarios
Regarding soil conservation, as compared to the baseline level, spatial differences in soil export were not obvious enough to be captured in both RWB and PT scenarios. However, in RW and ID scenarios, the reductions in soil export were observed in subwatershed 2 and 5, and this was mainly because frequent conversions of unused land to cropland took place in these regions and generated more positive effects on soil erosion control.

Comparison of AEE Among Three Scenarios
The values of AEE in RWB, PT, and RW scenarios were calculated in Table 3. ID and FR scenarios are excluded due to the limitation of AEE. The AEE WY in RW scenario (−1629.33) was significantly lower than that in RWB (−931.84) and PT (−589.73) scenarios, which meant per unit area of the conversion from unused land to cropland in RW scenario exhibited a stronger negative effect on water yield than the corresponding conversions in both RWB and RW scenarios. Notes: WY = water yield; NE = nitrogen export; SE = soil export.
For water purification, the value of AEE NE in the RWB scenario was −11.98 kg/ha, 76.70% lower than that in the PT scenario, indicating that the afforestation efforts in riparian buffer performed better on the removal of nutrients than that in a steep slope. The value of AEE NE in the RW scenario reached 5.58 kg/ha, which represented that per hectare of the substitute of unused land by cropland would lead to an extra 5.58 kg nitrogen output.
Regarding soil conservation, the values of AEE SE in RWB, PT, and RW scenarios were −2.00, −2.65, and −9.34 t/ha, respectively. The RW scenario performed best in AEE SE compared with PT and RWB scenarios, and significantly improved soil erosion control after wasteland reclamation. The result suggested that in unit area, planting trees on steeper slopes could more effectively decrease soil loss than that in the riparian buffer.

Linking Land Use Change and WRESs in the Past 18 Years
Land use change affected WRESs mainly by altering the properties, processes, and components of the ecosystem (e.g., evapotranspiration, nutrient interception, and sediment retention capacity) [37]. Indeed, the potential mechanism between land use change and WRESs was more complex and still poorly understood [38]. In this paper, annual average values (e.g., precipitation and evapotranspiration) were used to eliminate the extremes of meteorological data and ensured land use change as the only influencing factor (similar to previous studies [22,39]), because we focused more on the impact of land use change on WRESs. During the past 18 years, the BAHP experienced drastic woodland expansion (+602.61 km 2 ) and grassland shrinkage (−520.57 km 2 ), including highly relative urbanization growth (43.16%), which indicated that afforestation projects in the BAHP conducted by the local policy-makers had made great progress and generated positive effects on water yield (a 0.40% increase), water purification (a 0.93% decline in nitrogen export), and soil conservation (a 14.05% decline in soil export). Most notably, these positive impacts on WRESs were purely derived from land use change in the BAHP. The study [20] pointed out that compared to land use change, climate change has a stronger influence on water yield and weaker impacts on soil conversation and nitrogen export. Another study [40] also demonstrated that climate factors had a highly positive relationship with water yield. The change levels among WRESs in this study were basically consistent with the above-mentioned study results. Here, it could be acknowledged that land use change has a stronger impact on soil conservation than that on water yield and water purification in the BAHP, an arid and semi-arid area.

Implications of WRESs in the BAHP Based on Scenario Analysis and AEE
In order to understand the impacts of one or several special land use conversions on WRESs, we created one baseline and five alternative scenarios to compare the output results of these scenarios by the InVEST model. In view of the total value changes of WRESs, RWB and PT scenarios slightly enhanced water purification and soil conservation but reduced water yield because of the expansion of woodland (similar to previous studies [22,23]). In addition, we found that the RW scenario, similar to the agricultural expansion scenario proposed in previous studies [22,23,41], generated an opposite conclusion-the RW scenario significantly enhanced soil conservation due to a significant reduction (15.95%) in soil export compared with the baseline level. This was primarily because the RW scenario did not transform vegetation types (i.e., woodland and grassland) to cropland, namely, a kind of special agricultural development strategy taking ecology protection into account. Although the RW scenario greatly improved agricultural production (a 9.29% increase in cropland), as a consequence, it aggravated water quality deterioration (a 7.70% increase in nitrogen export) because of high nitrogen loading in the cropland [42]. The ID scenario resulted in the largest reductions in both water yield (−5.78%) and soil export (−17.90%), and an increase of 1.41% in nitrogen export, which was basically accepted because more foods were provided at the expense of relatively slight water quality degradation. Notably, water yield reduction did not mean the total lessening of ESs provision, because existing studies demonstrated revegetation (including afforestation) could reduce water yield [43,44] and even generate positive effects on the hydrological cycle by accelerating the exchange of water vapor [45]. Finally, the additional FR scenario showed a massive reduction (−18.68%) in nitrogen export and significantly strengthened soil conservation, which meant that the FR scenario could be applied to the ID scenario and effectively made up for the deficiency of water purification, thus deserving to be the optimal measure to enhance WRESs provision in the BAHP.
In this study, in view of total value changes of WRESs, RWB and PT scenarios, similar to scenarios used in previous studies [23,41], generated relatively small impacts on WRESs. In reality, this was mainly because there were quite limited areas suitable for the implementations of RWB and PT scenarios in the BAHP. Besides, the PT scenario performed greater on WRESs due to the higher increase of woodland in the PT scenario (+7.28%) than that in the RWB scenario (+4.10%). To make the two scenarios comparable, not just depending on the only indicator such as total value change of WRESs, AEE was proposed to provide new insights on assessing the impacts of different scenarios on WRESs. The final results showed that the RWB scenario performed better than the PT scenario in per unit area of water yield reduction and water purification totally different from the results derived from total variations of WRESs in the previous section. This finding suggested an important message that inconsistent results may be produced based on different evaluation indicators when assessing one or more kinds of ESs. Such a suggestion could be reflected in this study in which afforestation around the riparian buffer was more appropriate in views of its improvement on water purification and hydrological cycle. Accordingly, afforestation on steep slopes was more beneficial to soil erosion control than that on the riparian buffer.
This study made contributions in highlighting the importance of spatial heterogeneity on issues about geography science and ecology science, and revealing a necessity for cautious assessment on ESs. Exactly, land use structure and topography often varied in different study areas, thus determining different land uses and ESs provisions from the same scenario. For example, the study [41] showed that cropland as the main land use type made agriculture expansion scenario (similar to RW scenario in this study) perform most greatly on WRESs. Moreover, when the study area [22,46] was featured with low average slopes and cropland was not the primary land use type, the agriculture expansion scenario also exhibited a stronger impact on WRESs than other scenarios. The abovementioned study, consistent with our results, jointly indicated spatial heterogeneity was worth considering in similar issues. Besides, multiple evaluation indicators could result in diverse conclusions because of spatial heterogeneity, for instance, there were distinct differences between AEE and total value variations in the explanations on how much various scenarios impacted WRESs. A similar result was also reasonably concluded based on another study [47], which exhibited that the PT scenario performed evidently better than the RWB scenario on soil conservation due to reductions of 50.64% and 30.05% in soil exports, respectively, but showed a narrow gap on enhancing soil conservation in views of AEE SE in PT (0.74 t/ha) and RWB (0.80 t/ha) scenarios. As available data was quite limited in other similar studies, more verification of AEE needs to be conducted in the future.

Limitations and Prospects
The InVEST model simplified hydrological processes, therefore there were many uncertainties for realistic simulation. For instance, the water yield model did not consider the annual variation of runoff generation and concentration. However, the generation of hydropower depends on the flow regulation, the nutrient retention model did not involve any biochemical reaction, and the sediment delivery ratio model also ignored the sediment sources from landslides and river courses [33]. Despite the limitations, the InVEST model was widely identified as a suitable tool for multiple ecosystem services assessments due to less input data and relatively high accuracy [48,49]. Due to the absence of data in the BAHP, a further field study needed to be carried out to validate the results calculated by the InVEST model in the future. Even so, approaches and results in this study still supported a more comprehensive understanding of the driving mechanisms of how land use change impacted WRESs and informed decision-makers and stakeholders of alternative land use planning to serve optimal WRESs provision.
ESs assessment was the frontier and hot spot of ecology science and geography science, and it was the bridge and link between natural processes and human processes [12]. However, the challenges on ecosystem service science still exist, including that (1) formulating a scientific framework used for explaining the relationship between ecosystem processes and services under different land use change driven scenarios; (2) developing multiple indicators to comprehensively assess ESs provision, trades off and synergies, driver factor and so on; and (3) integrating ESs assessments into spatial planning to accomplish the combination of science-policy because all are difficult and require new interdisciplinary knowledge.

Conclusions
This study evaluated land use change and three WRESs in the BAHP from 2000 to 2018, including water yield, soil conservation, and water purification (in terms of soil export and nitrogen export). In addition to scenario analysis, a new indicator was put forward in order to further compare the impacts of different scenarios on WRESs. The results indicated that from 2000 to 2018, the BAHP mainly experienced a significant increase in woodland (602.61 km 2 ) converted primarily from grassland (534.96 km 2 ) and a rapid urbanization expansion (+43.16%) at the expense of cropland and grassland, ultimately leading to the corresponding alterations in WRESs. From 2000 to 2018, WRESs in the BAHP were characterized by a slight improvement in water yield (+0.40%), significant enhancement in soil conservation (a 14.04% decrease in soil export), and a moderate degeneration in water purification (a 6.69% decrease in soil export). For scenario analysis, in views of the total value variations of WRESs, as compared to the RWB scenario, the PT scenario performed better in water purification and soil conservation but had worse performance on water yield. In addition, compared with the baseline, both PT and ID scenarios performed significantly well in soil conservation because of obvious declines (more than 15.95%) in soil export but performed relatively poorly in water yield and water purification. We found that the FR scenario effectively mitigated water quality deterioration as a result of an obvious decline (−18.68%) in nitrogen export, which suggested that nutrient loading parameters in the InVEST model had strong impacts on output results of nutrient export. Notably, the new indicator, AEE, indicated the opposite results. AEE indicated that the RWB scenario actually performed greater on water yield and water purification compared with the PT scenario, as both AEE WY (−931.84) and AEE NE (−11.98) in the RWB scenario were smaller than AEE WY (−589.73) and AEE NE (−6.78) in the PT scenario, which was different from conclusions drawn by total value variations of WRESs in the previous section. Based on the above-mentioned analysis, three important points were made: (1) Afforestation in the riparian buffer could improve relatively higher ESs provision than that on steep slopes in views of the unit area; (2) Integrating the FR scenario into the ID scenario could effectively mitigate more nitrogen exports in the ID scenario; (3) Multiple-scale indicators should be developed into a comprehensive assessment on ESs, so as to avoid a one-sided conclusion.  Data Availability Statement: Data sharing not applicable. No new data were created or analyzed in this study. Data sharing is not applicable to this article.