Study on the Interaction Effects of Landscape Pattern on the Synergistic Trade-Offs of Ecosystem Services Based on Multi-Model Fusion: A Case Study of Chengdu-Chongqing Economic Circle

Abstract

A deep understanding of the spatiotemporal changes in ecosystem services (ESs) under the influence of urbanisation, as well as clarifying the trade-offs and synergies between different services and their driving factors, is crucial for sustainable regional development and the formulation of rational urban expansion policies. Dramatic changes in landscape patterns, driven by the interplay of human activities and natural processes, can have profound effects on regional ESs. Existing research primarily discusses the synergistic trade-offs between ESs, with less focus on the interactions among different landscape patterns and the synergies among ESs. In the present study, we established a multi-model fusion method for ES analysis (PLUS-InVEST-Trade-offs/Synergies-Geographical Detectors (GDs)) to explore the synergistic trade-offs of ESs and their driving factors in the Chengdu-Chongqing Economic Circle from an urban agglomeration perspective. Our findings indicated the following. (1) The distribution of synergistic/trade-offs relationships among ESs exhibited spatial variability. The varying responses of different urban clusters to these policies, along with their unique topography and landforms, are the reasons behind the differences in the synergistic/trade-offs relationships of ESs among these urban clusters. (2) In the Chengdu-Chongqing Economic Circle, the driving factors of the synergistic/trade-offs effects among ESs displayed spatial differentiation. In a certain range, the degree of landscape agglomeration interacts with the degree of landscape fragmentation to promote synergistic/trade-offs effects of ESs. Our findings will provide a new analytical perspective for policymakers in the region and serve as a valuable reference for formulating targeted policies in different sub-regions.

1. Introduction

Ecosystem services (ESs) bridge natural ecosystems and human society, contributing directly and indirectly to human well-being [1]. Recently, ESs have become a prominent research hotspot in geography and ecology, attracting widespread attention within the academic community [2,3]. However, due to irreversible global climate change, intensive land use, and unsustainable resource exploitation, the decline of ESs has been confirmed at both global and regional scales [4]. To promote efficient ecosystem management, there is an urgent need to develop methods for quantifying ESs and identifying the factors that drive them [5,6]. ESs are influenced by both natural and human drivers [7]. Changes in landscape patterns are widely recognised as a significant driving force behind the variation in ecosystem service functions, as they can affect the structure and function of ecosystems by altering surface biophysical parameters, thereby impacting the provision of ESs [8,9,10,11]. Therefore, accurately understanding the driving mechanisms of landscape pattern changes in ESs is crucial for improving ecosystem service functions and making effective landscape management decisions.

ESs are not only significantly influenced by land use but also exhibit complex interrelationships [12] and potential driving mechanisms [13]. Coupled model methods, which integrate various methods, can address issues that a single model cannot resolve by leveraging the strengths of each model. Zhang et al. employed the SD-FLUS-InVEST method to predict the spatiotemporal changes and interrelationships of various future ESs [14]. Guo et al. used the LUH2-PLUS-InVEST method to project land use and carbon storage under different scenarios in China [15]. Wu et al. analysed the evolution of land cover and habitat quality in the Guangdong-Hong Kong-Macao Greater Bay Area using the PLUS-InVEST-Geographical Detectors (GDs) method [16]. Xiong et al. applied the PLUS-PLS-SEM method to assess ecosystem value in their study area [17]. A few researchers have used the InVEST-Trade-offs/Synergies method to evaluate ESs [18,19]. However, in existing studies, few researchers have investigated the driving mechanisms of ES trade-offs/synergies from a “history-present-future” perspective. In this study, the existing models are coupled, and the PLUS-InVEST-Trade-offs/Synergies-GD multi-model fusion method is used to further explore the interaction relationship between landscape pattern index and ES trade-offs/synergies.

With societal development, ESs are playing an increasingly important role in the formulation of future policies. Therefore, there is an urgent need to quantify the value of these services [20]. Studies on ESs have demonstrated that landscape characteristic indices, climate change, and various land use policies significantly impact the spatial distribution and changes in ESs [18,21,22,23,24]. However, these studies have primarily focused on the interaction between driving factors and individual elements within ESs. Driving factors not only affect the spatial distribution of ESs but also influence their trade-offs and synergy mechanisms [25]. Research by Feng et al. indicated that elevation, sand composition, silt, and vegetation type can significantly impact the trade-offs between soil moisture and soil erosion control [26]. Failing to identify the driving factors of ecosystem trade-offs and synergies may result in poor management practices and threaten the supply–demand balance of ecosystems [27]. As spatiotemporal and multi-scale analyses of ecosystems deepen, there remains a gap in exploring the interactions between landscape pattern indices and the synergistic/trade-offs relationships within ESs.

Urbanisation is an inevitable process driven by the continuous improvement of social productivity. With the large-scale expansion of urban construction land and fundamental changes in land cover types, there is an intensifying pressure on the regional ecological environment [28,29]. In recent years, the Chengdu-Chongqing Economic Circle has emerged as a new growth engine for China’s economy, which is attributed to its unique foundations and advantages [30]. In line with Sichuan Province’s 14th Five-Year Plan and the Chongqing 14th Five-Year Plan, along with related research [31,32,33], this study divides the Chengdu-Chongqing Economic Circle into four urban agglomerations: the Chengdu Plain urban agglomeration, the Southern Sichuan urban agglomeration, the Nanchong-Suining-Guangan (NC-SN-GA) urban agglomeration, and the Chongqing-Dazhou (CQ-DZ) urban agglomeration.

This study focuses on the Chengdu-Chongqing Economic Circle as the research area, using the PLUS model to identify patterns of land use change from historical to future scenarios. Simultaneously, the InVEST model was employed to evaluate four ecosystem services: carbon storage (CS), soil conservation (SC), water yield (WY), and habitat quality (HQ) within the study area, and the synergistic/trade-offs relationships among these ESs were quantified using the Trade-offs/Synergies model. Finally, the GD model was applied to analyse the intrinsic driving effects between landscape pattern indices and the synergistic/trade-ofsf relationships of ESs. Based on this, the main objectives of the study are as follows: (1) to propose a multi-model fusion method for ES analysis (PLUS-InVEST-Trade-offs/Synergies-GDs); (2) to explore the synergistic and trade-offs effects of ESs within the Chengdu-Chongqing Economic Circle and its urban agglomerations; and (3) to investigate the influence of intrinsic driving forces on the synergistic and trade-offs effects between landscape pattern indices and ESs during the urbanisation process. The findings can guide adjustments to regional land use policies and urban development, while also offering new insights for analysing regional ESs. Overall, our work may aid policymakers in developing regional landscape management schemes. This study will also serve as a valuable methodological reference for future research.

2. Materials and Methods

In this study, we propose the ES analysis method of PLUS-InVEST-Trade-offs/Synergies-GDs. Through this method, the driving force of the ES synergistic trade-offs effect in Chengdu-Chongqing economic circle is discussed. The flowchart is shown in Figure 1.

2.1. Study Area

The Chengdu-Chongqing Economic Circle is located in southwestern China (101°57′ E–108°56′ E, 27°40′ N–32°19′ N), centred around Chengdu and Chongqing, and includes 15 prefecture-level cities in Sichuan Province and 27 districts and counties in Chongqing, covering a total area of 185,000 km². It is located in the hinterland of the Sichuan Basin and the upstream of the Yangtze River, with low terrain on its western side and higher terrain on its eastern side, featuring diverse and intricate landform types. The research area is shown in Figure 2. The Chengdu-Chongqing Economic Circle has an elevation ranging from 50 to 5732 m, primarily characterised by basins and hills, surrounded by mountains to the southwest and east, and enjoys a subtropical monsoon climate with high temperatures and rainfall in summer. In 2021, the Chengdu-Chongqing Economic Circle accounted for 6.8% of the national population and 6.5% of the total national economy, serving as an important platform for China’s western development strategy and a key demonstration area for promoting new urbanisation. It is an important urbanisation area in the country, with the advantage of bearing east and enlightening west, connecting north and south. It has excellent natural endowments, a strong comprehensive carrying capacity and a relatively sound transportation system. While the Chengdu-Chongqing Economic Circle boasts the highest urbanisation rate in southwestern China and has 25 national-level nature reserves, the rapid urbanisation process has increasingly highlighted the conflict between urban land expansion and the ecological boundaries of nature reserves. Exploring the interaction between the synergistic/trade-offs relationships of ESs and landscape pattern indices from the perspective of the Chengdu-Chongqing Economic Circle and its urban agglomerations is an urgent need for maintaining environmental security in the region.

2.2. Data Sources

The data used in this study primarily include land use, socioeconomic, and topographical data, as shown in

Table 1. Data sources used in this study.

Data Type	Data Name	Data Source	Resolution
Environmental Data	Land Use	Chinese Academy of Sciences Resource and Environment Science Data Center (https://www.resdc.cn (accessed on 5 April 2024))	1 km
	Soil	Harmonized World Soil Database (HWSD)	1 km
	Precipitation	WorldClimv2.1 (https://www.worldclim.org (accessed on 5 April 2024))	1 km
	Temperature
	Evapotranspiration	National Tibetan Plateau Data Center	1 km
	DEM	Geospatial Data Cloud	30 m
Socioeconomic Data	Population	Resources and Environmental Science Data Center (https://www.resdc.cn (accessed on 17 March 2024))	1 km
	GDP
	Distance to the Road	National Geomatics Center of China (https://www.ngcc.cn (accessed on 15 March 2024))	30 m
	Distance to the Highway
	Distance to Buildings
	Distance to the Railroad
	Distance to the River

. Specifically, the land use data cover three periods (2000, 2010, and 2020) from the Chengdu-Chongqing Economic Circle land use type dataset provided by the Chinese Academy of Sciences. The land use types are classified according to the latest national classification standards, with a data resolution of 30 m. The socioeconomic data encompass various factors, including Gross Domestic Product (GDP), population, nature reserves, roads, railways, buildings, and distances. The environmental data include various factors such as the Digital Elevation Model (DEM), rainfall, and soil. In this paper, the grid data input had a resolution of 30 m × 30 m.

2.3. Methodology

2.3.1. PLUS-Based Land Use Simulation

To achieve a precise simulation of urban land use type conversion at a fine scale over specific periods, Liang et al. improved the rule extraction method based on land expansion analysis (LEAS) and the cellular automata model (CARS) using a multi-type random seed mechanism, thus constructing the PLUS model [34]. The conversion probability of the PLUS model can be expressed using the following formula:

$$P_{O,i,k}^{d=1,t}=P_{i,k}^{d=1}\times {\Omega }_{i,k}^t\times D_{i,k}^t$$

(1)

where $P_{O,i,k}^{d=1,t}$ is the conversion probability of cell $i$ at iteration time $t$ transitioning from the initial land type to the target land type $k$. $P_{i,k}^{d=1}$ is the suitability probability of land type $k$ for cell $i$ at iteration time $t$ as output by the random forest. ${\Omega }_{i,k}^t$ represents the neighborhood effect of land type $k$ for cell $i$ at iteration time $t$. $D_{i,k}^t$ is the adaptive inertia coefficient for land type $k$ for cell $i$ at iteration time $t$.

2.3.2. Multi-Scenario Settings

To analyse the future development trends of the Chengdu-Chongqing Economic Circle, this study combined the current status and relevant policies of the area to construct three major scenarios: the Natural Development Scenario (NDS), the Economic Development Scenario (EDS), and the Ecological Priority Scenario (EPS). These scenarios simulated land use and ecosystem service changes in the Chengdu-Chongqing Economic Circle by 2050. The scenarios are defined as follows:

(1)

Under the NDS, land use in 2050 follows a natural development trend, simulated using the Markov chain module in the PLUS model, based on observed land use change trends from 2010 to 2020.

(2)

The objective of the EPS is to protect areas with high ecological quality. In this scenario, the conversion of forest land, grassland, and water bodies into arable, construction, and unused lands is restricted to ensure ecological land development.

(3)

The goal of the EDS is to ensure economic development by promoting the growth of arable and construction land while limiting the development of ecological land.

2.3.3. Landscape Pattern Analysis

The calculation of landscape indices was conducted using FRAGSTATS 4.2, a widely employed method for landscape pattern analysis [35]. We studied the landscape pattern changes across 12 different research areas. These indices were calculated using a moving window approach, which illustrated the spatial distribution patterns of landscape fragments and their internal variations [36]. After multiple tests, we applied a 1000 m × 1000 m window to analyse the landscape patterns of the Chengdu-Chongqing Economic Circle. Based on existing studies [37,38], the 12 landscape pattern indices were categorised into three types: landscape fragmentation index, landscape diversity index, and landscape aggregation index, as shown in

Table 2. Landscape indices used in this study.

Type	Landscape Indices
Landscape Fragmentation Index	Landscape Division Index (DIVISION)
	Total Edge (TE)
	Edge Density (ED)
	Number of Patches (NP)
	Patch Density (PD)
	Splitting Index (SPLIT)
	Landscape Shape Index (LSI)
Landscape Diversity Index	Shannon’s Evenness Index (SHEI)
	Shannon’s Diversity Index (SHDI)
Landscape Aggregation Index	Largest Patch Index (LPI)
	Aggregation Index (AI)
	Contagion Index (CONTAG)

.

2.3.4. Ecosystem Service Assessment

The InVEST model spatially expresses the value and quantity of ESs under different land cover scenarios, enabling the quantification and visualisation of these services. It is widely used in the field of ES assessment.

(i) Carbon storage refers to the total amount of carbon stored in an ecosystem in various forms over a specific period, primarily depending on four carbon pools: aboveground biomass, belowground biomass, dead organic matter, and soil organic carbon. The calculation formula is as follows:

$$C_{total}={\displaystyle \sum }_{i=1}^n\left(\left(C_{abovei}+C_{belowi}+C_{deadi}+C_{soili}\right)\times S_i\right)$$

(2)

where $C_{total}$ represents the total carbon storage. $C_{abovei}$ is the aboveground biomass carbon storage. $C_{belowi}$ is the underground biomass carbon storage. $C_{deadi}$ is the dead organic matter carbon storage. $C_{soili}$ is the soil carbon storage. $S_i$ is the land area of land-use type $i$. According to the existing research results, this study set the carbon density of various land use types in Chengdu-Chongqing Economic Circle.

(ii) Water yield is described as the amount of water remaining after subtracting evapotranspiration from precipitation within a certain range, taking into account the soil characteristics. The calculation formula is as follows:

$$Y\left(x\right)=\left(1-{\displaystyle \frac{AE{T}_x}{P_x}}\right)\times P_x$$

(3)

where $Y\left(x\right)$ is the water yield at grid cell $x$. $AE{T}_x$ is the annual evapotranspiration at grid cell $x$. $P_x$ is the annual precipitation at grid cell $x$. In this paper, normal yearly precipitation and evapotranspiration during the study period were used (Table 1). The biophysical parameters used in this study were drawn from related papers [39,40].

(iii) Soil conservation is primarily characterised by the reduction in soil erosion and sediment retention within the region. The reduction in soil erosion is represented by the difference between potential soil erosion and actual soil erosion, while sediment retention is represented by the product of upslope sediment input and the sediment retention rate. The calculation formula is as follows:

$$Soi{l}_i=R_i\times K_i\times L{S}_i\times \left(1-C\times P\right)$$

(4)

where $Soi{l}_i$ is the soil retention at grid cell $i$. $R_i$, $K_i$, and $L{S}_i$ are the precipitation erosion factor, soil erodibility factor, and slope factor at grid cell $i$, respectively. $C$ is the cover management factor. $P$ is the conservation practice factor. The biophysical parameters utilised in this model were described in related papers [41,42].

(iv) Habitat quality refers to the potential of an ecosystem to provide the necessary conditions for the survival and reproduction of individuals and populations within a specific time and space. It is an important indicator of the degree of habitat fragmentation and biodiversity in a region. The calculation formula is as follows:

$$Q_{xj}=H_j\times \left(1-\left({\displaystyle \frac{D_{xj}^z}{D_{xj}^z+k^z}}\right)\right)$$

(5)

$$D_{xj}={\displaystyle \sum }_{r=1}^R{\displaystyle \sum }_{y=1}^{Y_r}\left({\displaystyle \frac{w_r}{{{\displaystyle \sum }}_{r=1}^R{w}_r}}\right)\times r_y\ast i_{rxy}\ast {\beta }_x\ast S_{jr}$$

(6)

where $Q_{xj}$ is the habitat quality at grid cell $x$ for land use type $j$. $H_j$ is the habitat adaptability at grid cell $x$ for land use type $j$. $D_{xj}$ is the habitat degradation at grid cell $x$ for land use type $j$. $k$ is the half-saturation constant. $z$ is the model default constant. $w_r$ is the weight of the threat source. $R$ is the number of threat sources. $Y_r$ is the number of threat source grid cells. $r_y$ is the stress value at grid cell $y$. $i_{rxy}$ is the effect of threat source $r$ on the habitat of grid cell $y$ from grid cell $x$. ${\beta }_x$ is the degree of habitat protection. $S_{jr}$ is the relative sensitivity of habitat type $j$ to threat factor $r$. This paper considered unused land, cropland, and construction land as threat sources and set the weight and maximum stress distance of each threat source. Different land use types have different sensitivity to threat sources and their own habitat suitability, which will further affect the quality of the habitat. Relevant parameters were set according to existing research [43].

2.3.5. Trade-Offs/Synergies Between ESs

Pearson correlation analysis was used to determine the relationships between ES indicators. A key step in Pearson correlation analysis is identifying the variables to statistically analyse the significance and trade-offs of the correlation coefficient ($r$). The formula is as follows:

$$r={\displaystyle \frac{\sum \left(X_i-\overline{X}\right)\left(Y_i-\overline{Y}\right)}{\sqrt{\sum {\left(X_i-\overline{X}\right)}^2}\sqrt{\sum {\left(Y_i-\overline{Y}\right)}^2}}}$$

(7)

where $X_i$ and $Y_i$ represent the values of the two variables and $\overline{X}$ and $\overline{Y}$ represent the means of $X_i$ and $Y_i$, respectively.

2.3.6. Influences on Future ESs

The GD is an effective spatial statistical method for analysing the geographical and spatial variation of variables [44]. It includes four detectors: factor detector, risk detector, ecological detector, and interaction detector. The model can measure spatial heterogeneity, detect explanatory factors, and analyse the interactions between variables. This study utilised the factor detection and interaction detection tools of the model to examine the explanatory power ($q$) of the landscape indices on the ESs of the Chengdu-Chongqing Economic Circle and revealed the interactions between different driving factors. The formula is as follows:

$$q=1-{\displaystyle \frac{1}{N{\sigma }^2}}{{\displaystyle \sum }}_{h=1}^L{N}_h{\sigma }_h^2$$

(8)

where $q$ represents the explanatory capability of the landscape indices and $h$ = 1, 2, 3, …. $L$ is the classification or stratification of the landscape indices. $N$ and ${\sigma }^2$ are the classification or stratification of the landscape indices. $N_h$ and ${\sigma }_h^2$ are the sample size and variance of layer $h$. In the geographical detector, the factor detector can detect the extent to which a factor X explains the spatial differentiation of attribute Y, and the larger the value of q, the stronger the explanatory power of the independent variable X to attribute Y. In the extreme case, a q-value of 1 indicates that factor X completely controls the spatial distribution of Y, and a q-value of 0 indicates that factor X has nothing to do with Y. The interaction detector can evaluate whether the combination of factors X1 and X2 will increase or decrease the explanatory power of the dependent variable Y, and it judges the relationship between the two factors by comparing the size relationship of q(X1), q(X2) and q(X1 ∩ X2).

3. Results

3.1. Spatial and Temporal Changes in Land Use

As shown in Figure 3, arable and forest land are the main land use types in the Chengdu-Chongqing Economic Circle, together accounting for over 80% of the total land area. From 2000 to 2020, arable land and grassland showed a decreasing trend, while construction land continued to increase. Construction land expanded by 4363.845 km², primarily converted from arable land (4528.2 km²). In contrast, forest land increased from 50,839.873 km² in 2000 to 52,322.073 km² in 2020, mainly converted from arable land and grassland; grassland decreased from 13,610.399 km² in 2000 to 10,766.339 km² in 2020, primarily transitioning to arable and forest land.

From 2020 to 2050, under the NDS, both arable land and grassland are projected to decline, with construction land increasing by 1001.301 km². The rate of construction land expansion slows but remains unchecked, continuing to encroach upon grassland and forest land. In the EDS, which prioritises food security and economic development, arable land increases by 530.688 km² and construction land rises by 2764.524 km², while forest land decreases by 645.392 km² and grassland by 2291.553 km², sacrificing ecological land to ensure economic development, which is detrimental to regional ecological protection and sustainable development. Under the EPS, policies promoting the conversion of cropland to forest lead to an increase in forest land by 1043.968 km², thereby protecting ecological land. Construction land rises by 577.094 km², with a more restrained expansion compared to the NDS and EDS, benefiting regional sustainable development.

3.2. Spatial and Temporal Changes in ESs

In this study, four ES values of different time nodes in the study area were calculated, and the calculation results are shown in Figure 4.

3.2.1. Spatial and Temporal Changes in WY

From 2000 to 2020, WY increased from 7.759 × 10¹⁰ m³ to 14.150 × 10¹⁰ m³. Compared to other ESs, WY exhibited different spatial distribution characteristics in 2000 and 2020, with a higher amount in 2020, likely linked to increased rainfall in the Chengdu-Chongqing Economic Circle that year, resulting in a higher WY around the Sichuan Basin than in 2000. Under the EPS, WY is projected to reach 8.988 × 10¹⁰ m³ by 2050. In the NDS and EDS, WY is expected to reach 8.998 × 10¹⁰ m³ and 9.122 × 10¹⁰ m³, respectively, primarily attributed to minimal construction land expansion under the EPS, which enhances runoffs capacity and increases water yield.

3.2.2. Spatial and Temporal Changes in CS

In 2020, the total CS was 18.388 × 10⁸ t, an increase of 6.601 × 10⁶ t compared to 2000, primarily due to the expansion of forested areas. High CS zones are mainly distributed in mountainous areas, while low CS zones are concentrated in cities such as Chengdu and Chongqing. The mountainous areas experience less human interference, resulting in better vegetation growth, whereas urban areas like Chengdu and Chongqing have seen rapid construction land expansion, leading to sparse vegetation. By 2050, the CS distribution characteristics are expected to be similar under all three scenarios. Under the NDS and EPS scenarios, CS is projected to reach 18.385 × 10⁸ t and 18.387 × 10⁸ t, respectively. The lowest CS is forecasted under the EDS scenario, at 18.349 × 10⁸ t, as ecological land is sacrificed to prioritise economic development, resulting in a decrease in CS.

3.2.3. Spatial and Temporal Changes in SC

The SC increased from 5.616 × 10⁹ t in 2000 to 7.851 × 10⁹ t in 2020, a change mainly driven by rainfall erosion. In 2020, rainfall in the Chengdu-Chongqing Economic Circle was significantly higher than in 2000, although factors such as vegetation cover and land management also played a role. Mountainous areas within the Chengdu-Chongqing Economic Circle have higher SC, while the plains have lower SC. By 2050, under the EDS and EPS, SC is projected to be 5.983 × 10⁹ t and 5.984 × 10⁹ t, respectively, with the lowest SC occurring under the NDS at 5.980 × 10⁹ t. Across all three scenarios, the spatial distribution of SC shows very similar characteristics.

3.2.4. Spatial and Temporal Changes in HQ

In 2020, the average HQ was 0.445, a decrease of 0.027 compared to 2000. Owing to urban expansion, the areas with declining HQ were mainly concentrated in urban clusters such as Chengdu and Chongqing. The spatial distribution of HQ is lower in urban cluster regions and higher in mountainous areas. By 2050, under the NDS and EPS scenarios, the average HQ is projected to be 0.443 and 0.444, respectively, showing only slight changes from 2020. The lowest average HQ (0.436) is expected under the EDS scenario, where urban expansion is most extensive and ecological land protection is weaker, leading to a decline in HQ.

3.3. Changes in Trade-Offs/Synergies Between ESs

Figure 5 shows that CS_SC, HQ_SC, and HQ_CS all exhibit synergistic effects, with the positive correlation of CS_SC continually strengthening, while the positive correlations of HQ_SC and HQ_CS are steadily weakening. Notably, WY_SC demonstrates a synergistic effect, with the positive correlation first weakening and then strengthening. In contrast, WY_CS and WY_HQ show trade-offs effects, where the negative correlations first increase and then decrease. Compared to the historical period, by 2050, under three different scenarios, CS_SC and HQ_SC continued to exhibit synergistic effects, but the positive correlation weakens to its lowest point, while WY_CS and WY_HQ display trade-offs, with the negative correlation reaching its maximum. Additionally, the synergy of HQ_SC remains largely unchanged from the historical period, whereas WY_SC shifts from a synergistic effect to a trade-offs effect.

Based on relevant national plans and policies, further research will be conducted at the urban agglomeration level to examine the synergies and trade-offs among various ecosystem services in different city clusters within the Chengdu-Chongqing Economic Circle across various historical periods and future scenarios.

As shown in Figure 6, both CS_SC and HQ_SC exhibited synergistic effects during different historical periods, with the strongest positive correlation observed in the Chengdu Plain urban agglomeration. HQ_CS also demonstrated synergistic effects, while WY_CS and WY_HQ displayed trade-offs. In 2000, their strongest correlation was found in the Southern Sichuan urban agglomeration, whereas in 2010 and 2020, the strongest correlations appeared in the NC-SN-GA urban agglomeration. Notably, in 2000, WY_SC exhibited synergistic effects in the CQ-DZ urban agglomeration, the Chengdu Plain urban agglomeration, and the NC-SN-GA urban agglomeration while showing trade-offs effects in the Southern Sichuan urban agglomeration. After this, WY_SC shifted to a trade-offs relationship in the NC-SN-GA urban agglomeration.

As shown in Figure 7, compared to historical periods, HQ_CS exhibits a synergistic effect under different future scenarios, with the strongest positive correlation occurring in the Southern Sichuan urban agglomeration. WY_CS, WY_HQ, and WY_SC all display trade-offs, with the strongest negative correlations also found in the Southern Sichuan urban agglomeration. Notably, the synergistic effect of WY_CS in the CQ-DZ urban agglomeration and the Chengdu Plain urban agglomeration shifts to a trade-offs relationship.

3.4. Impacts of Driving Factors on ES Trade-Offs/Synergies

3.4.1. Single Factor Detection

A study was conducted using the GD model to explore the driving relationships between the trade-offs and synergies among various ecosystem services and urban landscape indices in the Chengdu Plain urban agglomeration and CQ-DZ urban agglomeration under different scenarios in 2050. The results of the Southern Sichuan urban agglomeration and the NC-SN-GA urban agglomeration are shown in Appendix A, Appendix B and Appendix C.

As shown in Figure 8, under the EDS scenario in 2050, within a 95% confidence interval, no significant relationship between CS_HQ, CS_SC, CS_WY, HQ_WY, and any driving factors were observed in the Chengdu Plain urban agglomeration. TE has the strongest explanatory power for HQ_SC, with a q-value of 0.375. SHDI has the strongest explanatory power for SC_WY, with a q-value of 0.327. In the CQ-DZ urban agglomeration, HQ_SC, CS_WY, and CS_HQ show no significant relationship with any driving factors. SPLIT, SHEI, and LPI have the strongest explanatory power for HQ_WY, SC_WY, and CS_SC with a q-value of 0.626, 0.439, and 0.437, respectively.

As shown in Figure 9, under the EPS scenario in 2050, within a 95% confidence interval, CS_WY and SC_WY in the Chengdu Plain urban agglomeration show no significant relationship with any driving factors. LPI has the greatest impact on HQ_WY, while CONTAG is the main influencing factor for both HQ_SC and CS_SC. TE is the most significant factor affecting CS_HQ. In the CQ-DZ urban agglomeration, HQ_SC, CS_WY, CS_SC, and CS_HQ show no significant relationship with any driving factors, with LSI being the primary influencing factor for both HQ_WY and SC_WY.

Figure 10 indicates that under the NDS scenario in 2050, within a 95% confidence interval, CS_WY and SC_WY in the Chengdu Plain urban agglomeration show no significant relationship with any driving factors. SHEI has a relatively high impact on HQ_WY, with an explanatory value of 0.392. CONTAG significantly influences HQ_SC and CS_SC, with explanatory values of 0.396 and 0.420, respectively. NP and PD are the key factors affecting CS_HQ, with an explanatory value of 0.455. In the CQ-DZ urban agglomeration, HQ_SC, SC_WY, CS_WY, CS_SC, and CS_HQ show no significant relationship with any driving factors. ED and LSI are the key factors influencing HQ_WY, with an explanatory value of 0.646.

3.4.2. Interactive Factor Detection

As shown in Figure 11, under the EDS scenario in 2050, within a 95% confidence interval, HQ_WY in the Chengdu Plain urban agglomeration is the least sensitive to the interaction of various factors. For CS_HQ, the interaction between PD and LPI has the highest explanatory power, reaching 0.633. For CS_SC, the interaction between AI and CONTAG has the highest explanatory power, reaching 0.55. For CS_WY, the highest explanatory power arises from the interaction between ED and LPI, at 0.484. For SC_WY, the interaction between PD and SPLIT has the highest explanatory power, at 0.64. For HQ_SC, the interaction between LPI and TE has the highest explanatory power, at 0.581.

As shown in Figure 12, under the EDS scenario in 2050, within a 95% confidence interval, all the paired ecosystem services in the CQ-DZ urban agglomeration exhibit sensitive interacting factors. For CS_HQ, the interaction between PD and LPI has the highest explanatory power at 0.58. For CS_SC, the interaction between LPI and CONTAG has the highest explanatory power at 0.592. For CS_WY, the highest explanatory power is derived from the interaction between PD/NP and SHDI, at 0.435. For SC_WY, the interaction between AI and SHDI has the highest explanatory power at 0.47. For HQ_SC, the interaction between SPLIT and SHDI has the highest explanatory power at 0.411. For HQ_WY, the interaction between DIVISION/LPI and CONTAG has the highest explanatory power at 0.587.

As shown in Figure 13, under the EPS scenario in 2050, within a 95% confidence interval, CS_HQ in the Chengdu Plain urban agglomeration is the least sensitive to the interaction of various factors. For CS_SC, the interaction between DIVISION and CONTAG has the highest explanatory power at 0.7. For CS_WY, the interaction between ED and CONTAG has the highest explanatory power at 0.682. For SC_WY, the highest explanatory power arises from the interaction between SHEI and CONTAG, at 0.534. For HQ_SC, the interaction between DIVISION and CONTAG has the highest explanatory power at 0.784. For HQ_WY, the interaction between TE and LPI has the highest explanatory power at 0.667.

Figure 14 shows that under the EPS scenario in 2050, within a 95% confidence interval, all the paired ecosystem services in the CQ-DZ urban agglomeration exhibit sensitive interacting factors. For CS_HQ, the interaction between LSI and LPI has the highest explanatory power at 0.527. For CS_SC, the interaction between AI and CONTAG has the highest explanatory power at 0.58. For CS_WY, the highest explanatory power arises from the interaction between PD/NP and LPI, at 0.559. For SC_WY, the interaction between PD/NP and CONTAG has the highest explanatory power at 0.478. For HQ_SC, the interaction between LSI and DIVISION has the highest explanatory power at 0.446. For HQ_WY, the interaction between DIVISION and CONTAG has the highest explanatory power at 0.641.

As shown in Figure 15, under the NDS scenario in 2050, within a 95% confidence interval, all the paired ecosystem services in the Chengdu Plain urban agglomeration exhibit sensitive interacting factors. For CS_HQ, the interaction between CONTAG and SPLIT has the highest explanatory power at 0.634. For CS_SC, the interaction between LPI and CONTAG has the highest explanatory power at 0.812. For CS_WY, the highest explanatory power arises from the interaction between TE and CONTAG, at 0.681. For SC_WY, the interaction between SHDI and CONTAG has the highest explanatory power at 0.667. For HQ_SC, the interaction between TE and CONTAG has the highest explanatory power at 0.831. For HQ_WY, the interaction between AI and SPLIT has the highest explanatory power at 0.655.

As shown in Figure 16, under the NDS scenario in 2050, within a 95% confidence interval, CS_HQ and SC_WY in the CQ-DZ urban agglomeration are the least sensitive to interacting factors. For CS_SC, the interaction between LPI/SHEI and CONTAG has the highest explanatory power at 0.451. For CS_WY, the interaction between AI and SHDI has the highest explanatory power at 0.535. For HQ_SC, the interaction between AI and LPI has the highest explanatory power at 0.493. For HQ_WY, the interaction between SHEI and TE has the highest explanatory power at 0.673.

4. Discussion

4.1. Changes in the Synergistic/Trade-Offs Relationships of ESs Under Different Administrative Scales

This study used the Pearson correlation coefficient method to explore the changes in the coordination and trade-offs relationships between ESs in the Chengdu-Chongqing Economic Circle at both the macro pattern and urban agglomeration levels. From the perspective of the Chengdu-Chongqing Economic Circle, SC, CS, and HQ exhibit synergistic relationships between each pair. WY shows a trade-offs relationship with HQ and CS, while it initially has a synergistic relationship with soil retention, which later transforms into a trade-offs. This finding is consistent with previous research [45,46]. The series of ecological protection policies implemented in this region by China has, to some extent, increased the area of forest and grassland, thereby enhancing the region’s carbon sequestration and soil retention capacities [47]. Additionally, there exists a threshold for achieving the restoration of local biodiversity through artificial ecological restoration projects [48]. This explains the decreased correlation between HQ and SC, as well as between HQ and CS, and the increased correlation between CS and SC. The combined effects of ecological restoration policies, rapid urbanisation, and climate conditions have led to a decline in regional water yield [49,50,51], which provides strong evidence for the wave-like changes in the correlation between WY and SC, WY and CS, and WY and HQ.

At the urban agglomeration level, the highest correlations between CS and SC, and HQ and SC occur in the Chengdu Plain urban agglomeration, while the highest correlations between HQ and CS, WY and CS, WY and HQ, and WY and SC occur in the Southern Sichuan urban agglomeration. The trade-offs effects of ESs vary based on geographic location and over time. The “Grain for Green” (reforestation) policy has improved regional soil retention and CS capacities while reducing water yield [52]. The high level of urbanisation and limited cultivated land in the Chengdu Plain urban agglomeration and the CQ-DZ urban agglomeration show that historical WY and SC in these areas exhibit a synergistic relationship. In contrast, the Southern Sichuan urban agglomeration and the NC-SN-GA urban agglomeration exhibit a trade-offs relationship. The Chengdu Plain urban agglomeration, characterised mainly by flat terrain, has a high level of urbanisation and contains most of the national-level nature reserves within the Chengdu-Chongqing Twin City Economic Circle, which in turn enhances the synergy between CS and SC, as well as HQ and SC. In contrast, the Southern Sichuan urban agglomeration, dominated by hilly and mountainous terrain, has high vegetation coverage, which intercepts rainfall [53], resulting in significant trade-offs among WY, CS, HQ, and SC. The Southern Sichuan urban agglomeration has actively participated in the Yangtze River upper and middle reaches shelter forest construction and the Grain for Green policy, resulting in the conversion of other land types to forest and grassland, thereby improving the synergy between HQ and CS in the Southern Sichuan urban agglomeration.

The implementation of ecological restoration policies, rapid urbanisation, and climate conditions is the main factor influencing the trade-offs/synergistic of ecosystem services in the Chengdu-Chongqing economic circle. The varying responses of different urban clusters to these policies, along with their unique topography and landforms, are the reasons behind the differences in the synergistic/trade-offs relationships of ESs among these urban clusters.

4.2. Driving Mechanism of ES Synergistic/Trade-Offs Effects and Landscape Pattern Index

Changes in landscape pattern indices caused by rapid urbanisation significantly impact the ES synergies and trade-offs [54,55]. This study further explores the potential driving mechanisms between the synergies and trade-offs of ESs and landscape pattern indices in the Chengdu-Chongqing Economic Circle using the GD method. The study argues that changes in the sensitive factors affecting the synergies and trade-offs of ESs under different scenarios are due to the scenario parameter settings, which lack representativeness. Therefore, this research only discusses cases where the synergies and trade-offs of ESs and the sensitive factors remain consistent across different future scenarios. The results show that in the Chengdu Plain Urban Agglomeration, the synergy of CS_SC and HQ_SC is most sensitive to landscape aggregation, while the synergy of HQ_CS is most sensitive to landscape fragmentation. In the CQ-DZ Urban Agglomeration, the trade-offs of WY_HQ is most sensitive to landscape fragmentation.

In the Chengdu Plain Urban Agglomeration, the significant combination of driving factors with high explanatory power for the trade-offs of WY_CS and the synergy of HQ_SC is landscape fragmentation and landscape aggregation. In the CQ-DZ Urban Agglomeration, landscape aggregation is the significant combination of driving factors for the synergy of CS_SC. This aligns with existing studies [56], which indicate that land use changes are the main cause of variations in regional ESs [57] and that the changes in landscape pattern indices have a notably significant impact on ESs [58]. An increase in the fragmentation of vegetated areas negatively affects both aboveground and belowground biomass carbon density, leading to a reduction in CS [59]. Additionally, the increased fragmentation of urban construction land and natural habitat land results in a continuous decline in habitat quality [60]. The increase in the aggregation of areas with the same land type not only improves habitat quality but also enhances soil retention and carbon sequestration capacity [21,61]. The fragmentation of cultivated and forested land has been shown to increase water yield; however, it also reflects vegetation damage, indicating a decline in habitat quality [62]. Therefore, it is evident that the synergies and trade-offs of ESs in the Chengdu-Chongqing City Economic Circle result from the combined effects of landscape fragmentation and landscape aggregation and that there is significant spatial variability in the relationship between the synergies and trade-offs of ESs and landscape pattern indices.

Landscape pattern indices, to some extent, affect the trade-offs/synergistic effects of ESs in the Chengdu-Chongqing economic circle. Landscape aggregation significantly influences the trade-offs/synergistic effects of ESs in the study area, interacting with landscape fragmentation within a certain range, thereby jointly promoting the trade-offs/synergistic effects of ESs.

5. Conclusions

This study proposes a new method for regional ES analysis, integrating the PLUS model, InVEST model, trade-offs/synergy analysis, and the GD model. We aim to explore the synergies and trade-offs among various ESs in the study area during different historical periods and under future scenarios, as well as to analyse the potential driving factors behind these synergies and trade-offs. Additionally, this study is one of the few to use urban agglomerations as a starting point to examine the relationships between ESs within different urban clusters in conjunction with relevant planning and policies. Our research findings indicate the following. (1) The distribution of synergies and trade-offs among ESs exhibits spatial variability. The regions with better synergies and trade-offs among ESs are mainly concentrated in the Chengdu Plain Urban Agglomeration and Southern Sichuan Urban Agglomeration, which are jointly influenced by the topography and policy responsiveness of the different urban clusters. (2) There is spatial differentiation in the driving factors behind the synergies and trade-offs of ESs in the Chengdu-Chongqing City Economic Circle. From the perspective of interactions, landscape aggregation significantly affects the trade-off/synergistic effects of ESs in the study area, interacting with landscape fragmentation within a certain range to jointly promote these trade-offs/synergistic effects. The findings of this study not only provide a new analytical approach for policymakers involved in regional planning but also offer references and guidance for the formulation of specific policies in different subregions.