Land use and climate change-based multi-scenario simulation of ecosystem service trade-offs/synergies: A case study of the central Yunnan urban agglomeration, China

TL;DR

This study uses simulations to analyze how land use and climate change affect ecosystem services in the central Yunnan urban area, revealing trade-offs and synergies for sustainable planning.

Contribution

The study introduces a multi-scenario simulation framework combining land use and climate change to assess ecosystem service interactions in a specific urban region.

Findings

Ecosystem services in CYUA show high spatial heterogeneity by 2030, with distinct regional patterns for water yield, soil conservation, habitat quality, and carbon stock.

Trade-offs and synergies between services like WY-HQ and HQ-SR were identified, showing complex spatial relationships.

Across scenarios, consistent patterns of weak synergies and trade-offs were observed, particularly for SR-HQ, SR-CS, and WY-CS.

Abstract

Exploring Land use and climate change-based multi-scenario simulation of ecosystem service trade-offs/synergies is of great importance to regional ecological security and sustainable development. Taking the Central Yunnan Urban Agglomeration (CYUA) as a case study, six different scenarios of LULC-RCP were established to quantitatively assess four key ecosystem services(ESs) of water yield (WY), carbon stock (CS), soil conservation (SR) and habitat quality (HQ) with multiple objective programming and patch-generating land use simulation(MOP-PLUS) and integrated valuation of ecosystem services and tradeoffs (InVEST) models. The ESs were revealed regarding spatio-temporal trade-offs/synergies using Spearman correlation and geographically weighted regression (GWR). It was found that: (1)the ESs in CYUA is characterized with high spatial heterogeneity in 2030; specifically, the distribution of WY and SR was low in the northwestern region and high in the southeastern region, while the distribution of HQ and CS was high in the western region and the periphery, and low in the eastern and central regions; (2) the trade-offs between WY-HQ, and WY-CS, and the synergies between WY-SR, HQ-SR, HQ-CS, HQ-CS, and HQ-SR; (3) under the six different scenarios, the spatial distribution of trade-offs/synergies between the four ESs was consistent: the SR-HQ, SR-CS, and WY-CS showed an overall weak synergistic relationship; the HQ-CS showed an overall weak trade-offs; the HQ-WY, CS-WY showed an overall weak synergistic relationship in the northern and southern areas and an overall weak trade-off relationship in the center. The findings of this study may provide a theoretical foundation for ecosystem management in CYUA and offer technical support for the evaluation of national land space.

Figures

Research location map.

Study technical route.

Land use transfer matrix for CYUA, 2000-2020 (hm2).

Spatial distribution of different land use scenarios for CYUA 2030.

Spatial distribution of ESs in CYUA ESs in 2030 under various scenarios.

Spatial distribution of ES trade-offs/synergies under different scenarios.

Spatial distribution of precipitation and temperature under different climate models.

Equations

Full text

Introduction

Ecosystem services (ESs) are the various direct or indirect products and services that humans obtain from the structure, functions, and processes of ecosystems, and is a bridge connecting human society and natural ecosystems [1]. Thus, the investigation on ESs is of great practical significance for improvement in human quality of life and regional sustainable development [2,3]. Due to the complexity of ecosystem services, spatial heterogeneity, and the diversity of human demands, the interactions between ecosystem services exhibit dynamic changes, forming trade-offs characterized by one increasing while the other decreases, or synergistic effects where both increase or decrease together [4,5]. The ESs changes are directly or indirectly affected by many factors of climate, soil, vegetation, and land type [3], among which, the changes in climatic environments and land-use patterns are two main factors driving changes in ESs and transitions between different types of ESs [6]. As a major driver of climate change, changes in atmospheric composition, e.g., increased concentrations of carbon dioxide, and other greenhouse gases such as nitrogen-containing compounds due to combustion and fertilization, can indirectly contribute to global warming by directly affecting biochemical processes and nutrient cycling in surface organisms [7]. Furthermore, sea level rise and changes in vegetation types due to global warming can in turn have a range of indirect effects on the bio-physical systems of the planet. Studies have demonstrated that climate change mainly featured global warming, has resulted in plenty of environmental issues like water scarcity, degradation of ESs, and increased risk of catastrophes [8], for example, about 60% of ESs have deteriorated because of climate change [9]. Since climate change has already severely affected the environment on which human beings depend and may further deteriorate in the future, it is becoming a significantly challenge for the sustainable development of human societies [10]. In addition, regional sustainable development is being seriously affected by deteriorating regional ESs due to intensification of human activities and ecological degradation. Overall, exploring the influences of land use (LC) and climate change on regional ESs and ESs trade-offs/ synergies to improve ecosystem management has become a global concern.

Land use/cover (LULC) change as a primary factor affecting ecosystem services [11], LULC change mainly affects the provision of ESs with changes in the internal structure and function of regional ecosystems by altering the physical properties of soils, vegetation cover, and LC [4]. Currently, the existing domestic and foreign researches on ESs trade-offs/synergies have been carried out with focus on the spatio-temporal differentiation, drivers, and scale effects [12,13]. For example, Mengba Liu et al argued that at different spatial scales, the impacts of LULC change and its coupling effects on ESs differed significantly among different levels of urban agglomerations(12); Chen et al explored the relationship between intensity of ESs and LULC change of counties in China from a spatial perspective, and the results showed significant spatial dependence and heterogeneity between the intensity of ESs and LULC change [14]. Climate change is a critical factors influencing ESs [15], which directly affects ecosystems with respect to the structure and function by altering the temperature and the frequency and intensity of precipitation, thereby affecting the capacity of providing ESs as well as the interactions and feedback mechanisms within ecosystems [16]. For example, in arid land biodiversity conservation services may be sacrificed in certain areas to maintain water availability [13]. It have been demonstrated the impacts of different climatic factors on ESs, among which the reduced precipitation and increased temperature are the main ones [17]. Consequently, it is vital to evaluate the influences of both climate change and land use on ESs and the ES trade-offs/synergies. Nevertheless, the main focuses of existing is the impacts of single factors (LULC or climate change) on ESs from a historical perspective, and is lack of attention on coupling effects of climate change and LULC, particularly such effects on ESs interactions under multiple scenarios in the future.

The influence of land use change (LUC) on ESs is manifested in the differences in ecosystem structure and function of ecosystems due to land use intensity and the direction of LUC. Currently, many researchers have investigated the influencing mechanisms of LUC on changes in ESs, using the tools like linear regression model, land transfer matrices, and scenario analyses [18–20]. Using a spatial linear regression model, Zhang et al investigated the spatial dependence of ecosystems on LUC in urban areas at different scales, providing practical guidance for urban development planning and environmental protection [21]. Using a transfer matrix in conjunction with a coupled coordination degree model, Sun et al investigated the impacts of LUC on ESs in terms of supply and demand in China and assessed the degree of coupling coordination between LUC and ESs at provincial scale [22]. In the course of these studies, many ES analysis models have emerged, such as the Ecosystem Service Value (ESV) [23,24], the Environmental Quality Index (EQI), and the Integrated Valuation of Ecosystem Services and Trade-offs (InVEST) model [25]. Among these models, the InVEST model enables quantitative analysis on spatial distribution and dynamic change of multiple ecosystem services, and thus has been widely employed for assessment on regional ESs [26]. However, most of the existing studies on the relationship between LUC and ESs were carried out with focus on single ES or overall ESs [27,28], and with less attention to multiple ESs. In addition, the dual pressure of global climate change and rapid urbanization has brought more uncertainty to ESs. Therefore, it is still an urgent issue to investigate the influencing mechanisms of LUC on ESs under multiple scenarios in the future, taking into account regional development needs.

In answer to global climate change, various global future scenario pathways have been devised by the Intergovernmental Panel on Climate Change (IPCC) to envisage a range of climate outcomes for human societies under different development approaches. For example, in the earlier Special report on emissions scenarios (SRES), Representative concentration pathways (RCPs), and Shared socioeconomic pathways (SSPs), a variety of trends that are likely to emerge during the development of societies and ecosystems in the 21st century have been depicted [29]. The establishment of the model provides researchers with criteria for prediction. In this paper, based on the latest IPCC Future Development Scenarios Framework (SSP-RCP), the future LC scenarios were combined with RCP scenarios to construct the MOP-PLUS-RCP coupled model to simulate future changes in ESs under various scenarios considering both climate and land use, which can help to estimate the uncertainty of changes in ESs and the ESs trade-offs/synergies, thereby facilitating the identification of hotspots by policymakers regarding the changes in ESs to formulate more effective environmental responses, and ultimately realizing green and healthy sustainable development.

Plateau urban agglomerations, which are usually featured with unique geography, climate, and ecology, play an vital role in China’s strategic layout and can act important nodes connecting domestic and foreign markets with their special geographical location and resource endowment, facilitating the socio-economic development of the regions along the plateau regions; and therefore they are an important fulcrum for implementing many national strategies in China like the Western Development, the Belt and Road Initiative and so on. At the same time, they are conducive to optimization of the urbanization structure and promotion of balanced regional development. However, they are facing many challenges for their special natural environment, climate change, and intensified human activities, including intricate water-sand dynamics, fragile ecosystems, and imbalance between human and environmental systems. Therefore, predicting changes in ESs and ESs trade-offs/synergies in the future for plateau urban agglomerations based on LC and climate change-based multi-scenario simulations is of great significance for formulating overall regulation policies for sustainable development to achieve sustainable development in plateau regions with harmony between people and the land. Given this, this study uses the Central Yunnan Urban Agglomeration (CYUA) as a case study and employs the constructed MOP-PLUS-RCP coupled model to simulate land use changes (LUC) under three scenarios for 2030: Natural Increase (NIS), Ecological Protection (EPS), and Economic Development (EDS). Additionally, six different LULC-RCP scenarios were set up in conjunction with two climate scenarios (RCP4.5 and RCP8.5). Secondly, the InVEST model was used to quantitatively assess future changes in four key ecosystem services (ESs): water yield (WY), carbon storage (CS), soil retention (SR), and habitat quality (HQ). Thirdly, Spearman correlation analysis and geographically weighted regression (GWR) methods were employed to deeply analyze the spatio-temporal trade-offs and synergies between various ecosystem services. Finally, the main development pathways and management policies for future ecosystem management were proposed.

The study aimed (1) to simulate the changes in and spatial distribution of land use demand in 2030 under various scenarios by coupling the MOP-PLUS model, (2) to evaluate the spatio-temporal changes of four key ESs under six LULC-RCP scenarios with the InVEST model in combination with future LC and climate scenarios, (3) to investigate the trade-offs/synergies between certain ESs under various LULC-RCP scenarios using Spearman correlation and GWR, (4) to propose main development pathways and management policies for ecosystem management in the future in CYUA, according to the changes in ESs and ESs trade-offs/synergies under certain scenarios. The study is to creatively consider both climate change and LC, with focus on the changes in multiple ESs and the trade-offs/synergies between them, to offer policy recommendations for determining preferential development pathways and improving the overall ecosystem quality. The results of the study offer new insights, theoretical basis, and strategic directions, which may contribute to the solution to climate change and the achievement of sustainable development in CYUA.

Study area and data source

Study area

In the early 20th century, geographers and urban scholars began to focus on the interactions between cities, introducing concepts such as urban clusters to explain the spatial distribution patterns of cities and their interactions. In China, with the deepening of urbanization and regional economic development, urban agglomerations have become an important research field and practical issue for promoting national and local economic growth. CYUA is the city cluster in the most developed area of Yunnan Province (Fig 1), with a total area of 111,400 km^2^. It is one of the 19 city clusters in China with that the state focuses on cultivating, and the core area of productivity layout and production development in Yunnan Province. CYUA is dominated by mountains and basins, and belongs to the high altitude areas at low latitude. The region is rich in biological resources, with a forest cover of more than 50%, vegetation primarily dominated by subtropical evergreen broadleaf forests and scrublands & grasslands, and well-developed river systems, making it an important regional ecological barrier. Recently, as a result of accelerating urbanization, the demand on natural resources is growing, intensifying the conflict between environmental protection and urbanization as well as human activities, and the sustainable development in the region are facing multiple pressures and severe tests in terms of regional land and the ecological environment.

Research location map.

Data sources and data processing

This study utilizes land use maps, digital elevation models, soil data, climate data, and socio-economic data. Table 1 summarizes the resolution and sources of each data type. Data preprocessing was carried out using ArcGIS software. Based on the 25 land use categories available in the data, we reclassified the land use data into cropland, forest, grassland, water bodies, built-up areas, and unused land. All data underwent preprocessing steps such as projection transformation, clipping, mask extraction, and resampling. ArcGIS was used to unify the coordinate system of all raster data, and the spatial resolution was standardized to 30 meters.

Table 1: Data sources.

Research methodology

Projections of future climate and land use change scenarios

This study takes CYUA as a case study and constructs the MOP-PLUS-RCP coupled model to predict different climate and land use scenarios for 2030. The Integrated Valuation of ESs and InVEST model was used to quantitatively assess the future changes in four key ESs: water yield (WY), carbon storage (CS), soil retention (SR), and habitat quality (HQ). Spearman’s correlation analysis and geographically weighted regression (GWR) were applied to explore the spatiotemporal trade-offs and synergistic effects among different ESs. Finally, the study proposes key development pathways and management strategies for future ecosystem management. The overall research framework is shown in Fig 2.

Study technical route.

Climate data and emission scenarios.

To predict the distribution of future climate forcings, this study uses two Shared Socioeconomic Pathways (SSPs) from the latest CMIP6 updates: RCP 4.5 and RCP 8.5. RCPs are quantitative descriptions of alternative environmental development pathways for the 21st century, combining climate change projections with assumptions about mitigation or adaptation policies to generate comprehensive climate change scenarios. The RCP 4.5 scenario (4.5 W/m²) represents a moderate path for greenhouse gas emissions. In this scenario, environmental systems face some degree of degradation, and climate protection measures are assumed to be implemented [29]. The RCP 8.5 scenario (8.5 W/m²), on the other hand, represents the extreme case of rapid socioeconomic development, leading to the highest path for future greenhouse gas emissions [30]. In this scenario, global markets become increasingly integrated, fostering innovation and technological progress. However, socioeconomic development heavily relies on intensified exploitation of fossil fuel resources, with a higher proportion of coal, and an energy-intensive lifestyle continues globally [31]. This study selects RCP 4.5 and RCP 8.5 because they represent the moderate and worst-case scenarios for greenhouse gas emissions, respectively.

Land use change scenarios.

MOP is an open and flexible method, which mainly consists of three parts: decision variables, objective function, and constraints, and takes into account the manager’s expectations of relevant planning by setting proper objective optimization functions and constraints [32]. To investigate the trade-offs/synergies between ESs on the basis of simulated LUC under different scenarios, three LC scenarios of NIS, EPS and EDS were constructed herein in the light of the characteristics of LUC and compliance to the requirements of macro-policy regulation.

NIS: in this this scenario, the demand on LC by 2030 followed the same trend as the historical data, is not affected by policies, and serves as a reference baseline for other scenarios. The amount of various types of LC in 2030 under this scenario was calculated using the LUC matrix (2010–2020), which was generated by Markov model, and with various types of LC in 2020 as the baseline.EPS: this scenario indicates that priority is given to the protection of ecological land (woodlands, grasslands, waters). In this scenario, the ecological benefits provided by each type of LC are maximized through planning and management, and the function for assessing ecological benefits is expressed as follows.

$$[eqn]$$

In equation (1): f(x) denotes the total ecological benefit of all types of LC, ESV_i_ denotes the ecological benefit per hectare of the type i land use (in million CNY), and x_i_ denotes the area of the type i LC. In this paper, the economic value was obtained by correcting the economic value generated by the yield per unit area of the main crops in CYUA, with reference to the revised Equivalent Value per Unit Area of ESs in China by Xie et.al as the basis [33]. The objective optimization function of the EPS is: $[eqn]$ .

EDS: this scenario represents a great effort on economic development effort, so as to maximize the economic benefits to humans from various types of land, and the function for assessing economic benefits is given in equation (2) below.

$$[eqn]$$

In equation (2), f(x) denotes the total economic benefit, EC_i_ denotes the economic benefit per hectare of type i LC (in million CNY), and x_i_ denotes the area of type i LC. The average land economic value is sourced from the Yunnan Statistical Yearbook, where the total value of output for forestry, animal husbandry, and fishery corresponds to the economic benefits achieved by woodland, cropland, grassland, and waters, respectively. The GDP of the secondary and tertiary sectors corresponds to the economic benefits of built-up land, and the unused land was set based on the previous study (Ling et al., 2024). The objective optimization function of the EDS scenario is: $[eqn]$ .

The constraints are as follows.

1)Total land area constraint: It is necessary to ensure that the total area of all land use types remains unchanged. As shown in equation (3) below.

$$[eqn]$$

2)Total population constraint(P): this constraint ensures that the total population living on the land should not exceed the land’s carrying capacity for human activities. The upper limit of population was set according to the policy laid down in the 《Central Yunnan Urban Agglomeration Development Plan》(https://www.yn.gov.cn/), and the existing population is the lower limit. Therefore, the mathematical expression for the constraint on the total population can be written as equation (4).

$$[eqn]$$

In the formula, P represents the total population of the study area, measured in ten thousand people.

3)Cropland area constraint: the food yield of farmland should be at least the food demand of the population. Herein, the per capita food demand was assumed as 141.2 kg; the food self-sufficiency rate is 114.11 percent; the food production per hectare of the cropland was 4316.24 kg; the proportion of crops planted was 67 percent; and the multiple cropping index is 1.11. The cropland constraint can be written as follows:

$$[eqn]$$

4)Forest area constraint: The area of forest land in the CYUC showed a downtrend from 2000 to 2020, which was expected to continue in the future. According to the ‘Yunnan Province Ecological Civilisation Construction Leading Plan (2021-2025)’ and the ‘Development Plan for the Central Yunnan Urban Agglomeration (CYUC)’, the local government scheduled to implement more ecological protection plans. Herein, using the woodland area in 2020 as the baseline, the land area (2020) ±1% was determined as the upper and lower bounds, respectively; and then the constraint can be written as:

$$[eqn]$$

5)Grassland area constraint: a downtrend of grassland area was identified during the study, so using the grassland area projected by the NIS scenario as a baseline, the projected area ±1% was determined as the upper and lower bounds, respectively; and then the constraint can be written as:

$$[eqn]$$

6)Water area constraint: Large areas of water in CYUA have been reclaimed for urbanization and crop production. To prevent further shrinkage of lake area, the local government has carried out the Dianzhong Diversion Project for water resources protection. As a result, the water area will be on an upward trend in the future. Therefore, the water area in 2020 was set as the lower bound of constraint, and the water area in 2030, which was calculated according to the natural growth rate of the water area from 2015 to 2020, was considered as the upper bound of the constraint, then the constraint can be written as:

$$[eqn]$$

7)Built-up land constraint: the built-up land should satisfy the population demand and the per capita building land in CYUA was 0.0109 ha. Therefore, the constraint can be written as:

$$[eqn]$$

However, the built-up area was projected to be 332,112.951 hectares after 2020 and 367,072.209 in 2030, based on the current trend of growth in built-up land in CYUA. Then the constraint can be detailed as:

$$[eqn]$$

8)Unused land area constraint: since CYUA has been reclaiming unused land for long time due to socio-economic development, it is expected a smaller area of unused land in 2030 than that of 16,049.16 ha in 2020. However, since CYUA is a plateau zone, the utilization rate of unused land is low. Therefore, the unutilized land will not be less than 15849.3 ha. Then, the constraint can be written as:

$$[eqn]$$

9)Other constraints: in view of the ecological protection measures proposed in the ‘Yunnan Province Ecological Civilisation Construction Leading Plan (2021-2025)’, the study concluded that the area of ecological land (woodland land, grassland, and water) in CYUA will be at least 75% of the total area by 2030. Therefore, the constraints can also be written as:

$$[eqn]$$

PLUS model.

As a raster data-based meta cellular automata (CA) model, the PLUS model can be used for simulation of patch-scale change in LULC [34]. The model integrates a land expansion analysis strategy-based rule mining framework (LEAS) and a multi-type random patch seeding-based CA model (CARS). Using the random forest classification (RFC) method, the proper probability for various types of LULC expansion can be obtained by mining the driving factors. It consists of two modules: (1) a rule mining framework based on Land Expansion Analysis Strategy (LEAS); (2) a Cellular Automata (CA) model based on multi-type random patch seeding (CARS). The steps are mainly as follows.

By using LEAS, the change patches of different land use types over a period can be analyzed to obtain the transition rules for all land use types. The transition rules for each land use type are mined and transformed into a binary classification problem, which allows for the calculation of the change probability and inertia probability for each land use type. Additionally, LEAS can employ the Random Forest Classification (RFC) algorithm to explore the relationship between the growth of different land use types and multiple driving factors.The CARS module is a Cellular Automata (CA) model based on multi-type random patch seeding. The CA model is a scenario-driven land use simulation tool that integrates both “top-down” (i.e., overall land use demand) and “bottom-up” (i.e., regional land use competition) effects. During the simulation, land use demand influences land use competition through adaptive coefficients, thereby driving the land use amount to meet future demand (obtained through the MOP model in this study).

LULC-RCP multi-scenario model.

According to the needs of the study, six modelling scenarios, as shown in Table 2, were established by coupling the future climate models (RCP4.5, RCP8.5) and the three land use scenario data (NIS, EPS, EDS) in 2030. Among them: in Scenario 1, the data on LC under NIS scenario and the climate of RCP4.5 scenario in 2030were used as input data; in Scenario 2, the data on LC under NIS scenario and the climate of RCP8.5 scenario in 2030 were used as input data; in Scenario 3, the data on LC under EPS scenario and the climate of RCP4.5 scenario in 2030 were used as input data; in Scenario 4, the data on LC under EPS scenario and the climate of RCP8.5 scenario in 2030 were used as input data; in Scenario 5, the data on LC under EDS scenario and the climate of RCP4.5 scenario in 2030 were used as input data; and in Scenario 6, the data on LC under EDS scenario and the climate of RCP8.5 scenario in 2030 were used as input data.

Table 2: Future LULC-RCP multi-scenario models.

Quantifying ESs

1. WY.

The WY of the waters characterizes the regional water supply capacity, and during the experiment, in the InVEST model, the WY module was used for WY evaluation. The WY module is on the basis of the law of water balance, taking into account the factors of vegetation cover, climate, and soil texture, etc., and the primary formula is as follows:

$$[eqn]$$

Where: $[eqn]$ denotes the W of each raster x (mm); $[eqn]$ and $[eqn]$ denote the actual annual evapotranspiration and annual precipitation of each raster x (mm), respectively [17].

2. CS.

In the InVEST model, the CS module divides the CS of EC into four parts: above-ground CS, below-ground CS, soil CS and dead organic matter CS. Using the module, the CS in the watershed was calculated herein with the following equation:

$$[eqn]$$

Where: $[eqn]$ denotes the total CS (t/hm^2^); $[eqn]$ denotes the above-ground CS (t/hm^2^); $[eqn]$ denotes the below-ground CS (t/hm^2^); $[eqn]$ denotes the soil CS (t/hm^2^); and C_dead__ denotes the dead organic matter CS (t/hm^2^). Herein, the data on carbon density for the various types of LC was sourced from the National Ecosystem Science Data Centre and certain literature [28,35,36]. The carbon density data were sourced from studies on various regions in China, thus the selected literatures should have a study area as close to or similar to that in this study as possible to prevent significant error due to data gap. In addition, carbon density varies depending on many factors like climate, soil properties and LC. Therefore, appropriate corrections should be made in accordance with the types of LC and the climate characteristics in different periods in CYUA. It has been reported a positive correlation of carbon density with annual precipitation and a weak correlation between carbon density and average annual temperature [37]. The specific calculation process refers to equations (3)–(9) in the literature [38], and the carbon densities in each period are shown in Annex 1.

3. HQ.

In the InVEST model, the HQ module has been employed for calculating the HQ of upper basin of the Yangtze River, using the following formula. In this paper, cropland and built-up land were used as habitat stress factors, and the values of the required parameters were assigned with reference to the existing literature [16].

$$[eqn]$$ $$[eqn]$$

where $[eqn]$ denotes the HQ of raster x; $[eqn]$ denotes the total threat level of raster x; K and Z denote scaling factors; $[eqn]$ denotes the HQ of LULC type j; R denotes the coercion factor; $[eqn]$ denotes the number of rasters occupied by the stress factor r; $[eqn]$ denotes the weight of the stress factor, within a range of 0–1; $[eqn]$ denotes the value of the coercion factor for raster y; $[eqn]$ denotes the coercion ( $[eqn]$ ) of the raster y to the level of stress for raster x; $[eqn]$ denotes the level of accessibility for raster x; and $[eqn]$ denotes the sensitivity of type j habitat to the stress factor r.

4. SC.

Herein, the erosion model of Revised Universal Soil Loss Equation was used to estimate soil conservation in CYUC, using the following formula [39]:

$$[eqn]$$

Where: A denotes the soil conservation per unit area; R denotes the factor of rainfall erosivity; K denotes the factor of soil erodibility; LS denotes the factor of slope length; C denotes the factor of vegetation cover and management; P denotes the factor of soil and water conservation measures.

ESs trade-offs/ synergies

Spearman correlation analysis.

To detect the correlation between the ESs on the time scale, Spearman’s correlation analysis and significance test were performed based on the R platform using the ‘Hmisc’ package [40]. Positive and negative correlations indicate synergies and trade-offs between ESs, respectively; and the larger the absolute value of the correlation coefficient, the greater the relationship. The calculation formula is as follows [41]:

$$[eqn]$$

Where, $[eqn]$ represents the Spearman rank correlation coefficient, which indicates the correlation between variables $[eqn]$ and $[eqn]$ , with a value range of [−1, 1]. $[eqn]$ , represents positive correlation, $[eqn]$ , represents negative correlation, $[eqn]$ , represents no correlation. $[eqn]$ and $[eqn]$ represent the values of variables $[eqn]$ and $[eqn]$ for the i-th data point, where i = 1,2,…,n. $[eqn]$ and $[eqn]$ are the sample means of variables $[eqn]$ and $[eqn]$ , respectively.

Geographically Weighted Regression (GWR).

The ESs trade-offs/synergies were quantified spatially using GWR, and the formula is as follows [42]:

$$[eqn]$$

Where, $[eqn]$ denotes the geographic coordinates of the i^th^ raster; $[eqn]$ denotes the measured value of i^th^ dependent variable Y at the spatial location; $[eqn]$ denotes the measured value of the independent variable dataset $[eqn]$ at the spatial location; k denotes the total number of independent variables; $[eqn]$ denotes the constant at $[eqn]$ ; $[eqn]$ denotes the value of continuous function at the i^th^ grid point; and $[eqn]$ denotes the random error.

Results

Multi-scenario simulation of LC

MOP-based forecasts of LC demand.

The analysis of LC data in CYUA between 2000 and 2020 indicated that the general trend of LUC in CYUA during the 20 years was a reduction in the area of cropland, woodland and grasslands, and a rise in the area of waters and built-up land, and a stable area of unused land (Supplementary Material, Appendix1 S1 Fig and Fig 3). According to the calculations following the aforementioned equations (1)–(13), the area of each type of land use under various scenarios in 2030 was obtained and compared with those in 2020, and it was found that: cropland decreases under the NIS and EPS scenarios, and increased in the EDS scenario; woodland land decreased under the NIS and EDS scenarios, and increased under the EPS scenario; grassland decreased to varying degrees under all three scenarios, with the most obvious decrease up to 69799.8 ha under the EDS scenario; waters increases in all three scenarios; and built-up land showed a similar trend to that of waters, with highest increase of 100826.6 ha under the EDS scenario, followed by an increase of 83346.93 ha under the NIS scenario. All the data is detailed in Table 3.

Table 3: Area of different types of land use (hm2) in CYUA in 2030 under different scenarios.

Land use transfer matrix for CYUA, 2000-2020 (hm2).

Multi-scenario simulation of LC based on MOP-PLUS model.

With various data from multiple sources, including data on LC, topography and socio-economy in CYUC in 2000, 2010 and 2020, the simulation of land use in CYUA in 2020 was performed by the PLUS model, and then a comparison between the simulation data (Supplementary Material, Appendix1 S2b Fig.)and the actual data in 2020(Supplementary Material, Appendix1 S2a Fig.) was carried out. The simulation data showed an overall accuracy of 91.97% and a kappa coefficient of 0.876, indicating that the model meets the accuracy requirements of this study.

Based on the land use maps from 2000 to 2020, we predicted the land use spatial distribution for 2030 under the NIS, EPS, and EDS scenarios (Fig 4). It is found that the dominated type of land use in CYUA under all the three scenarios is woodland, which is largely distributed in the highland mountainous regions in the western and southern areas, and is an important ecological barrier in Yunnan Province, with a proportion of approximately 49.08% of the total area; followed by grassland with an area proportion of approximately 26.84%, which is largely distributed in Kunming City, Qujing City, and Honghe Prefecture; cropland is largely distributed in blocks in the plains of CYUA, with an area proportion of approximately 20.23%; waters are mainly dominated by the six highland lakes distributed in the center of CYUA; built-up land is mainly centered in Kunming City, and the rest is distributed in the form of blocks in the urban centres and their surroundings; and the unused land accounts for a very small proportion of 0.2%, which is mainly dominated by the high altitude snow-covered mountains. Furthermore, the simulations of LC under the various scenarios indicated an increase in built-up land and waters in 2030 in comparison to that in 2020; specifically, the expansion of built-up land is mainly due to encroaching on cropland and forests, and the expansion of waters is due to spreading centered around the lakes.

Spatial distribution of different land use scenarios for CYUA 2030.

Changes in ESs and trade-off/ synergies for multi-scenario simulations

Analysis of changes in CYUA ESs for multi-scenario simulations.

The MOP-PLUS-RCP and InVEST models were coupled to predict the future (2030) changes of WY, HQ, SR, and CS services in the CYUA under six different scenarios, and the results are provided in Table 4 and Fig 5 and described as follows:

Table 4: Number of simulated ESs in CYUA in 2030under various scenarios.

Spatial distribution of ESs in CYUA ESs in 2030 under various scenarios.

(1)WY and characteristics of its spatial distribution: under different scenarios, Scenario 6 was projected to have the highest WY (5312.13 m^3^/ha), followed by Scenario 2 (5300.90 m^3^/ha) and Scenario 4 (5288.22 m^3^/ha). For climate scenarios, the WY of RCP8.5 is higher than that of RCP4.5, which may be related to the high precipitation and low temperature of RCP8.5. For land use scenarios, the WY of EDS scenario was larger than that of the EPS and NIS2 scenarios, and the main explanation to this is the increase in cropland under the EDS scenario, which mitigates the impacts of reduced woodlands and grasslands. With respect to spatial distribution, it was roughly the same under the same climate scenario. For RCP4.5 (in Fig 5 S1, S3 and S5), the WY showed a trapezoidal distribution, with gradual increase from north to south; for RCP8.5 (in Fig 5 S2, S4 and S6), the distribution of WY was low in the northwestern region, and high in the eastern and southern regions. It should be noted that the lower reaches of the Yuan River in the south were always the main water-producing area in all scenarios.(2)HQ and characteristics of its spatial distribution: From the average HQ of various scenarios, the effect of land use on HQ was greater than that of climate. Therefore, HQ showed no significant change under different climate scenarios. While for HQ under LC scenarios, it was ranked from high to low as follows: EPS (0.610)> NIS (0.606)> EDS (0.604), all of which are of high HQ, and this is closely correlated with the high forest coverage rate of the region. Therefore, due to more attention paid to ecological protection, which resulted in increased woodland, the EPS scenario showed the best HQ. With respect to spatial distribution, all six scenarios showed a same pattern of high in the western region and low in the centre. Changes in HQ were largely concentrated in the center, where is densely populated and highly urbanised, indicating that HQ was directly affected by the increase in built-up land.(3)SR and characteristics of its spatial distribution: among the various scenarios, S4 showed the best SR, followed by S2 and S6. Compared with that in 2020, the SR of all the six scenarios has increased, with the EPS scenario showing the most significant increase. With respect to spatial distribution, the high SR was mainly located in the lower reaches of the Yuan River in the south, while the rest of the study area showed different low values, especially in the central area, with a downtrend from west to east. This indicates that the high vegetation cover in the south plays a positive role in soil conservation, while in the west, due to the high altitude difference and dense steep slopes, the soil and water conservation capacity is compromised even if there is a high vegetation cover. Although the terrain in the central region is relatively flat, most of the region is an urban area under development, so the soil conservation capacity is low, indicating that continued urban expansion has led to encroachment on ecological land.(4)CS and characteristics of its spatial distribution: the order of CS under different scenarios is as follows: S4 > S2 > S6 > S3 > S1 > S5. It can be seen that the CS capacity under the RCP8.5 is significantly greater than that under the RCP4.5, and that the EPS scenario has the most prominent CS capacity among the various LC scenarios. From the perspective of their spatial distributions, the CS capacity and HQ of S6 showed a highly consistent overall pattern of high in the west and low in the centre. In addition, the CS reduction mainly occurs in the periphery of the central cities. The above findings suggest that urban development in the future will occupy a large number of lands with high CS, such as agricultural land and forests, resulting in decreased CS in the region.

Changes in ES trade-offs/synergies under different scenarios.

Using the Spearman correlation, the trade-offs/synergies between ESs under various scenarios was analyzed to explore the changes in relationship between ESs under various scenarios in 2020 and 2030. As shown in Table 5, the results indicated a positive correlation between WY-SR, HQ-SR, HQ-CS and SR-CS, reflecting a synergistic relationship, and a negative correlation between WY-HQ and WY-CS, indicating a trade-off relationship, in both 2020 and 2030 across different scenarios. Compared to 2020, under the RCP4.5 climate model (S1, S3 and S5), the trade-off between WY-HQ and WY-CS in 2030 was weakened, while the synergistic relationship between WY-SR was strengthened, particularly in S5, where the increase was approximately 0.22. The synergistic effects between CS-SR and CS-HQ are weakening and increasing, respectively; under the RCP8.5 climate model (S2, S4 and S6), the synergistic relationship between CS-SR and CS-HQ were weakened and increased, respectively, which are in exactly the opposite way to that of the RCP4.5 model. It is noteworthy that the HQ-SR synergies showed different degrees of weakening in all six future scenarios. Under the same climate scenarios, the trade-offs between WY-HQ and WY-CS and the synergies between WY-SR, HQ-CS, and SR-CS were higher in the EDS scenario than in the EPS and NIS scenarios, and the HQ-SR synergy was highest in the EPS scenario.

Table 5: ES trade-offs/synergies under various scenarios.

For the purpose of further elucidating the significant spatial heterogeneity of ESs trade-offs/synergies due to geographic differences in highland mountainous regions, the spatial layout of such trade-offs/synergies was obtained at the raster scale using GWR analysis. As presented in Fig 6, the results indicated that the trade-offs/synergies between four pairs of ESs were consistent in spatial distribution under all the six scenarios; SR-HQ spatial effects showed a weak synergy in general, but showed a weak trade-off in the highly urbanised region in the central area; SR-WY showed a synergy in areas distributed in blocks in Qujing City, Chuxiong City, and Honghe Prefecture, and showed a trade-off relationship in areas mainly distributed in Yuxi City; SR-CS showed a weak synergistic relationship in the whole area excepting the southwestern region, in where it showed a weak trade-off relationship; HQ-WY showed a strong synergy in areas largely located in the north and the south, and a weak trade-off relationship in areas mainly distributed in the west and the northeast and particularly evident in the central region; HQ-CS showed a weak synergy in the western and northeastern regions, and a strong synergistic relationship in the rest; WY-CS showed a weak synergy in the northern and northeastern regions and a relatively strong synergy in the south of Honghe Prefecture, and a weak trade-off in the southwest of Yuxi City, Kunming City, and the east of Qujing City.

Spatial distribution of ES trade-offs/synergies under different scenarios.

Discussion

The purpose of this study was to systematically understand the influencing mechanisms of climate and land use changes on ecosystem services (ESs). This includes analyzing their future trajectories and the trade-offs/synergies between them. To achieve this, a coupled MOP-PLUS-RCP model was established. Using this model, land use change (LUC) in CYUA for 2030 was simulated under three different scenarios: NIS, EPS, and EDS. Next, the simulations were combined with two future climate models, RCP4.5 and RCP8.5. This led to the creation of six LULC-RCP scenarios. Under these scenarios, the future changes in four key ESs (WY, CS, SR, HQ) were quantitatively assessed using the InVEST model. Additionally, spatio-temporal trade-offs and synergies in ESs were explored through Spearman correlation analyses and Geographically Weighted Regression (GWR) at the grid scale.

Multi-scenario simulation of LC

The assignment of LC planning is a stepwise process from top to bottom, which is manifested by focusing on the scale and proportion of LC and neglecting the overall spatial layout and regional integration, thereby resulting in generally the problem of spatial conflict between the structure-quantity synergy of top-down LC and the bottom-up layout of LC, when simulating the future LC [43]. In this study, the MOP-PLUS coupled model was constructed for simulating the quantity-structure and spatial distribution of LC of multiple scenarios in the future, and according to the principle of maximising economic and ecological benefits, EDS and EPS scenarios were set up; moreover, NIS scenario was set up as a reference scenario based on the historical trend of development. The MOP algorithm can predict the future quantity demand of LC from top to bottom according to the governmental policy requirement and the actual situation of CYUA. In addition, the PLUS model combines the RF and the bottom-up CA allocation mechanism to reasonably allocate the quantities of LC under each scenario, and to simulate the types of LC under various scenarios. Therefore, the combination of the MOP algorithm and the PLUS model can solve the synergistic problem in simulating quantities and spatial patterns of LC under various scenarios, which lays the foundation for subsequent researches.

Influences of changes in LC and climate on ESs and trade-offs/synergies between them

The LUC and climate have been recognised as critical influencing factors on changes in ESs, and such changes will increasingly affect the dynamics of ESs in the future [44]. Based on the spatio-temporal distributions of ESs under two climate models (Fig 7) and six scenarios(Fig 5), it was found that the average values of WY, SR, and CS in the RCP8.5 climate model were higher than those in the RCP4.5 climate model; whereas, the precipitation and air temperature gradually increased spatially from northwest to southeast, as well as the WY and SR. The explanation for this is that regions with higher precipitation and air temperature typically have denser vegetation cover, which more effectively reduces the risk of soil erosion and directly influences the recharge of surface runoff and groundwater flow [45]. Moreover, the effect of precipitation on WY may be partially offset by a large water consumption by high vegetation cover due to transpiration and by high actual evapotranspiration attributed to water retention and evapotranspiration through the forest canopy [46]. Thus HQ-SR, WY-SR showed synergistic effects. According to previous studies, there is a positive correlation of carbon density with precipitation and a weak correlation of carbon density with air temperature [37]. Since in the two climate models of RCP4.5 and RCP8.5, the average annual precipitation was 851.55 mm and 1080.60 mm, respectively, and the average annual temperature was 16.55°C and 16.47°C, respectively, the carbon density of various types of land use in the RCP8.5 is remarkably higher than that in the RCP4.5, This result is consistent with the findings of Zhang et al. in the Xijiang River Basin [31]. The higher the vegetation cover, the higher the CS, and in addition, a higher vegetation generally means a richer biodiversity and higher soil nutrients, and Similar to the findings of Wang et al [47]. Therefore, HQ-CS and SR-CS showed a relatively high synergistic relationship. However, precipitation promotes both the growth of vegetation and continuous loss of nutrients from the soil, especially in the northwestern part of the region, which has large elevation differences and steep slopes, leading to serious soil nutrient runoff despite high vegetation cover, while the eastern part belongs to the typical karstic ground, with sparse vegetation and serious rainfall erosion [48]. Therefore, WY-HQ and WY-CS showed a weak negative correlation in these regions.

Spatial distribution of precipitation and temperature under different climate models.

In term of LC, under the same climatic conditions, the main changes of the four ESs in the NIS, EPS, and EDS scenarios were concentrated in Kunming City and its periphery, such as Yuxi City, Qujing City, and Chuxiong City, etc., which are regions where the increasing population and accelerated urbanisation have resulted in the sustained expansion of built-up land into the surrounding croplands, woodland, grasslands, and waters, especially in the EDS scenario, which may bring a number of ecological issues, such as the reduction of biodiversity and carbon sink capacity due to the reduction of forests and grasslands; the weakened ability of urban ecosystems to purify air and regulate the hydrological cycle due to reduction of waters [49]. The study also found that areas with expanding built-up land tend to be areas with more severe trade-offs between ESs. These findings can help decision makers to implement ecological projects in the future.

This study can strengthen stakeholders’ understanding on the mutual effects between ecosystem performance indicators under various scenarios. For future policies on land use regulation, it is necessary to take proactive interventions to steer land use towards sustainability. The main recommendations are as follows: ① in term of climate, strengthening climate monitoring and developing appropriate management policies for different patterns of global climate change with use of prediction models; ② in terms of ecological protection, enforcing the closure of high energy-consuming and high-polluting factories, prohibiting grazing, and preventing deforestation activities of land reclamation and mining, etc. In addition, it is recommended to implement actively restorations to general ecoregions such as reforestation and afforestation, thereby strengthen ecological protection.

Shortcomings and prospects

In terms of land prediction, although the MOP method can integrate the policy planning and the actual situation of CYUA into the simulation of future land use, it is unable to constrain the natural geographic conditions during the LUC, such as precipitation, temperature and other climatic conditions that may change the conversion of land use; in the computation and selection of biophysical parameters for various modules of the InVEST model, some of the required data and parameters are sourced from existing literature and may differ from the actual situation of CYUA, leading to difference between the results and the actual situation or existing research results. Although different scenarios were set up herin to explore the influences of land use and climate in the future on ESs, the probability of future development under various scenarios of climate and land use could not be fully explained, and it is necessary to investigate the dynamic changes of trade-offs and synergies between ESs and the relevant driving mechanisms under various climatic scenarios and at different scales in the future research, so as to strengthen the comprehensive evaluation on the trade-offs and synergies between ESs in CYUA, and to offer proposals on ES zoning and regulation of anthropogenic activities.

Conclusion

(1)In simulations under all the scenarios, the area of both waters and built-up land increased in the CYUA. For waters, the highest increase rate of area was in the EPS scenario, followed by the NIS and SD scenarios; for built-up land, the highest increase rate was in the EDS scenario, followed by the NIS scenario.(2)The spatial heterogeneity of ESs was strong in CYUA. Specifically, the spatial distribution of WY and SR was in general low in the northwestern region and high in the southeastern region; and for HQ and CS, it was high in general in the western and surrounding regions and low in the eastern and central regions. In all the land use scenarios, the ESs of RCP8.5 climate model were higher than those of RCP4.5 climate model. While under the two different climate models, the order of WY was EDS > NIS > EPS, and the order of HQ, SR and CS was EPS > NIS > EDS.(3)The changes in ES trade-offs/synergies varied from different scenarios. Among them, WY-SR, HQ-SR, HQ-CS, SR-CS presented a positive correlation, and WY-HQ and WY-CS presented a negative correlation. Compared with 2020, WY-HQ, WY-CS, HQ-SR, and SR-CS in 2030 under Scenarios 1, 3, and 5 showed a weakening trend, and WY-SR showed a strengthening trend; while the behaviors under Scenarios 2, 4, and 6 were exactly the opposite.(4)Under all the six different scenarios, the spatial distribution of ES trade-offs/synergies showed no difference for the four ESs. The spatial effects of SR-HQ, SR-CS, and WY-CS were generally weak synergistic relationships; the spatial effects of HQ-CS were generally weak trade-offs; and the spatial effects of HQ-WY and CS-WY showed weak synergies in the northern and southern region of CYUA, and weak synergies in the central region.

Supporting information

S1 FigLand use data from 2000 to 2020.(TIF)

S2 FigSimulated land use data for 2020.(a) versus actual land use data (b).(TIF)