Future Impacts of Climate Change and Land Use on Multiple Ecosystem Services in a Rapidly Urbanizing Agricultural Basin, China

Abstract

Ecosystem services (ESs) in rapidly urbanizing agricultural basins are vulnerable to environmental changes. Adequately understanding the driving forces and the dynamics of ESs related to water quantity and quality can provide a basis for making sound management decisions on the development of basins. Here, we explored the impacts of future land use and climate changes on four ESs: nitrogen and phosphorous purification, water supply, and soil retention services in the Taihu Basin region of eastern China. Spatially explicit methods, a cellular automata-Markov (CA-Markov) model and the delta downscaling method were used to quantify the ESs, simulate land use changes, and project future climate changes, respectively. We built a business-as-usual land use scenario, representative concentration pathways (RCPs) scenarios for climate change, as well as a combined land use and climate change scenario to analyze the changes in the drivers and the responses of ESs. The results showed the following: (1) future land use changes would significantly enhance the nitrogen purification service while reducing the phosphorus purification service compared to other services; (2) climate change would have substantial effects on water supply and soil retention, but these impacts would vary with different RCPs scenarios during three future periods; and (3) the combined scenarios of both drivers would obviously influence all ESs and lead to a nitrogen purification service that was different from the other three services. Moreover, the policy implications of the results were discussed. The findings can help guide the creation of policies for land structure and patterns, climate change adaptation, and ecosystem-based management to promote the sustainable development of watersheds at the regional scale.

1. Introduction

Ecosystem services (ESs) have drawn considerable attention since the publication of the research by Costanza (1997) [1], and numerous studies have been carried out to quantify and map the supply of ESs [2,3]. The circulation of material and flow of energy on Earth are altered with climate change, population growth, technological progress, and urbanization [4], which complicate the spatial patterns and temporal changes in ESs. Understanding the mechanisms that drive these changes is essential to maintain ecosystems as dynamic and adaptive systems to promote the supply of ESs [5]. Moreover, this understanding is a key step in integrating ESs into actual management that takes driving forces as the control factors [6,7].

The driving mechanisms underlying the supply of ESs are complex due to the various driving forces at different spatial and temporal scales [8,9,10]. Climate change and land use change are commonly considered two significant driving forces that influence ESs [11,12]. Most studies on climate change have used climate models [13,14,15,16] to analyze the factors of temperature, precipitation, and greenhouse gas emissions [17,18,19,20] at multiple spatial scales [21,22,23]. The results indicated that there were mainly negative impacts of climate change on most types of ESs, such as food provision [24], water yield [11,25], and soil conservation [4], but carbon sequestration was highly variable across studies and ecosystem types [26,27]. Land use change impacts multiple ESs by altering ecosystem patterns and processes. Recent studies have shown that land use change caused by agricultural expansion, urbanization, and grazing may lead to the reduction in carbon sequestration, habitat quality, and water-related services [28,29,30,31], while increases in protected areas, woodlands, wetlands, and other ecological lands can promote ESs [32,33,34].

Integrating regional land use with climate change is critical for understanding the mechanisms driving ESs. Assessments of historical data and analyses of future scenarios have been applied to analyze the combined driving forces of multiple ESs at different spatial scales [11,35,36,37] and in different ecosystems [9,12,18]. Most results have demonstrated that a great number of ESs are highly sensitive and vulnerable to environmental changes. Although previous studies have examined the potential effects of land use and climate change on ESs, most predominant drivers such as climate change, land use change, and management policies remain in question, especially at a regional level [18,38]. Moreover, the policy implications of the relevant studies have not been sufficiently analyzed. Therefore, the integrated effect of land use and climate changes on ESs remains an important area for future studies [20,39].

As areas with intricate hydrologic processes and intensive human activities, watershed ecosystems are dramatically affected by land use and climate changes, thereby influencing the provision of water-related ESs [40]. In China, many watersheds suffer from critical water issues caused by water shortages and water quality degradation due to the negative consequences of land use and climate changes. It is crucial for watershed managers, planners, and scientists to better understand the responses of multiple ESs to a single driving force along with the interaction of both drivers. Nonetheless, limited information is available on the driving forces of ESs at a regional scale. Previous studies have quantitatively assessed the impacts to water quantity and quality service from climate change or land use separately [25,29,41], whereas only a few studies have compared the combined impacts of both drivers on water-related services [37,40,42]. Consequently, research on the combined drivers of water-related services in a regional watershed should be strengthened to provide supporting information to guide human activities, improve the supply of multiple ESs, and promote sustainable development of watershed ecosystems [43,44].

Given the demand for ES studies in regional watersheds, we examined the responses of four ESs (nitrogen purification, phosphorus purification, water supply, and soil retention) to future land use and climate changes. This study addresses the following research questions: what are the potential influences on multiple ESs of future land use changes, different climate scenarios, and the combination of these drivers? In addition, we presented a sufficient analysis of a feasible policy to improve ESs. We used multiple spatially explicit models to quantify the four ESs, a CA-Markov model to predict land use changes in 2040, and the delta downscaling method to project climate factors under three representative concentration pathways (RCPs) scenarios until 2100. We selected the Jiangsu section of the Taihu Basin region as our study area because it is a typical agricultural area accompanied by rapid urbanization that is similar to many parts of the world and has become one of the most developed regions in China. We hope our study will provide scientific support for freshwater management and sustainable use of ESs in this region.

2. Methods

2.1. Study Area

Jiangsu section of the Taihu Basin region (30°45′36″–32°15′04″ N, 119°03′14″–121.17′56″ E) covers a total area of 19,289.52 km² and is located in Yangtze River Delta of Eastern China (Figure 1). This area includes Suzhou, Wuxi, Changzhou, and a small part of Zhenjiang and Nanjing City, which belong to an independent water conservancy region delineated by the Chinese Ministry of Water Resources. The region has a typical subtropical monsoon climate with an annual mean temperature of 15–17 °C and annual mean rainfall of 1180 mm [10,45]. The study area is one of the most populous and developed regions of China, with a population of more than 22 million. Moreover, this area produced 5.2% of China’s gross domestic product (GDP) in 2010 according to the 2010 China census. The dominant land cover is cultivated land, covering 45.97% of the study area. Due to rapid urban sprawl, the area of agricultural land has been decreasing since 1985 [46,47]. Forest and grassland, accounting for 5.55% and 0.13%, respectively, are mainly located in the mountainous areas in the south and around lakes (Figure 1).

As the largest lake in the study area, Taihu Lake (2238.1 km²) is also the third largest freshwater lake in China [48]. Taihu Lake plays an important role in water provision for Shanghai Municipality as well as the surrounding cities. However, because of global warming, population increases, and rapid urbanization, this region has faced serious degradation of freshwater ecosystems together with frequent outbreaks of eutrophication since the 1990s, which are primarily caused by nonpoint source pollution [49,50]. In general, the study area is one of the most typical regional watersheds that are faced with competition between clean water demands and insufficient water-related ESs in China.

2.2. Data Sources and Processing

The data used for the ES model and the analysis of the driving forces can be divided into six types as follows: (1) Climate data, including daily and annual temperature, precipitation, evaporation, and solar radiation, which were derived from the China Meteorological Data Sharing Service System (http://cdc.nmic.cn/). We applied ordinary kriging interpolation to convert the site data into spatial raster maps. The Coupled Model Intercomparison Project phase 5 (CMIP5) multi-model dataset was collected, analyzed, and provided by the National Climate Center. (2) Land use vector data from 2000, 2005, and 2010 were provided by the National Earth System Science Data Sharing Infrastructure (http://www.geodata.cn), and GIS data transformation tools were used to process the data. (3) Digital elevation model (DEM) data with a resolution of 30 m × 30 m were downloaded from the International Scientific Data Service Platform (http://www.gscloud.cn) and filled using the ArcGIS hydrology model to remove sinks. (4) The 100 m × 100 m soil data, root restricting layer depth, and plant available water content (PAWC) were provided and processed by Nanjing Agricultural University and the Institute of Soil Science, Chinese Academy of Sciences [51]. (5) Water quality, hydrologic, and other statistical data were acquired from the local environmental protection, water conservancy, and statistical departments of each county. (6) The biophysical parameters for the ES model were obtained from regional related publications, local government documents, as well as the Integrated Valuation of Ecosystem Services and Trade-offs (InVEST) user guide [52], and the details were specified in the ES quantification. All spatial data were unified to a 30 m spatial resolution grid using ArcGIS 10.2 (Environmental Systems Research Institutes, Inc., Redlands, CA, USA).

2.3. Ecosystem Services Quantification

Due to the major threats to water quality and water yields in the Taihu Basin region, we studied four important ESs, including nitrogen and phosphorus purification services and water supply and soil retention services. A set of biophysical indicators that represent the ecosystem functions that benefit humans were used as ES metrics [53]. Biophysical models and empirical estimations were applied to quantify and map the ES indicators. The details are introduced as follows.

2.3.1. Water Supply

Hydrologic cycles in ecosystems generate water supply (WS) services for human consumption or hydropower production. As the indicator of the WS service, annual water yield was quantified by the water yield model in the InVEST tool, which considers only surface water yield and ignores the influence of groundwater. The water yield in each pixel of the region WS(x) is expressed by Equation (1):

$$WS\left(x\right)=\left(1-\frac{AE{T}_x}{Pr{e}_x}\right)\times Pr{e}_x$$

(1)

where $Pr{e}_x$ and $AE{T}_x$ represent the annual precipitation and actual evapotranspiration in pixel x (mm); AET_x/Pre_x is expressed based on the Budyko curve [54,55] as shown in Equation (2):

$$\frac{AE{T}_x}{Pr{e}_x}=\frac{(1+{\omega }_x+R_x)}{(1+{\omega }_x+\raisebox{1ex}{$1$}\!\left/ \!\raisebox{-1ex}{$R_x$}\right.)}$$

(2)

where ${\omega }_x$ is a nonphysical parameter, and $R_x$ represents the aridity index. They are calculated by Equations (3) and (4):

$${\omega }_x=Z\times \left(\frac{AW{C}_x}{Pr{e}_x}\right)$$

(3)

$$R_x=\frac{k_j\times E{T}_0}{Pr{e}_x}$$

(4)

where Z is an empirical constant that represents local precipitation patterns and hydrogeological characteristics; AWC_x is the volumetric PAWC (mm) estimated by Equation (5). k_j is the evapotranspiration coefficient of each land type, which was acquired from previous studies [56,57,58] as well as the InVEST guidebook. $E{T}_0$ is the potential evapotranspiration (mm/d) quantified by the modified Hargreaves equation [59] with meteorological station data.

AWC_x = Min(Rest.layer.depth, root.depth)·PAWC

(5)

where Rest.layer.depth is the raster data of soil depth (mm) at which root penetration occurs; root.depth is the depth of vegetation rooting for each land type (mm), and the parameter values were assigned according to the results of Liu [60], Huang [61], and Lv [62].

The performance of the water yield model was tested in two subwatersheds (Meilin and Lihe) for which dense and consistent hydrological monitoring data have been available over the last decade, and Nash–Sutcliffe model efficiency coefficients were used to assess the performance. Previous studies [63,64] have indicated that the accuracy of the input data such as precipitation and potential evaporation data largely determines water yield results. Therefore, different interpolation methods were used to model the spatial distribution of input data in order to improve the performance of the model.

2.3.2. Nutrient Purification Services

The ecological factors in terrestrial ecosystems can contribute to water quality improvements by removing nutrient pollutants from runoff, and nutrient loadings of the ecosystem are widely used as inverse proxies of water purification services. In this study, annual nitrogen loading (NL) and phosphorus loading (PL) were assessed by the nutrient purification model in the InVEST tool [53,65], and low loading values corresponded to high nutrient purification services. The model can be expressed as shown Equation (6):

$$Ex{p}_x=AL{V}_x{\displaystyle \prod }_{y=X+1}^X\left(1-E_y\right)$$

(6)

where Exp_x is the nutrients exported from the upper grid x (kg/ha); X represents the nutrient transport route; E_y is the filtration efficiency downstream at pixel y (%); and $AL{V}_x$ represents the nutrient loading at pixel x (kg/ha), which can be calculated by Equation (7):

$$AL{V}_x=HS{S}_x\times po{l}_x$$

(7)

$$HS{S}_x=\frac{{\mathsf{\lambda }}_x}{{\mathsf{\lambda }}_w}$$

(8)

$${\mathsf{\lambda }}_x=\log ({\displaystyle \sum }_U{Y}_u)$$

(9)

where $po{l}_x$ is the nutrient export coefficient at pixel x; $HS{S}_x$ represents the sensitivity value of hydrology calculated by Equation (8), where ${\mathsf{\lambda }}_w$ is the average runoff coefficient of the entire region and ${\mathsf{\lambda }}_x$ is the runoff coefficient at pixel x, which can be quantified by Equation (9). ${\displaystyle \sum _U}{Y}_u$ is the total water yield flowing into pixel x, which can be calculated by the water yield model in the InVEST tool.

The major parameters of the model are the nitrogen and phosphorus export coefficients and filtration efficiency of each land use type as well as the flow accumulation threshold, which were derived from published field studies and local estimates [49,66,67,68].

The model was calibrated by comparing our results with study results of Huang [69] and Wang [70] in the region based on previous studies [71,72]. Because nutrient export coefficients and filtration efficiency are the sensitive factors of the model [73,74], we adjusted these parameters within a reasonable range until the deviation was within 20%.

2.3.3. Soil Retention Service

Ecological factors (such as vegetation cover) in ecosystems strengthen the function of preventing the erosion and input of soil to streams, helping to maintain the capacity of the soil to filter pollutants and regulate water quality. The annual potential reduction in soil loss was used as the indicator of the soil retention (SR) service, which was quantified by the revised universal soil loss equation (RUSLE) [75]. As the model may generate deviations from the actual conditions of the study area, we adjusted the model based on experimental data from the subwatersheds. The SR model can be expressed as Equation (10):

$$SR\left(x\right)=R\times K\times LS\times \left(1-C\times P\times A\right)$$

(10)

where R is the rainfall erosivity index (MJ·mm·(ha·h·a)⁻¹) calculated by Equation (11); K is the soil erodibility (t·ha·h·(ha·MJ·mm)⁻¹) calculated by Equation (12); LS is a dimensionless parameter representing the slope length and steepness index, which can be quantified by Equation (13); C is the cover-management factor quantified by Equation (14); P represents the support practice factor with a range from 0 to 1, which was determined based on the study by Zhang [76]; and A is the adjustment factor of the model, and was set as 1.3 in our study area.

$$R={\displaystyle \sum }_{i=1}^{12}1.735\times {10}^{[1.5\times \lg \left(\frac{p_i^2}{p}\right)-0.8188]}$$

(11)

$$K=\left\{0.2+0.3\times \exp \left[-0.0256S_a\left(1-\frac{S_i}{100}\right)\right]\right\}{\left(\frac{S_i}{C_i+S_i}\right)}^{0.3}\times \left[1-\frac{0.25\times C}{C+\exp \left(3.72-2.95C\right)}\right]\times \left[1-\frac{0.75\times S_n}{S_n+\exp \left(2.95\times S_n-5.51\right)}\right]$$

(12)

$$LS=\{\begin{array}{c}{\left(\frac{F_a\times C_s}{22.13}\right)}^n\times {\left[\frac{\sin \left(s\times 0.01745\right)}{0.09}\right]}^{1.4}\times 1.6, s\le 25°\\ 0.08\times {\beta }^{0.35}{s}^{0.6}, s>25°\end{array}$$

(13)

$$C=\{\begin{array}{c}1,c=0\\ 0.6508-0.3436\times \lg c, 0<c\le 78.3\%\\ 0,c\ge 78.3\%\end{array}$$

(14)

In Equation (11), p is annual average precipitation (mm), and p_i is monthly average precipitation (mm). In Equation (12), S_a, S_i, C_i, and C are the values of sand, silt clay, and organic carbon contents (%), respectively. In Equation (13), F_a is the flow accumulation threshold, C_s is the grid size, s is slope, and n and $\beta$ are the associated parameters of the slope and flow direction, respectively. In Equation (14), C is a nondimensional parameter with a range from 0 to 1, and c is vegetation coverage that can be can be determined by the normalized difference vegetation index (NDVI) [77].

For the accuracy assessment, we compared our estimates with field measurements of soil erosion in the Meilin subwatershed as well as the results of previous similar studies [78,79] in this region.

2.4. Driving Scenario Simulation

2.4.1. Land Use Change

We used a cellular automata-Markov model to predict land use changes in 2040 based on land use status. The model is the integration of Markov chain analysis (MCA) and a cellular automata (CA) model with the advantages of both long-term predictions and complicated spatial variation simulations, which can help to understand the temporal and spatial changes in land use [80]. MCA applies the occurrence probability of random events in Markov theory to predict future secular changes, and the model can be expressed by Equation (15):

$$R\left(t+1)=R\right(t)\cdot P$$

(15)

This model indicates that the next state $R(t+1)$ can be predicted when the original state $R\left(t\right)$ and transition probability P of land use are known. However, MCA cannot predict the spatial patterns of land use changes because it lacks the description of spatial information. However, the CA model can narrow this gap due to its calculation principle as a dynamic network model that assumes that limited and discrete states of the variable will change in space under different rules. The key factors of the CA model are cell, state, lattice, neighbor, and rule, and their action mechanisms can be expressed by Equation (16):

$$S_{t+1}=f\left(S_t,N\right)$$

(16)

When predicting land use changes, a cell is the image grid; S is the spatial attribute (state set) of the cell, including spatial position, type, size, and perimeter; N is a neighbor of the cell whose scope can be defined by filters at different scales, such as $3\times 3$ and $5\times 5$; f is a state change rule, i.e., the constraint function of land use change determined by suitability layers of distance, elevation, ecological protection area, etc.; t and t + 1 represent different times.

We used IDRISI software to run CA-Markov model [81,82] due to its accurate prediction of quantitative structure and spatial patterns of land use. Raster data of land use in 2000, 2005, and 2010 were used to generate transition area and probability matrixes as well as a series of transition probability maps for eight land use types that were reclassified from original land use map. Then CA model was used to predict likely changes of spatial distribution for each land types. The transition of each pixel gives priority to contiguous suitable areas defined by $5\times 5$ spatial filters and suitability layers. There were three constraint conditions when constructed suitability layers: (1) Based on related water governance policies in Taihu Lake Basin, such as the “Water Pollution Control Regulations of Taihu Lake in Jiangsu Province”, we classified regions as unsuitable for development if they within 5 km from Taihu Lake, 1 km from rivers flowing into Taihu Lake and 200 m from a common water body. (2) According to the stipulations of principal function regionalization, ecological protection, and cultivated land policies, the forbidden development zones, ecological red line regions, and prime cropland areas were considered unsuitable for development. (3) The areas with slopes greater than 25° were also considered unsuitable for development in light of national and local regulations on soil and water conservation. Finally, we chose 2010 land use map at initial stage of the prediction and 40 iterations to simulate land use in 2040. Overall, the CA-Markov model makes the assumption that the change trends of land use in all time periods are same, and land use in 2040 was considered as the business-as-usual scenario of land use change.

To verify the accuracy of the CA-Markov model, we used the land use maps from 2000 and 2005 to simulate the 2010 land use and then compared it with the 2010 actual land use map. The precision of the model was quantified using Equations (17):

$$e_i=\frac{|x_{im}-x_{in}|}{x_{in}}$$

(17)

where $e_i$ is the error precision of the i land use type, and a lower value means a higher prediction accuracy; $x_{im}$ and $x_{in}$ are the simulated and actual areas of the i land use type, respectively.

2.4.2. Climate Change Scenarios

As the latest greenhouse gas (GHG) emission scenarios were shown in the Intergovernmental Panel on Climate Change (IPCC) Fifth Assessment Report (AR5) report, RCPs were used to project future climate changes in the study area. RCPs scenarios (RCP2.6, 4.5, 6.0, and 8.5) are built based on different atmospheric radiation intensities based on diverse hypotheses of population growth, economic development, technical progress, environmental conditions, and other factors. We selected the RCP2.6, 4.5, and 8.5 scenarios in our study because they represent the minimum, middle, and maximum emissions levels, respectively, and reflect three important changes in integrated driving forces.

The Coupled Model Intercomparison Project phase 5 (CMIP5) multi-model dataset was used to analyze future climate changes under the RCP2.6, 4.5, and 8.5 scenarios. The dataset has a resolution of 1° × 1° and is the result of 21 global climate models (GCMs) that were simulated through interpolation and downscaling processes, including average, maximum and minimum temperatures and monthly precipitation during a historical period (1971–2005) as well as a simulated period (2006–2100). Since the spatial resolution of the dataset is too coarse to represent the local climate characteristics in the study area, we applied a simple delta downscaling method [40,83] to increase the resolution of the data. The methods for downscaling the precipitation and temperature data are expressed in Equations (18) and (19), respectively:

$$P{\prime }_m=P_m\times \left(\frac{P_{m,future}}{P_{m,current}}\right)$$

(18)

$$T{\prime }_m=T_m+(T_{m,future}-T_{m,current})$$

(19)

where $P{\prime }_m$, $P_m$, $T{\prime }_m$, and $T_m$ are future mean monthly precipitation (mm), current mean monthly precipitation (mm), future mean monthly temperature (°C), and current mean monthly temperature (°C) of each meteorological station in the study area; $P_{m,current}$, $T_{m,current}$, $P_{m,future}$ and $T_{m,future}$ are the simulated mean monthly precipitation and temperature under current and future climate scenarios, respectively. Kriging interpolation was used to calculate the future climate throughout the region based on meteorological station data.

According to the above procedures, we simulated climate factors from 2006–2010 and compared our simulated results with the actual data to assess the accuracy of the climate change scenarios.

3. Results

3.1. Calibration of Ecosystem Service Model

We calibrated the ESs model using the methods described above, and the results showed that our ES assessments were reasonably consistent with the field data or relevant studies. For the WS service, we used kriging interpolation to produce a 30-m pixel raster for the input data for the water yield model by comparing modeled data with observation data in subwatersheds. The Nash–Sutcliffe model efficiency coefficients ranged from 0.52 to 0.85, which indicated that there was consistency between our estimate and the averaged water yields of the subwatersheds. The deviation between the results simulated by the nutrient purification model and previous studies [69,70] was controlled within 20%, demonstrating that the model was fit for the study area. The simulated NL result was particularly accurate as the precision reached 87.6%. Regarding the SR service, the field-measured soil erosion values were higher than the simulated results by 9.4% in the Meilin subwatershed. Moreover, the spatial distribution of our estimates was consistent with the findings of Zhang [78] and Zheng [79], further showing the feasibility of the model.

3.2. Future Changes in Land Use and Climate Factors

3.2.1. Land Use Change in 2040

The verification results of the CA-Markov model indicated that the overall accuracy of the predicted land use patterns in 2010 was 87.2%. The simulation accuracies of cultivated land, forestland, and water bodies were 90.5%, 89.6%, and 91.1%, respectively, while the simulation accuracies of urban construction land and rural residential land were relatively low at 85.7% and 86.4%, respectively. Overall, the simulated land use in 2010 met the minimum mapping accuracy requirements [84] and showed that the CA-Markov model was able to predict land use patterns in the future.

In view of the uncertainty of human activities and policy changes in the future, we simulated the land use patterns in 2040 to represent the medium-term changes in land use in the study area. Figure 2 showed the predicted results in 2040, and the total area of each land use type was given in

Table 1. The areas of land use types in 2010 and 2040.

Land Use Types	Land Use in 2010		Land Use in 2040
	Area (hm2)	Percentage (%)	Area (hm2)	Percentage (%)
Cultivated land	886,783.14	45.97	652,889.25	33.85
Forest land	107,012.43	5.55	138,185.46	7.16
Grass land	2520.27	0.13	2384.91	0.12
Water body	438,084.45	22.71	472,858.92	24.51
Wetland	4518.54	0.23	4008.06	0.21
Urban construction land	316,377.90	16.40	442,777.50	22.95
Rural residential land	168,790.50	8.75	208,786.14	10.82
Other land	4864.59	0.25	7061.58	0.37

. The results indicated that the area ratios of cultivated land (33.85%), water bodies (24.51%), and urban construction land (22.95%) were larger than those of other land use types. From 2000 to 2040, the area of cultivated land prominently decreased while urban construction land obviously increased, and there were different levels of increase in forestland, water bodies, as well as rural residential land.

3.2.2. Climate Change in Future Periods

According to the predicted climate changes, we simulated several climate factors under the RCP2.6, 4.5, and 8.5 scenarios from 2011–2100. Three periods of 2011–2040 (short-term change), 2041–2070 (medium-term change), and 2071–2100 (long-term change) were summarized to assess the changes in ESs. We calculated the average monthly and annual value for each climate factor in the three periods as well as in the historical period from 1981–2010, which was used as the baseline scenario. Annual precipitation and temperature (Figure 3) were taken as examples to analyze future climate changes.

The average annual precipitation results showed a significant reduction from 2011–2040, small changes from 2041–2070, and an increase from 2071–2100 compared with the baseline scenario. From the view of the different climate scenarios, the average annual precipitation in RCP2.6 and RCP8.5 exhibited limited differences in the periods of 2011–2040 and 2041–2070, while the average annual precipitation was slightly higher under the RCP4.5 scenario than that in the other scenarios in the two periods. The average annual precipitation increased gradually with the increase in carbon emissions from 2041–2070. Overall, there were distinct temporal changes, and distinct variations appeared due to the multiple carbon emission sources from human activities, particularly in the long term.

The results indicated that the average annual temperature would significantly and rapidly increase in the future under the high emission scenario. Compared with the baseline scenario (16.062 °C), the temperature would increase by 0.759–0.86 °C in 2011–2040, 1.199–2.385 °C in 2041–2070, and 1.132–4.075 °C in 2071–2100. There were fewer differences among the three emission scenarios in 2011–2040, but the differences were distinct in 2041–2070 and 2071–2100. In the RCP2.6 scenario, the average annual temperature slowly increased and fell in the long term. However, in the RCP8.5 scenario, the temperature rapidly increased and reached 20.137 °C in 2100, exceeding the values in the RCP4.5 and RCP2.6 scenarios by 2.023 °C and 2.943 °C, respectively.

3.3. Changes in Ecosystem Services under Different Driving Forces

3.3.1. Impacts of Future Land Use Changes

The changes in ESs under land use changes in 2040 were simulated by setting the climate factors to the average values from 1980–2010 and leaving the other factors unchanged. The ESs in 2040 were viewed as the business-as-usual scenario over time and compared with the baseline scenario (ESs in 2010) to analyze the impacts of land use change on ESs in the future (

Table 2. Total amounts of ecosystem services and their changes under the land use patterns in 2010 and 2040.

	Nitrogen Loading (t/year)	Phosphorus Loading (t/year)	Water Supply (mm)	Soil Retention (t/year)
ESs of 2010	24,207.19	2142.87	728.18	5,655,058.58
ESs of 2040	22,116.20	2228.97	715.46	5,741,420.19
Change rate (%)	−8.64	4.02	−1.75	1.53

and Figure 4).

The results showed that under the influence of land use change, NL decreased by 8.64% while PL increased by 4.02%, which indicated that there were obvious differences between nitrogen purification services (enhanced by 8.64%) and phosphorus purification services (reduced by 4.02%). WS and SR slightly decreased by 1.75% and increased by 1.53%, respectively. The spatial distribution indicated distinct spatial heterogeneity among the ESs. The middle NL and PL values presented a trend of extensive distribution, and the WS value varied significantly across the whole region, while the SR value exhibited no evident difference.

3.3.2. Impacts of Future Climate Changes

Assuming that the land use patterns remained the same, we simulated the annual average ESs in three periods under the RCP2.6, RCP4.5, and RCP8.5 climate change scenarios. The annual average ESs are shown in Figure 5, and the change ratios compared with the baseline scenario are listed in

Table 3. Change ratios of ecosystem services under the future climate change scenarios (%).

ESs	2011–2040			2041–2070			2071–2100
	RCP2.6	RCP4.5	RCP8.5	RCP2.6	RCP4.5	RCP8.5	RCP2.6	RCP4.5	RCP8.5
NL	0.250	0.187	0.240	0.180	0.219	0.339	0.151	0.196	0.475
PL	0.200	0.147	0.191	0.135	0.164	0.267	0.109	0.144	0.369
WS	−7.531	−5.521	−7.273	−2.168	−2.310	−4.123	−0.687	0.581	−3.092
SR	−15.441	−11.801	−14.946	−7.900	−8.386	−9.716	−7.842	−4.011	−2.942

.

The results indicated that compared to the baseline scenario, NL and PL would increase slightly, while WS and SR would decrease significantly, which means climate change had a negative influence on these ESs. For NL and PL, the average annual changes were smallest under RCP4.5, while the impacts of RCP2.6 and RCP8.5 exhibited little differences in 2011–2040. NL and PL increased from the low emission scenario to the high emission scenario in 2041–2070 and 2071–2100. From the view of the temporal changes in the emission scenarios, both NL and PL showed limited decreases under the RCP2.6 scenario whereas they showed obvious increases under the RCP8.5 scenario. Although WS decreased compared with the baseline, it gradually increased over time, and the increase was especially evident in the RCP2.6 and RCP4.5 scenarios. The WS in the RCP8.5 scenario was inferior to that in other scenarios in 2041–2070 and 2071–2100. SR also increased from 2011 to 2100 but exhibited distinct changes under the different scenarios in the three periods. The SR values of RCP4.5, RCP2.6, and RCP8.5 were maximum in 2011–2040, 2041–2070, and 2071–2100, respectively, which showed good consistency with the average annual precipitation under the climate change scenarios in 2011–2040 and 2071–2100.

3.3.3. Impacts of Future Land Use and Climate Changes

We simulated the combined effects of land use and climate changes on ESs in the medium term (2041–2070) because climate change had few impacts on land use during this period [85]. Due to the uncertainty of the physical environment and human activities in the future, the annual average climate factors in 2041–2070 together with land use in 2040 were used as input variables for the ESs model. The annual average amounts of ESs under the RCP2.6, RCP4.5, and RCP8.5 scenarios, as well as the change ratios compared with the baseline scenario, were shown in

Table 4. The annual average amounts and change ratios of ecosystem services under the combined effects of land use and climate changes in 2041–2070.

ESs	RCP2.6		RCP4.5		RCP8.5
	ES Amount	Change Ratio (%)	ES Amount	Change Ratio (%)	ES Amount	Change Ratio (%)
NL (kg/year)	22,147,373.39	−8.509	22,168,400.05	−8.422	22,197,212.21	−8.303
PL (kg/year)	2,232,678.69	4.191	2,232,806.33	4.197	2,235,257.03	4.311
WS (mm/year)	687.36	−5.605	696.67	−4.327	682.74	−6.240
SR (t/year)	5,306,416.54	−6.165	5,278,306.64	−6.662	5,201,598.08	−8.019

and Figure 6.

The results showed that PL increased compared to the baseline scenario, which indicated that the integration of land use and climate change would reduce the phosphorus purification service, and there was not much difference among the three emission scenarios as the increases ranged from 4.191% to 4.311%. Similarly, WS and SR would decrease under the combination of land use and climate changes, but there were obvious changes among the different scenarios, and RCP8.5 had the greatest effect on both services. The change in NL was most evident with reductions of 8.303–8.509%, which indicates that the nitrogen purification service exhibited the largest increase among all services. In general, the differences between the nitrogen purification service and other three services were distinct, and RCP8.5 had the most negative impact on ESs under the integration of land use and climate change.

4. Discussion

4.1. The Combined Effects of Multiple Driving Forces on Ecosystem Services

In this study, we quantified ESs changes compared with the baseline in 2010 to reveal the future effects of land use and climate changes and the combination of the two factors on ESs in the Taihu Basin region. Land use played a more important role than climate change in the nutrient services in the future as the variation increased when land use changes were added to the climate scenarios in the medium term (Table 2 and Table 4). Moreover, the change rates of nutrient services under the influence of land use changes were compared with those under the influence of climate change to further validate the results (Table 2 and Table 3). These suggested that land use change is the primary force impacting nutrient services, indicating the critical role of anthropogenic activities [86,87]. In addition, the changes in the two nutrient services were opposite in the land use alone and combined conditions, which suggested that distinct differentiations existed in the mechanisms that altered the nitrogen and phosphorus services under these conditions. The change of land use type, patterns, and intensity affects ESs supply through the corresponding ecological processes changes [4]. Previous studies have indicated that cultivated land has a primary positive impact on nitrogen loading while construction land mainly causes the increase of phosphorus loading [49,88]. In our study, the total amount and spatial pattern of cultivated land, urban construction land, and rural residential land changed obviously and differently (Table 1, Figure 1 and Figure 2), which led to the different changes of two nutrient services.

For the WS and SR services, the effects of climate change alone and the combined impacts in the future (Table 3 and Table 4) were more remarkable than the effects of land use alone (Table 2), indicating that climate change is the primary driver of the changes in these services. Climate change influences the hydrological process and moisture-energy distribution directly or indirectly [4], resulting in the two services changes. In terms of temporal variation, SR change was similar to the precipitation change trend. It is well known that greater soil loss often occurs in regions with strong precipitation and runoff variability [89], which indicates that SR would be influenced by climate change through erosion and sediment transport processes. There is the same change trend between WS and precipitation when temperatures have slight difference among three RCPs scenarios during 2011–2040, while the trend becomes inconsistent when temperatures have significant difference during 2041–2070 and 2071–2100 (Figure 3 and Figure 5). Despite both the maximums of two climate factors in RCP8.5, WS service is the minimum, while RCP4.5 with the moderate precipitation and temperature has the maximum service, which indicates effective combination of climate factors can promote WS. Generally, it is essential to implement climate change regulation measures to increase WS and SR services. Moreover, although land use played a lesser role in the total amounts of the two services, it amplified and mitigated the negative impacts of climate change on WS and SR, respectively (Table 3 and Table 4), and obviously influenced the spatial changes in WS.

In this context, the studies by Fan [41], Hoyer [90], Chaplin-Kramer [91], and Gardner [92] reported that hydrological ESs were remarkably driven by land use and land cover patterns; in particular, nutrient purification services were more sensitive to the spatial distribution and temporal change in the land use patterns [29]. Meanwhile, studies in the Francolí River basin [25], Upper Merrimack River watershed [37], Tualatin and Yamhill basin [90], and Eastern Tibetan Plateau [4] reached a similar conclusion that climate change would remarkably influence the water yield service in the future, despite the differences in the sediment retention service due to decreased soil erosion as a result of the relatively flat topography. In addition, the studies by Martin [87] and Samal [37] suggested that despite the dominant influence of climate change on water yield, land use could amplify its effects on both the magnitude and timing of water yield and simultaneously mitigate the degradation of certain ESs. These changes occur because land use mediates the hydrological and biogeochemical responses of catchments that are governed by climate drivers through the impacts of vegetation cover, impervious surfaces, and water residence time.

4.2. Policy Implications for Improving Ecosystem Services in the Future

Our study found that land use change caused distinct effects on nutrient purification services while climate change negatively influenced all four services. To mitigate the negative impacts, it is essential to develop effective land use, climate change adaptation, and ecosystem management strategies. Land use changes via the variation in the composition and configuration of landscapes will indirectly affect ESs by altering ecological structures and processes [9], and the water quality service is especially sensitive to land use change [91]. Increases in agricultural land pose a threat to water quality services [29,91]; in particular, increases in agricultural land cause nitrogen enrichment and weaken ecosystem functions related to nitrogen purification [49,50,93]. In recent years, a trend of decreasing cultivated land has been observed, and the simulated results showed that cultivated land would decrease by 26.36% in 2040 (Table 2); therefore, nitrogen loading would also decrease. Although reductions in cultivated land area and intensive agricultural activities are effective methods to enhance the nitrogen purification service, these changes may reduce crop production and threaten regional food security [88]. Therefore, policies related to the water purification service should fully consider the impacts on other ecosystem functions and attempt to balance ecological protection and agricultural development. Improvements to agricultural technology such as soil testing for formulated fertilization, organic farming, crop rotations, and increases in the yield per unit area can be considered useful ways to reduce NL.

Previous studies have demonstrated that phosphorus loading (PL) was closely correlated with urban land and industrial land [49,88,94]. In our study, PL increased with the increase in urban construction land from 16.40% to 22.95%, and the spatial distribution of PL was coincident with the distribution of urban land in our study (Figure 1, Figure 2 and Figure 4). Thus, rapid urbanization may play a prominent role in the degradation of phosphorus purification services, and urbanization policies should be optimized to reduce PL and facilitate phosphorus purification services. First, it is crucially important to promote ecosystem-based green infrastructure, such as woodland patches, grassland greenbelts, and urban wetlands, to strengthen the filtering capability of urban ecosystems. Meanwhile, green infrastructures can regulate surface runoff to enhance water storage and flood regulation services. Second, the spatial patterns of landscape types should be reasonably arranged during the urbanization process based on landscape ecology theory [86], which will help balance the sources and sinks of nutrients at spatial and temporal scales. For instance, embedding vegetation zones in urban construction land and optimizing landscape patches, connectivity, and ecological networks can improve the nutrient-retaining functions of “sink” landscapes as well as ensuring the safety of the urban ecosystem. Furthermore, vegetation buffers around “source” landscapes and watercourses should be preserved through management practices to retain nutrients and sediment [91,95,96].

Climate change exerts substantial effects on WS and SR according to our study, although it is commonly seen as a driving force at the large spatial scale [20,22,23]. As a result, climate change adaptation measures should be implemented to reduce the negative influence. According to the IPCC, there are three main categories of adaptation actions in response to climate change: physical, social, and institutional [97]. Physical actions, such as improving water facilities and increasing plant coverage, are commonly considered effective measures to mitigate the negative impacts of the climate change at the regional watershed scale. As irrigation is the dominant factor that influences soil erosion in cultivated crop areas [98], the construction of reservoirs and irrigation facilities can regulate surface runoff and improve the SR function of agricultural ecosystems in a watershed. Plant coverage can also protect soil from rain wash and enhance the SR service. In addition, plant coverage, especially forest coverage, can influence evapotranspiration and precipitation through microclimatic interactions, capture moisture to be transferred to the soil, and strengthen the water infiltration and retention capacity of soil, which further improve the water yield service [29,31,99]. Despite the limited influence of climate change on local nutrient purification services, the spatial spillover effect of carbon emissions cannot be ignored. Our results indicate that nutrient loadings would increase in the high emission scenario over time. Because hydrologic processes influence the leaching of nutrients from soil to water [41], measures to regulate hydrologic cycles while reducing carbon emissions should be intensified to mitigate the impacts of climate change on multiple ESs.

Overall, the results indicated the importance of the impacts of land use and climate changes on the ESs of a watershed, and it is of great significance to implement comprehensive measures based on these two drivers. Ecosystem-based management (ESM) includes a set of adaptation measures to increase the supply of multiple ESs by optimizing the components, structures, and functions of ecosystems. Ecological zoning management can be considered ESM due to their similar features and integrated measures. In China, policy instruments used in Ecosystem Function Conservation Areas and ecological redline [100] are typical approaches that can be used to strengthen high-quality ecosystem management. Likewise, the Water Ecological and Environmental Function Zoning of 10 major basins is another important plan in China. Currently, Taihu Basin has been the first to apply this zoning practice with three types of management objectives, including water ecology, spatial governance, and species protection. There are three levels; the second level is based on land use, soil type, and climatic indicators, and the third level has 49 zonings to adopt differential policies. However, this zoning management is still at the initial stage because of the lack of adaptation projects. Although measures regarding industrial structures, land use, and ecological restoration have been designed, these measures do not consider the heterogeneity of zoning, and the implementation of policies across multiple scales is insufficient. Thus, managers should consider the different landscape patterns across different spatial scales and a combination of climatic characteristics at multiple temporal scales [9] to develop appropriate measures for the management of watersheds.

4.3. Strengths and Limitations

The creation of effective policies for watershed ecosystems requires a sound knowledge base that can be used to identify the drivers of ESs and understands how the changes in these drivers ultimately impact the provision of multiple ESs. In our study, we explored the relationship between the future changes in two main drivers and multiple ESs. To quantify the future changes in drivers and ESs, we used a series of simulated methods of driving factors and spatially explicit models. As the latest greenhouse gas emission scenario, RCPs were applied to project three possible future climate change scenarios, and the downscaling method was used to make the scenarios suitable for climate change at the regional scale. In addition, we simulated climate changes in three future periods by considering the temporal accumulation and uncertainty of climate factors. Because land use is directly controlled by human activities and exhibits relatively frequent changes, we simulated only one scenario that was most likely to happen in the future. The CA-Markov model was regarded to be suitable to display the spatial heterogeneity of land use and improved the accuracy of the simulation based on embedded land use policies. Furthermore, we elaborated on the policy implications from the impacts of the changes in the drivers and proposed targeted countermeasures to enhance watershed ecosystems.

Although we attempted to improve the reliability of our study, there are still several limitations to be highlighted. First, the selection and calibration of the ESs models included a wide range of uncertainties. For example, the water yield model used in our study mainly quantifies surface water yield. The nutrient models were calibrated based on existing related studies due to a lack of field measurements in most Chinese basins, which led to the potential uncertainties of this result. Besides, we consider the variations in only the model parameters related to land use and climate change, but other parameters might change in the future due to complex physical and social environments. These uncertainties will be reduced with improvements to the mechanism analyses and the further calibration of the model parameters in future research. Furthermore, physical model-based ES assessments rarely consider the socio-political dynamics that play a crucial role in changing land use and emission trajectories less. Therefore, in further research, high spatial and temporal resolution integrated assessment models that include biophysical and socio-economic factors should be developed [101]. Second, the quantification of the drivers of ESs also brings uncertainty. Although we simulated land use based on a range of conditions, urban expansion and agricultural land change may be derived from more complicated driving forces, such as land use management and stakeholder behavior, which is difficult for the CA-Markov model to account for, and an agent-based model [40] and integrated environmental modeling [102] should be used to model such complicated dynamics in future research. Multi-criteria analysis and integrated impact modelling framework also help to track and quantify land policy, environmental payments, human pollution, and other socio-economic drivers [103,104]. The delta method used to establish the climate scenarios at the local scale considers only the differences in the average, maximum, and minimum between the measured data and downscaling data, and this method assumes that the spatial patterns remain unaltered in future climate change scenarios [105], which resulted in certain deviations of the climate factors. Moreover, both data processing and the analysis of the drivers are scale dependent [88]. The main challenges are quantifying the different spatial and temporal scales of drivers and ecosystem processes, developing multiple scenarios to identify cross-scale interactions of drivers and ESs [106], examining hydrologic changes at varying scales relevant to the drivers [87], and exploring different spatial downscaling methods.

In fact, there is a complex interplay among drivers, and this interplay may have a greater influence on ESs than the drivers themselves. Despite the predominant influence of some drivers on ESs, policy decisions regarding other drivers may mitigate or magnify the effect. Further studies should, therefore, include multiple drivers of ESs and investigate the potential interactions of drivers as well as the relationships between the ESs under the multiple drivers. Furthermore, ecosystems may change when the effects of combined drivers exceed the tolerance thresholds of the ecosystem such that the ecosystem is unable to rehabilitate. However, cumulative impacts of drivers along with the resilience of ESs were not analyzed under future scenarios in our study despite being useful at revealing the greatest impacts on ESs [18,20].

5. Conclusions

The scarcity of freshwater has become one of the most important issues in the world. As land use and climate changes and their combined effects have become a pressing issue, our study quantified the future changes in these drivers and the corresponding responses of multiple ecosystem services (ESs) related to water quantity and quality under different scenarios in the Taihu Basin region. The results show that (1) future land use changes had greater effects on nutrient purification services than other services and would enhance the nitrogen service by 8.64% while reducing the phosphorus service by 4.02% compared to the baseline scenario in 2010; (2) climate change would significantly influence the water supply and soil retention services, and these ESs varied under the different future climate scenarios during three periods; (3) when the land use and climate change scenarios were combined, all ESs obviously changed, and the change in the nitrogen service was different from the change in the other three services, and RCP8.5 combined with future land use changes played the most negative roles in the ES changes.

The results indicated the effects of future land use and climate changes on multiple ESs under different scenarios, which helped to explore the predominant driving force of each ES in the future. Based on these results, we recommend that policy makers in the Taihu Basin region optimize the structure and pattern of land use, take climate change adaptation actions, and improve the regulation of ecological zoning management. The asymmetric and nonlinear interactions of biophysical and socioeconomic drivers intensify the changes in the structure, process, and function of watershed ecosystems and further influence ESs and human well-being. Consequently, an important area for future research is quantifying multiple drivers, identifying the uncertainty of the models, and building more scenarios to incorporate the complex drivers for effective and targeted management to ensure the sustainable development of a watershed.