Detecting the Complex Relationships and Driving Mechanisms of Key Ecosystem Services in the Central Urban Area Chongqing Municipality, China

Abstract

Ecosystem services (ESs) are highly vulnerable to human activities. Understanding the relationships among multiple ESs and driving mechanisms are crucial for multi-objective management in complex social-ecological systems. The goals of this study are to quantitatively evaluate and identify ESs hotspots, explore the relationships among ESs and elucidate the driving mechanisms. Taking central urban area Chongqing municipality as the study area, biodiversity (BI), carbon fixation (CF), soil conservation (SC) and water conservation (WC) were evaluated based on the InVEST model and ESs hotspots were identified. The complex interactions among multiple ESs were determined by utilizing multiple methods: spearman correlation analysis, bivariate local spatial autocorrelation and K-means clustering. The linear or nonlinear relationships between ESs and drivers were discussed by generalized additive models (GAMs). The results showed that during 2000–2018, except for CF that exhibited no obvious change, all other ESs showed a decrease tendency. High ESs were clustered in mountains, while ESs in urban areas were lowest. At administrative districts scale, ESs were relatively higher in Beibei, Banan and Yubei, and drastically decreased in Jiangbei. Multiple ES hotspots demonstrated clear spatial heterogeneity, which were mainly composed of forestland and distributed in mountainous areas with high altitude and steep slope. The relationships between ES pairs were synergistic at the entire scale. However, at grid scale, the synergies were mainly concentrated in the high-high and low-low clusters, i.e., mountainous areas and urban central areas. Five ESs bundles presented the interactions among multiple ESs, which showed well correspondence with social-ecological conditions. GAMs indicated that forestland and grassland had positive impact on BI and CF. Additionally, SC was mainly determined by geomorphological factors, while WC were mainly influenced by precipitation. Furthermore, policy factors were confirmed to have a certain positive effect on ESs. This study provides credible references for ecosystem management and urban planning.

1. Introduction

In the Anthropocene era, human activities have a dominant influence on the global environmental change [1]. The accelerated processes of global industrialization and urbanization make a significant contribution to socioeconomic development [2], while inevitably disturbing ecosystem structures, functions and processes, resulting in degradation of ecosystem services (ESs) [3]. The Millennium Ecosystem Assessment (MA) [4] showed that 15 of the 24 ESs have been degraded during the past 50 years. How to optimize and improve ESs in the processes of social and economic development is an urgent issue to be addressed for urban sustainable development [5]. In recent years, many platforms and initiatives have been launched to promote the research and applications of ESs. The International Ecosystem Service Partnership was established in 2008, which provided an exchange platform for global research and practice of ESs [6]. The Aichi biodiversity targets, proposed by the 10th Conference of Parties (COP 10) to the Convention on Biological Diversity [7], emphasized enhancing the benefits to all from biodiversity and ESs. In 2012, the Intergovernmental Science-Policy Platform on Biodiversity and Ecosystem Services was established, and its core function is to conduct multi-scale comprehensive assessment of biodiversity and ESs in response to requests from decision makers [8]. The 2030 Sustainable Development Goals have recognized that multiple ESs management play an indispensable role in sustainable consumption and production, ecosystems protection and restoration, biodiversity conservation, etc. [9]. The systematic researches of ESs have recently attracted wide attention in the fields of land spatial planning and ecological protection [10,11]. Many countries had integrated the ESs concept and approach into planning and policy-making [12,13]. China is one of the first countries to implement national policies on ESs, such as the ecological redline policy, which was incorporated into China’s Environmental Protection Law [14]. However, there is no comprehensive consideration of multiple ESs in most ecosystem management strategies, which are mainly focused on achieving a single objective [15]. Therefore, it is urgent to explore the complex relationships among multiple ESs and detect the driving mechanisms in a more comprehensive way, so as to provide guidance for more rational and effective management plans.

ESs are defined as the benefits that human beings derive from the ecosystem [4]. Daily [16] proposed one of the first lists of ESs, which provided a basis for later classifications. Costanza et al. [17] grouped ESs into 17 major categories. MA [4] methodically classified ESs into provisioning, regulating, cultural and supporting services, which was the most commonly used classification. Remote sensing with spatial explicit features provides basic data support for directly or indirectly quantifying and mapping ESs through a combination of existing ecological models and data [18]. A variety of spatial data based on remote sensing are necessary inputs to Integrated Valuation of Ecosystem Services and Tradeoffs (InVEST) [19] and Artificial Intelligence for Ecosystem Services (ARIES) [20] models, which can generate maps of provisioning and regulating services displaying results in biophysical units, as well as outputs of cultural and supporting services showing results in relative rankings [21]. Combining the spatial and non-spatial responses to public attitude and preference surveys, Social Values for Ecosystem Services (SolVES) model can evaluate the perceived social values for ESs [22]. Based on this, many studies have been made towards understanding the interactions of ESs [23,24]. The relationships vary among different types of ESs, in general, a trade-off relationship is dominant between provisioning services (e.g., crop production) and regulating services (e.g., carbon fixation, soil and water conservation) [25], whereas the relationships between regulating, cultural and supporting services are synergistically dominated [26,27]. However, there are differences in the relationships between the same ES pairs at different scales in different regions. Extending the previous results to other study areas might lead to erroneous conclusions [2]. Since previous studies have primarily focused on detecting the relationships between ES pairs at regional scale [24], it is necessary to pay more attention to exploring the spatial difference of interactions at local scale and the cluster of multiple ESs across space or time, which can quantify the complex interactions among multiple ESs and provide scientific basis for establishing reasonable territorial spatial planning.

The complex relationships among ESs are deeply affected by various natural environment and social economic driving factors [27]. Determining the driving mechanisms is crucial for understanding the formation of trade-offs and synergies, and thus provides a reference for maximizing synergies and minimizing trade-offs [24,25]. Yang et al. [28] found that ESs were tightly correlated with social-ecological conditions and regulating ESs were usually aggregated with high forestland cover. Lawler et al. [29] proposed that urban expansion policies may lead to a reduction in some ESs, such as carbon fixation and crop production. At present, there is increasing effort to clarify the linear relationships between ESs and driving factors [28,30], while lacking the exploration of the nonlinear relationships. The policy effects on ESs have also not yet been sufficiently investigated [2]. Therefore, we explored the linear or nonlinear relationships between ESs and driving factors, which had rarely been discussed. Generalized additive models (GAM) provide a flexible statistical approach to identify and characterize linear or nonlinear regression effects by non-parametric smoothing functions [31].

As one of the important regions of high-quality development in western China, the central urban area Chongqing municipality (CUACM) has experienced rapid urbanization and intense human disturbances, which have induced intractable environmental challenges with grave implications for ESs supply [2]. Thus, there is a desperate need for the optimized multiple ESs management to achieve a balance between ESs supply and urban development. Biodiversity maintenance, as a supporting service, is the basis for all other ESs, and regulating services were assumed to be critical to ESs supply [4]. The CUACM is a typical mountain city, and the green spaces in mountain regions provide suitable habitats for wildlife [32]. The luxuriant forests in the mountains have great potential for carbon sequestration [33]. Biodiversity maintenance and carbon fixation have become the two mainstreams in ecological civilization construction. In 2021, the CBD COP 15 will be held in China, with the theme of “Ecological civilization: building a shared future for all life on earth” [34]. Besides, China has also set a clear climate target of achieving carbon neutrality by 2060 [35]. The undulating mountains, combined with steep slopes, leads to a widespread and serious problem of soil erosion. Although the soil erosion prevention and control in CUACM has made some achievements, there are still severe challenges. Hence, taking the CUACM as a case study, four key ESs including one supporting service and three regulating services, namely biodiversity (BI), carbon fixation (CF), soil conservation (SC) and water conservation (WC), were selected for the following analysis.

In this study, we focused particularly on proposing an integrated measurement framework for the quantitative assessment and hotspots identification of ESs, exploration of complex relationships among ESs and clarification of driving mechanisms. The following research objectives were proposed: (1) to quantitatively investigate the spatial-temporal pattern of BI, CF, SC and WC during 2000–2018; (2) to identify ESs hotspots and reveal their spatial–temporal heterogeneity in different land use types and terrain gradients; (3) to determine the trade-offs and synergies between ES pairs and detect the complex interactions among multiple ESs; and (4) to elucidate the linear or nonlinear relationships between ESs and driving factors.

2. Materials and Methods

2.1. Study Area

The CUACM (29°8′2″N–30°7′37″N, 106°14′49″E–106°58′26″E), located in southwestern China and the center of the main urban metropolitan area of Chongqing, covers approximately 6.6% of the Chongqing municipality and contains 9 districts (Figure 1). The study area is characterized by a crisscrossing distribution of mountains and rivers, and mainly composed of mountains and hills in the eastern Sichuan parallel ridge-and-valley area. Four parallel mountain ranges, namely, the Jinyun, Zhongliang, Tongluo and Mingyue Mountains, extend from northeast to southwest. The Yangtze River and Jialing River run through the CUACM from west to east. The area has a subtropical humid monsoon climate characterized by hot summers and warm and foggy winters, and the annual average temperature and precipitation total are 18.3 °C and 1199 mm, respectively. According to Chongqing urban master planning in 2007–2020, the CUACM is under a multi-center groups development mode. The urban population more than doubled between 2000 and 2018, surging from 2.97 to 7.92 million, representing 25.53% of the total population of Chongqing in 2018. With the construction of the Chengdu-Chongqing economic circle, the urbanization level increased substantially from 55.27% to 90.51%, and the gross regional domestic product (GDP) increased more than 13-fold [36]. The CUACM plays an important role in the upper Yangtze River economic belt and regional ecological protection.

2.2. Data Sources and Processing

The following datasets were used in this study, including land use types, geomorphological, meteorological, soil, vegetation, distance factors and socioeconomic dataset. The detailed descriptions of the data sources and processing were presented in

Table 1. Data sources and processing.

Datasets	Data Sources	Resolution	Data Processing
Land use types	Geospatial Data Cloud Platform [37]	30 m	Based on Landsat 7/ETM (2000), Landsat 5/TM (2010) and Landsat 8/OLI (2018) images from May to September during the vigorous growing season of vegetation with low cloud cover, the land use data were interpreted by the object-oriented classification method. The overall accuracy (OA) obtained from the confusion matrix based on the sample points from Google Earth was used in the accuracy assessment, with OA were 90.40%, 90.46% and 93.77%, respectively.
Geomorphological dataset	Geospatial Data Cloud Platform [37]	30 m	The ASTER GDEM V2 digital elevation model (DEM) data were obtained for slope, aspect, relief, terrain ruggedness index (TRI), topographic position index (TPI), terrain niche index (TNI) and sub-watershed extractions by SimDTA and ArcGIS.
Meteorological dataset	National Meteorological Information Center [38]	30 m	Temperature, precipitation and sunshine duration data were interpolated by ANUSPLIN 4.3 with data from 28 meteorological stations in the study area and its surrounding zones.
Soil dataset	Harmonized World Soil Database v1.2 [39]	1 km	The reference soil depth, salinity, sand, silt, clay, gravel and organic carbon content were resampled to 30 m resolution by cubic convolution interpolation. The soil saturated hydraulic conductivity was calculated based on the Soil Water Characteristics module in Soil Plant Atmosphere Water (SPAW).
Vegetation dataset	MODIS13 and MODIS17 [40]	250 m	The normalized difference vegetation index (NDVI) and net primary productivity (NPP) were resampled to 30 m resolution by cubic convolution interpolation.
Distance factors	Land use data	30 m	The distance to cultivated land, forestland, grassland, waters and construction land were calculated using Euclidean distance.
Socioeconomic dataset	Resource and Environmental Science Data Center of the Chinese Academy of Sciences [41]	1 km	The gross regional domestic product (GDP) was resampled to 30 m resolution by cubic convolution interpolation.
	WorldPop Dataset [42]	500 m	The population (POP) was resampled to 30 m resolution by cubic convolution interpolation.
	NPP-VIIRS-like NTL Data [43]	500 m	The nighttime light data (NTL) was resampled to 30 m resolution by cubic convolution interpolation.
	Chongqing statistical yearbook

.

2.3. Methods

The framework of the methodology established in this study can be divided into five sections (Figure 2). First, the quantitative evaluation of ESs based on the InVEST model. Second, spatial–temporal change trend analysis of ESs by the least-squares regression model. Third, ESs hotspots identification and spatial–temporal heterogeneity analysis. Fourth, the determination of the complex relationships among ESs using spearman coefficients, bivariate Local Indicators of Spatial Association (LISA) and ecosystem services bundles (ESBs). Finally, driving mechanism analysis by GAM.

2.3.1. ESs Evaluation and Validation

Habitat quality as a proxy for the potential of regions in supporting BI was modelled, and it refers to the ability of the ecosystem to provide conditions appropriate for individual survival, reproduction and population persistence [44]. Areas with high habitat quality will better support all levels of biodiversity and have relatively intact structure and function [45]. Habitat quality is a function of each threat’s relative impact and weight, habitat suitability of different land use types and sensitivity to threat factors [46]. The CF of different land use types was estimated by aggregating the amount of carbon stored in aboveground biomass, belowground biomass, soil and dead organic matter [19]. SC can be described as the difference between the potential and actual soil erosion, which is mainly determined by precipitation, soil properties, topography, vegetation coverage and anthropogenic factors, such as agricultural activities [47]. WC was evaluated by the water yield module, which is based on the Budyko curve and annual average precipitation [48]. It was modified by taking the velocity coefficient, terrain index and soil saturated hydraulic conductivity into account [49]. The detailed methodologies for ESs evaluation were described in

Table 2. Ecosystem service evaluation methodologies.

Ecosystem Services	Calculation Formulas	Parameters
Biodiversity	$\begin{array}{l}{Q}_{xj}=H_j\times \left[1-\left(\frac{D_{xj}^z}{D_{xj}^z+k^z}\right)\right]\\ D_{xj}={\displaystyle \sum _{r=1}^R{\displaystyle \sum _{y=1}^{Y_r}\left(\frac{w_r}{{\displaystyle {\sum }_{r=1}^R{w}_r}}\right)}}{r}_y{i}_{rxy}{\beta }_x{S}_{jr}\end{array}$	Qxj is the habitat quality of grid cell x in land use type j; Hj is the habitat suitability of land use type j; Dxj is the threat level of grid cell x in land use type j. The k is the half-saturation constant, which is often set as half of the maximum value of Dxj; z is the scaling parameter. R is the total number of threats; Yr indicates grid cells on threat r’s raster map; wr is threat r’s weight; ry is the relative impact of threat r in grid cell y; irxy is the impact of threat r from grid cell y on habitat in grid cell x; βx is the level of accessibility in grid cell x; and Sjr is the sensitivity of land use type j to threat r.
Carbon fixation	$C=C_{above}+C_{below}+C_{soil}+C_{dead}$	C is total carbon fixation; Cabove is carbon fixation in aboveground; Cbelow is carbon fixation in belowground; Csoil is carbon fixation in soil; and Cdead is carbon fixation of dead organic matter.
Soil conservation	$$\begin{gathered} \begin{array}{l}SEDRET=RKLS-USLE+SEDR\\ RKLS=R\times K\times LS\hspace{1em}\hspace{1em}USLE=R\times K\times LS\times C\times P\\ R={\displaystyle \sum _{i=1}^{12}\left(1.735\times {10}^{\left(1.5\lg \left(P_i^2/P\right)-0.8188\right)}\right)}\times 17.02\\ \begin{array}{ll}{K}_{EPIC}& =\left\{0.2+0.3\exp \left[-0.0256SAN\left(1-SIL/100\right)\right]\right\} \\ & \times {\left[SIL/(CLA+SIL)\right]}^{0.3}\\ & \times \left\{1-0.25C/\left[C+\exp (3.72-2.95C)\right]\right\}\\ & \times \left\{1-0.7S{N}_1/\left[S{N}_1+\exp (22.9S{N}_1-5.51)\right]\right\}\\ & S{N}_1=1-SAN/100\\ & K=(-0.01383+0.51575K_{EPIC})\times 0.1317\end{array}\\ LS=S\left[\left({\left(A_{i-in}+D^2\right)}^{m+1}-A_{i-in}^{m+1}\right)/D^{m+2}{x}_i^m{22.13}^m\right]\end{array} \\ \begin{array}{l} S=\{\begin{array}{cc}10.8\sin \theta +0.03& \theta <9\%\\ 16.8\sin \theta -0.50& \theta \ge 9\%\end{array}\\ m=\{\begin{array}{ll}0.2\hfill & \theta \le 1\%\hfill \\ 0.3\hfill & 1\%<\theta \le 3.5\%\hfill \\ 0.4\hfill & 3.5\%<\theta \le 5\%\hfill \\ 0.5\hfill & 5\%<\theta \le 9\%\hfill \\ \beta /\left(\beta +1\right)\hfill & \theta \ge 9\%\hfill \end{array}\\ \beta =\sin \theta /0.0986/\left(3\sin {\theta }^{0.8}+0.56\right)\\ C=\{\begin{array}{ll}1& f_c=0\hfill \\ 0.6508-0.3436\lg f_c\hfill & 0<f_c\le 78.3\%\hfill \\ 0\hfill & f_c>78.3\%\hfill \end{array}\\ P=0.2+0.03\times \theta \end{array} \end{gathered}$$	SEDRET is the amount of soil conservation; RKLS is the potential soil loss; USLE is the actual soil loss; and SEDR is sediment retention. R is rainfall erosivity index; Pi is monthly precipitation; and P is annual precipitation. K is soil erodibility; KEPIC is in US customary units; SAN, SIL, CLA and C are the sand, silt, clay and organic carbon contents of soil, respectively. LS is slope length and steepness factor; Ai-in is the flow accumulation; D is the grid cell size; xi is the mean of aspect weighted by proportional outflow from grid cell i; θ is percentage slope; and m is length-slope exponent. C is the cover-management factor; fc is vegetation coverage. P is the support practice factor; and θ is percentage slope.
Water conservation	$$\begin{gathered} \begin{array}{l}WC=\min \left(1,\frac{249}{V}\right)\times \min \left(1,\frac{Ksat}{300}\right)\times \min \left(1,\frac{0.9TI}{3}\right)\times Y\\ TI=\log \left(DA/\left(D_s\times \theta \right)\right)\\ Y_x=\left(1-AE{T}_x/P_x\right)\times P_x\\ AE{T}_x/P_x=1+PE{T}_x/P_x-{\left[1+{\left(PE{T}_x/P_x\right)}^{{\omega }_x}\right]}^{1/{\omega }_x}\\ PE{T}_x=K_c{\ell }_x\times E{T}_{0x}\\ {\omega }_x=Z\times AW{C}_x/P_x+1.25\\ AW{C}_x=Min\left(Rest.dept{h}_x,root.dept{h}_x\right)\times PAW{C}_x\end{array} \\ \begin{array}{l}E{T}_0={\displaystyle \sum _{i=1}^{12}E{T}_i}\\ E{T}_i=13.97d_i{D_i}^2{W}_{ti}/100\\ W_{ti}=\{\begin{array}{ll}\frac{4.95e^{0.062T}}{100}\hfill & T_i\ge 0\hfill \\ 0\hfill & T_i<0\hfill \end{array}\\ \begin{array}{ll}PAWC=\hfill & 54.509-0.132SAN-0.003SA{N}^2\hfill \\ & -0.055SIL-0.006SI{L}^2-0.738CLA\hfill \\ & +0.007CL{A}^2-2.688C+0.501C^2\hfill \end{array}\\ Z=\frac{(\omega -1.25)\times P}{AWC}\end{array} \end{gathered}$$	WC is the amount of water conservation; V is velocity coefficient; and Ksat is soil saturated hydraulic conductivity. TI is terrain index; DA is the number of grids in catchment area; Ds is the depth of soil; and θ is percentage slope. Yx is the water yield for grid cell x; AETx is the annual actual evapotranspiration; Px is the annual precipitation; PETx is the potential evapotranspiration; ωx is a nonphysical parameter of natural climate-soil properties; Kcℓx is plant evapotranspiration coefficient for each land use type; Z is the seasonal parameter; AWCx is the volumetric plant available water content; Rest.depthx is the root restricting layer depth; root.depthx is the vegetation rooting depth; and PAWCx is plant available water content. ET0 is average annual reference evapotranspiration; ETi, di, Di, Wti and Ti reference evapotranspiration, the number of days, sunshine duration, saturated water vapor density and temperature of month i, respectively.

.

The vegetation coverage has a significant positive effect in conservation practice to enhance species richness, which is commonly assessed to derive conclusions on a region’s biodiversity [50]. Net primary production (NPP) is defined as the remaining part of the organic substance amount produced by vegetation during photosynthesis and respiration, which is a direct reflection of the CF capability [51]. Therefore, BI and vegetation coverage, CF and NPP were fitted with linear regression to verify the accuracy of the models. The squares of the correlation coefficients (R²) were 0.85 and 0.82, respectively, suggested better goodness of fitting and higher accuracy. The reliability of the SC simulation results was verified by referencing related research of Liu et al. [52]. The simulation results coincided well with the reference estimates, with the correlation coefficient was up to 0.91, which indicated the credibility of the simulation results. The WC was calibrated with observed data in Chongqing water resources bulletin [53]. There was a strong correlation between observed and simulated values, with the average relative error less than 10% and R² value greater than 0.95. The model validation results showed that the simulated values were reasonable and reliable.

2.3.2. Spatial-Temporal Change Trend Analysis

At grid scale, the least-squares regression model was used to detect the regional differences in ESs change trends [54]. Based on the precision of data, the grid size was selected as 30 m × 30 m. Furthermore, we made the statistics of ESs and their changes at administrative district scale. The least-squares regression model can be expressed as follows:

$$Slope=\frac{n{\displaystyle \sum _{i=1}^n{x}_i{a}_i-{\displaystyle \sum _{i=1}^n{x}_i{\displaystyle \sum _{i=1}^n{a}_i}}}}{n{\displaystyle \sum _{i=1}^n{x_i}^2-{\left({\displaystyle \sum _{i=1}^n{x}_i}\right)}^2}}$$

(1)

where slope is the change rate of ES; n is the number of years; x_i represents the year i; and a_i is the ES in year x_i. When slope > 0, ES displays an increasing trend and vice versa. The statistical significance was tested by F-test [24].

2.3.3. ESs Hotspots Identification

ESs hotspots represent statistically significant high-value spatial aggregations [55], and they were identified by the Z-scores and P value test of the Getis–Ord value G_i* based on a distance weight matrix in spatial correlation analysis. G_i* is defined as:

$$G_i^{*}=\frac{{\displaystyle \sum _{j=1}^n{W}_{ij}{x}_j}}{{\displaystyle \sum _{j=1}^n{x}_j}}$$

(2)

where x_j is the property value of factor j; W_ij is the spatial weight between factors i and j; and n is the number of factors.

For statistically significant positive Z-scores, the larger the Z-score is, the higher the degree of high-value spatial aggregation. The areas with a confidence level greater than 90% were considered as ESs hotspots [56]. The same ecosystem has the potential to provide multiple ESs. The distribution of multiple ES hotspots could be determined by overlaying the four ESs hotspots [49,57]. The region with four ESs hotspots was defined as a class 4 of ESs hotspots, similarly, class 3, class 2, and class 1 ESs hotspots were identified.

To comprehensively reflect the detailed characteristics of topographic spatial heterogeneity of multiple ES hotspots, terrain niche index (TNI) was calculated by combining slope and elevation [58], and it was classified into ten gradients using the natural breaks classification method. Furthermore, the distribution index was used to reveal the spatial distribution of multiple ES hotspots on TNI gradients. The formulas are as follows:

$$T=\lg \left[\left(\frac{E}{\overline{E}}+1\right)\times \left(\frac{S}{\overline{S}}+1\right)\right]$$

(3)

where T is the terrain niche index; E and S are the elevation and slope of the grid cell, respectively; and $\overline{E}$ and $\overline{S}$ are the average elevation and slope of the study area, respectively.

$$P=\left(\frac{S_{ie}}{S_i}\right)/\left(\frac{S_e}{S}\right)$$

(4)

where P is the distribution index; S_ie is the area of ESs hotspots on the TNI gradient e; S_i is the area of ESs hotspots; S_e is the area of TNI gradient e; and S is the total area of the entire study area. Generally, greater P values indicate a high distribution frequency of ESs hotspots on a certain TNI gradient. When P > 1, the TNI gradient represents the dominant distribution of ESs hotspots.

2.3.4. Investigating the Complex Relationships among ESs

In order to eliminate the influences of different units and scales, ESs were standardized using a minimum–maximum normalization to obtain comparable and dimensionless data ranging from 0 to 1 [59]. Considering the sensitivity of this standardization to extreme values, the ESs were first transformed by winsorization where ESs with values less than the 5% quantile and greater than the 95% quantile were assigned values of the 5% and 95% quantile, respectively [60,61].

Spearman rank correlation coefficients were used to illustrate the trade-offs and synergies between ES pairs at a regional scale since the ESs data did not conform to normal distribution [62]. The formula is as follows:

$$r=1-\frac{6{\displaystyle \sum _{i=1}^n{d}_i^2}}{n(n^2-1)}$$

(5)

where r is the correlation coefficient; n is the number of ES pairs; d_i is the rank difference. When r > 0, there is a positive correlation, which indicates that the paired ESs were synergistic; when r < 0, there is a negative correlation suggesting trade-offs between ES pairs. In accordance with Raudsepp-Hearne et al. [63], the strength of the correlations was described by using the following classification: strong correlation (|r| ≥ 0.5), moderate correlation (0.5 > |r| ≥ 0.3), weak correlation (0.3 > |r| ≥ 0.1) and no correlation (|r| < 0.1). T-test was performed to test for statistical significance [64]. * p < 0.05 was considered statistically significant.

Bivariate local Moran’s I was used to explore the trade-offs (spatial dispersion) and synergies (spatial aggregation) between ES pairs on spatial scale [65]. The formula is as follows:

$$I_{ij}=\frac{x_i^a-{\overline{x}}_a}{{\sigma }_a}{\displaystyle \sum _{j=1}^n\left(W_{ij}\frac{x_j^b-{\overline{x}}_b}{{\sigma }_b}\right)}$$

(6)

where I_ij is bivariate local spatial autocorrelation coefficient; n is the number of grid cells; W_ij is spatial weight matrix; x_i^a and x_j^b are the value of ES a and ES b in grid cell i and j, respectively; and ${\overline{x}}_a$, ${\overline{x}}_b$, σ_a and σ_b are the average and variance of ES a and ES b, respectively. The statistical significance of bivariate local Moran’s I was examined by permutation tests; and 999 permutations were used in this study [66]. The significance value for spatial correlation between ES pairs were classified as significant (* p < 0.05), very significant (** p < 0.01), and extremely significant (*** p < 0.001).

Positive spatial autocorrelation was when similar values clustered together on a map, and negative spatial autocorrelation was when dissimilar values clustered together. Therefore, the four cluster types of local spatial autocorrelation generated by bivariate LISA were reclassified as trade-offs and synergies. The high-high and low-low clusters were reclassified as synergies, the high-low and low-high clusters were reclassified as trade-offs.

ESBs can reflect the aggregation and combination of multiple ESs on spatial and temporal scales, which can quantify the complex interactions among multiple ESs. ESBs were identified by K-means clustering algorithm [33]. To minimize the total error sum of squares, the number of initialized repeating runs and iterations were set to 150 and 1000, respectively [11,63]. The optimal number of clusters was determined by the ratio of the between-cluster sum of squares and the total sum of squares. The higher the ratio, the greater the difference between clusters and the smaller the difference within clusters. After series of trials and analysis, when the number of clusters was set to five, the ratios were the highest, which were 80.16%, 81.12% and 80.07% in 2000, 2010 and 2018, respectively. The five ESB types were visualized using petal plots, and the length of each petal represents the average of ESs [23].

2.3.5. ESs Driving Mechanisms

GAM was employed to examine the linear or nonlinear relationships between driving factors and ESs in 5468 grid cells of 1 km × 1 km located across study area using the mgcv package of R [61]. GAM had advantage in developing realistic response curves because it fit non-parametric smoothers to the data without requiring priori specifications to describe nonlinearity [67]. All driving factors entered the GAM as smooth terms. The general expression of the model is as follows:

$$g(\mu )={\beta }_0+x_i\beta +{\displaystyle \sum _{j=1}^n{f}_j(x_{ij})}+\epsilon$$

(7)

where g is the link function; μ is the expectation of dependent variable Y; β₀ is the intercept term; x_i and β are the independent variables and parameters of fixed effects; f_j(x_ij) is non-parametric smoothing functions; and $\epsilon$ is the error term.

Based on the principles of scientificity, representativeness and accessibility, 32 natural environment and socioeconomic factors were selected to build the driving factors database, including: (1) seven geomorphological factors, (2) land use types, (3) NDVI and NPP, (4) four meteorological factors, namely, temperature, precipitation, rainfall erosivity index and reference evapotranspiration, (5) four soil factors, namely, soil types, erodibility, saturated hydraulic conductivity and organic carbon content, (6) five distance factors, namely, the distance to cultivated land, forestland, grassland, waters and construction land, (7) three socioeconomic factors, namely, POP, GDP and NTL, and (8) six policies, namely, CUACM overall urban planning in 2004–2020 (CMOUP), multi-center groups development strategies in Chongqing urban master planning in 2007–2020 (MGDS), regulations on developments and controls of the four parallel mountains (RDCFM), Chongqing ecological function regionalization (CEFR), Chongqing ecological redline (CER), and beautiful landscape city planning in CUACM (BLCP).

The collinearity was detected by the variance inflation factor (VIF) and driving factors with a VIF greater than 5 were sequentially removed [68]. Akaike’s information criterion (AIC) and adjusted determination coefficient (R²_adj) were applied to check the goodness-of-fit of the model [69]. The AIC and R²_adj of four GAMs between ESs and driving factors were −25,399.34, −21,417.56, −11,355.13, −25,203.94 and 0.869, 0.859, 0.603, 0.666, respectively, suggested good model performance and explanatory power. T-test and F-test were used for categorical and numerical variables, respectively. The statistical significances were classified as significant (* p < 0.05), very significant (** p < 0.01), and extremely significant (*** p < 0.001).

3. Results

3.1. Spatial Patterns of ESs

All the four ESs demonstrated clear spatial clustering on the study area and exhibited different change trends. High BI areas with high habitat quality were concentrated in Jinyun, Zhongliang, Tongluo, Mingyue, Huaying and southeastern mountainous areas (Figure 3(a1–a3)). The Yangtze River, Jialing River and other water bodies also had higher BI. The BI displayed a decreasing trend, with the habitat quality index declining from 0.5 to 0.45 during 2000–2018. Areas with decrease trends were mainly distributed in peri-urban expanded regions, and most of which (99.51%) was nonsignificant decrease areas (Figure 3(a4)). Similar to BI, CF also showed a high spatial pattern in mountainous areas (Figure 3(b1–b3)). The change of CF was not obvious overall, with a trend of decrease first and then increase. In most of the study area (89.42%), CF remained unchanged. Additionally, a similar nonsignificant decreasing trend was observed in peri-urban expanded regions (Figure 3(b4)). The high values of SC were not evenly distributed in mountainous areas (Figure 3(c1–c3)), showing obvious topographic heterogeneity. There was a trend of decrease in SC from 1289.56 t/hm² to 921.90 t/hm². Approximately 81.36% of the study area showed a decreasing trend in SC, of which 97.82% decreased was nonsignificant (Figure 3(c4)). The high values of WC were mainly located in mountainous areas, which was profoundly influenced by precipitation (Figure 3(d1–d3)). The average WC depth decreased from 192.43 mm in 2000 to 141.79 mm in 2018. The areas with decreasing WC accounted for 65.20%, and were mainly found in mountainous areas and the northern half of the study area, whereas the increases occurred in the southern portion (Figure 3(d4)). Additionally, most of the changes were nonsignificant. The change trend of WC from 2000 to 2018 had good consistency with precipitation change. Due to the frequent interference of human activities, all ESs in urban areas were extremely low.

At the administrative districts scale, Beibei had the highest level of BI, followed by Banan and Yubei (Figure 4a). The highest average values of CF were also detected in Beibei, followed by Shapingba, Jiulongpo and Yubei (Figure 4b). The highest average SC and WC occurred in Banan, followed by Beibei and Yubei (Figure 4c,d). All ESs in Yuzhong were the lowest and far lower than those in other districts. It is worth mentioning that there were some increases in CF in Nanan, Yubei and Yuzhong. In addition to this, the general trends of ESs for all districts were descending. The maximum reductions of BI and CF were observed in Jiangbei. The dramatic decreases in SC and WC were observed in Beibei, Yubei and Jiangbei.

3.2. Spatial Heterogeneity of ES Hotspots

The results of hotspots analysis suggested clear characteristics in temporal and spatial distribution for multiple ES hotspots (Figure 5). During 2000–2018, the total area of multiple ES hotspots exhibited a slightly increasing tendency, with the area proportion increased from 31.55% (1726.45 km²) to 33.19% (1816.22 km²). Specifically, increases in class 1 ES hotspots were strongest. Class 2 ES hotspots also behaved a general increasing trend of increased first and then decreased. In 2018, the proportion of both to all ES hotspots reached as high as 53.46% and 32.16%, respectively. Class 3 and class 4 ES hotspots had a tendency to decrease, accounting for only 12.80% and 1.58% of all ES hotspots in 2018, respectively.

There was considerable variability in ESs delivered by different land use types, which led to the spatial heterogeneity of multiple ES hotspots (Figure 5). Class 1 ES hotspots occurred mainly in cultivated land, forestland and waters, with the average proportion of 55.41%, 21.78% and 15.23%, respectively. Class 2 ES hotspots were mainly composed of forestland and cultivated land, covering 83.68% and 11.76%, respectively. The average contribution of forestland to class 3 ES hotspots reached 91.12%. Although the area of forestland showed a decreasing trend, the area proportion showed an increasing trend, which was caused by the decrease of class 3 ES hotspots. The most predominant land use type that composed class 4 ES hotspots was forestland, accounting for 94.10% in 2000, 91.90% in 2010 and 94.07% in 2018. The non-hotspots were primarily covered by cultivated land and construction land. Generally, forestland was a major contributor to ESs and also had a significantly higher area proportion of multiple ES hotspots compared to other land use types. Among them, the average area proportion of class 1, 2, 3 and 4 ESs hotspots were 17.09%, 45.42%, 27.48% and 4.69%, respectively. The area of grassland was small, but it had an important role in ESs supply, second only to forestland. During 2000–2018, multiple ES hotspots accounted for an average of 84.02% of grassland. The average area proportion of multiple ES hotspots of waters ranged from 80.24% to 58.70% in a decreased trend. Additionally, average 70.27% of waters was class 1 ES hotspots. Most of unused land, cultivated land and construction land were non-hotspots, which accounted for averages of 70.25%, 84.70% and 92.86% of such kinds of land use, respectively.

The distributions of multiple ES hotspots had distinct spatial heterogeneity in different TNI gradients (Figure 6). In the 1st TNI gradient and high TNI gradients (6th–10th), the average distribution index of class 1 ES hotspots were greater than 1, indicating that it was primarily located in two parts of the study area: water areas with low altitude and slope, such as the Yangtze River and Jialing River, and the transition zone between mountainous areas and flat areas with medium-high altitude and slope. The distribution index of class 1 ES hotspots tended to decrease in the 1st TNI gradient over time, while it increased in the 6th to 10th TNI gradients. Class 2 and class 3 ES hotspots were mainly distributed in high TNI gradients (6th–10th), with average distribution index of 2.18 and 2.52, respectively. In the 10th TNI gradient, an upward trend in distribution index was observed of class 2 ES hotspots, while the trend of class 3 ES hotspots was downward during 2000–2018. The distribution index of class 4 ES hotspots presented a tendency to increase with an increase in TNI gradients, especially after the 7th TNI gradient. Although the distribution index showed a decreasing trend in the 10th TNI gradient during 2000–2018, it was still much larger than that in other gradients, which indicated that the mountainous areas with high altitude and steep slope were always the predominant distribution areas of four ES hotspots.

3.3. Spatial-Temporal Trade-Offs and Synergies between ES Pairs

During 2000–2018, the significant positive correlations (Spearman coefficient: r > 0; p < 0.001) were found between all ES pairs in the entire study area (Figure 7). It demonstrated that all six pairs of ESs exhibited synergistic relationships. In 2000, there was a strongest positive correlation between BI and WC. The relationship between SC and WC was also strong synergistic. In 2010, the highest positive correlation was observed between BI and CF. Additionally, WC also had strong positive correlation with BI and SC. In 2018, the synergistic relationship between BI and CF remained strongest. Strong positive correlations were also found between SC with BI and CF. In addition, the correlation coefficients among BI, CF and SC showed a tendency to increase during 2000–2018, which indicated that the degree of synergies were gradually enhanced. With the relationship between BI and WC changed from strong correlation to weak correlation, the synergy was diminished. The relationships between WC with CF and SC showed a general weakening trend of slightly enhanced first and then weakened.

There were obvious differences in spatial distributions of trade-offs and synergies between ES pairs (Figure 8). For BI and CF, synergies were always the dominant relationship in the study area, with the area proportion ranging from 68.53% to 68.18%. Among them, significant synergy accounted for average 59.20%, which was mainly distributed in the flat areas between mountains. The extremely significant synergy showed a gradually increasing trend, with an average area proportion of 33.04%, which were mainly concentrated in the high-high and low-low clusters, i.e., mountainous areas and urban central areas. The extremely significant trade-offs accounted for 4.13% of the study area and were spatially aggregated in waters, especially in Yangtze River and Jialing River. For BI and SC, synergies and trade-offs accounted for average 35.13% and 8.94%, respectively. The synergies were found mainly in the intermountain flat areas and scattered in Jinyun, Zhongliang, Tongluo and Mingyue Mountain. The trade-offs were mainly distributed in Yangtze River and Jialing River, and also observed in Qiaoping mountain and southeastern mountainous areas. For BI and WC, the area of synergies was comparable to trade-offs, accounting for 27.33% and 20%, respectively. The extremely significant synergies were distributed in two dominant parts: one located in mountainous areas and the other in urban areas. During 2000–2018, the former exhibited a downward trend, while the latter showed the opposite trend. The trade-offs tended to increase, and the most dramatic increase was observed in significant trade-offs, which located in some mountainous areas and southeastern areas. The extremely significant trade-offs still mostly occurred in Yangtze River and Jialing River. Apart from Yangtze River and Jialing River changing from trade-offs to synergies, the trade-offs and synergies between CF with SC and WC followed a spatial pattern similar to the relationships between BI with SC and WC. In addition, similarity was also observed in general change trends. For SC and WC, the synergies were dominated by the extremely significant synergies, which mainly occurred in urban areas, mountainous areas with an elevation of 200~700 m and a slope of 5~15°, as well as waters such as Yangtze River and Jialing River. Both the area of synergies and trade-offs showed the same increasing trend, but the magnitude of changes was most pronounced in significant trade-offs. The significant trade-offs increased by 441.40 km², which were mainly located in some mountainous areas and southeastern intermountain hilly and flat areas.

3.4. ES Bundles among Multiple ESs

Five ESBs were determined by K-means cluster analysis (Figure 9). In 2000, ESB1 was spatially clustered in Jinyun, Zhongliang, Tongluo, Mingyue, Huaying and southeastern mountainous areas with an average elevation of 535.83 m and slope of 9.57°. It was predominantly covered by the highest proportion of forestland (92.93%), occupying 6.30% of the study area. ESB1 was especially characterized by the highest BI, CF, SC and the second-highest WC, which was only 0.07 lower than that of ESB2. ESB2 were interlaced with ESB1, mainly distributed in mountainous areas with an average elevation of 515.78 m and slope of 6.45°, and mostly covered by forestland (87.84%). It showed the highest WC, the second-highest BI and CF, and moderate SC. ESB3 was observed in the piedmont transition zone between mountainous and flat areas with an average elevation of 436.57 m and slope of 7.29°, contained 94.51% of cultivated land. SC was dominant in ESB3 and second only to that of ESB1. The mean values of BI and CF were below the regional averaged values, whereas the mean value of WC was slightly above the regional average value. Covering 63.45% of the study area, ESB4 were widely spread in the intermountain flat areas, the trough and valley at the top of anticline mountains, as well as water areas. The land use pattern of this bundle was characterized by larger cultivated land ratios (92.24%), followed by waters (4.43%). ESB5 was characterized by dense urbanization and concentrated in urban areas, with the proportion of construction land of 90.58%. All ESs in this bundle were at the lowest level.

In 2010, ESB1 was increased and clustered in mountainous areas with an average elevation of 520.04 m and slope of 7.46°. It supplied the highest BI, CF, WC and the second-highest SC. Different from ESB2 in 2000, ESB2 were mainly distributed in waters, such as the Yangtze River and Jialing River, with an average elevation of 281.64 m and slope of 3.26°. The proportion of waters in this bundle was up to 66.66%, followed by forestland (4.43%). ESB2 provided a relatively higher BI and lower CF, SC and WC. The characteristics and distributions of ESB3, ESB4 and ESB5 were similar to those in 2000. The slight difference was that ESB3 showed the highest SC compared to other ESBs.

In 2018, ESB1 contributed strongly to WC, provided a relatively higher BI and CF, and moderate SC. It was located in mountainous areas with an average elevation of 489.40 m and slope of 6.59°, mainly consisting of forestland (79.34%) and grassland (10.77%). ESB2 was spatially clustered in mountainous areas with an average elevation of 538.99 m and slope of 7.97°, covered mainly by forestland (92.36%). It was characterized by the highest BI and CF, the second-highest SC, which was second only to that of ESB3. However, the mean value of WC was below the regional averaged value. ESB3, ESB4 and ESB5 had quite similar distributions and combinations of ESs with those in 2000 and 2010. Notably, as urban sprawl, the area of ESB5 presented a dramatic increasing tendency from 444.08 km² to 947.33 km². This increase mainly occurred in new cities groups and the periphery of the original cities.

4. Discussion

4.1. Exploring the Driving Mechanisms of ESs

ESs were deeply influenced by natural environment, socioeconomic and policies factors [33]. Different from previous studies [70,71], we quantitatively investigated linear or nonlinear relationships between ESs and driving factors based on GAM. Although relevant studies had elaborated that land use was the primary driver of ES changes, most of the results were only presented in the form of qualitative description [32,71]. There was particularly a lack of focus on political drivers of ESs [30]. We included land use types and policies as fixed effects in GAM, and the parameter estimates illustrated the influence of different independent variables on ESs (

Table 3. Results of GAM characterizing the relationships between ESs and fixed effects.

Fixed Effects	Biodiversity		Carbon Fixation		Soil Conservation		Water Conservation
	Estimate	Pr(>|t|)	Estimate	Pr(>|t|)	Estimate	Pr(>|t|)	Estimate	Pr(>|t|)
Cultivated land	0.424	<0.001 ***	0.091	0.01 **	−0.218	<0.001 ***	0.335	<0.001 ***
Forestland	0.939	<0.001 ***	0.741	<0.001 ***	0.172	<0.001 ***	0.496	<0.001 ***
Grassland	0.896	<0.001 ***	0.264	<0.001 ***	−0.175	<0.001 ***	0.638	<0.001 ***
Waters	0.878	<0.001 ***	—	—	−0.269	<0.001 ***	—	—
Construction land	0.056	<0.001 ***	—	—	−0.221	<0.001 ***	—	—
MGDS	—	—	—	—	—	—	—	—
RDCFM	—	—	0.011	0.01 **	0.017	0.004 **	—	—
CEFR	—	—	—	—	—	—	—	—
CER	—	—	0.021	<0.001 ***	—	—	−0.024	<0.001 ***
BLCP	0.013	<0.001 ***	0.011	0.003 **	—	—	−0.023	<0.001 ***
Intercept	0.970	<0.001 ***	1.028	<0.001 ***	1.349	<0.001 ***	1.001	<0.001 ***

The statistical significance levels were very significant (** p < 0.01) and extremely significant (*** p < 0.001). —: There were not statistically significant. MGDS: multi-center groups development strategies; RDCFM: regulations on developments and controls of the four parallel mountains; CEFR: Chongqing ecological function regionalization; CER: Chongqing ecological redline; BLCP: beautiful landscape city planning in CUACM.

). Regarding the land use types, BI had extremely significant strong positive correlations with forestland, grassland and waters, the parameter estimates were up to 0.939, 0.896 and 0.878. CF were most enriched in forestland, followed by grassland. Strong and moderate positive correlation was found between WC with grassland and forestland, respectively. These parallels the majority of previous findings that, forestland and grassland can not only support high biodiversity by providing resources and habitats for species, but also provide high carbon sequestration and had low runoff potentials that contributed greatly to water conservation [61,72]. However, SC was shown to be weakly correlated with land use types. This may be because SC was more strongly affected by geomorphological factors, which was confirmed by Baró et al. [33].

The effects of policies on ESs, although significant, appear to be modest. The RDCFM, CER and BLCP had weak positive effects on BI, CF and SC, while they had weak negative effects on WC (Table 3). Two possible reasons may be able to account for these. First, lack of scientific and sufficient information on ESs to serve as comprehensive evaluation criteria for delineation of protected areas resulted in the co-existence of over-protected and conservation gaps in policies [73]. Second, the policies were formulated based on the conditions at a certain time, while ecosystems were dynamic, and follow changes in the natural and socioeconomic environment [74]. Despite effectiveness limitations and time lag of policies, our results indicated that policies had positive effects on the BI, CF and SC. However, WC was mainly determined by climate factors and fluctuated with precipitation variation. The study of Cui et al. [70] had demonstrated that the spatial change of WC was consistent with that of precipitation.

The results of GAM for ESs revealed the relative effect and importance of various driving factors, and the smooth functions showed the linear or nonlinear relationships between ESs and driving factors. The distance to forestland was the main influencing factor for BI with the highest F-value of 17.05, followed by the distance to cultivated land (F = 12.51) and elevation (F = 9.174) (Table S1, Figure 10). There was a decreasing trend generally in BI with the increasing distance to forestland. The observed appearance of the crest was closely related to waters with a higher BI, which were crisscrossed and interlaced with mountains. The CUACM was characterized by northeast–southwest trending parallel mountain ranges and west–east trending rivers, interspersing with flat valleys and urban groups. The complex spatial mosaic of landscapes resulted in a strongly nonlinear relationship between BI and the distance to cultivated land, waters and construction land. BI showed a general increasing trend of decrease first and then increased with increasing elevation and TNI. This was because waters were primarily located in the low TNI with low altitude and slope. Relevant studies had shown that waters created diverse and stable habitats for waterfowl and fish to forage, refuge and breed [75]. In addition, the same relationships were also observed between BI with NPP and soil organic content. NDVI had a very significant positive effect on BI, and an estimated degree of freedom equal to one confirmed the linear relation. There was a general negative relationship between BI and precipitation. One likely explanation for this may be that a large amount of precipitation had the potential to induce geological disasters in mountainous areas, which had a certain negative impact on BI. A similar negative correlation had been discovered in a study of the Yangtze River Economic Belt [76]. BI decreased nonlinearly with increasing population density. The ever-increasing population led to an expanding demand for natural resources and produces [77], which caused dramatic interference with BI [78].

The strongest effects on CF were soil organic content (F = 113.68), distance to forestland and grassland (Table S1, Figure S1). CF showed an extremely significant positive correlation with soil organic content, which had been suggested to be one of the main determinants of plant species distributions and ecosystem functioning [79]. As the distance to forestland increased, CF was overall decreased. As it can be observed, crest clearly appeared because grassland and cultivated land had a certain potential for carbon sequestration [80]. Consistent with previous studies [81], forestland was the major component of carbon pool in the study area, with a contribution of 76.92%. Thus, CF showed a trend of decreased first and then increased with increasing distance to grassland. The higher elevation, TRI and TNI, the higher forestland coverage and the higher CF. Affected by the trend of the mountains from northeast to southwest, the area proportion of forestland was relatively higher on east and west slope. The CF was fluctuated with aspect variation and peaked approximately at 100° and 275°. The smooth fitting curve between CF and TPI appeared a V-shape. The areas, with the TPI equal to zero, were mostly flat areas, such as rivers and lakes, which had the lowest CF. The higher or lower the TPI, the closer it was to the ridge or valley and the higher the CF. The CF had a very significant linear increase trend with increasing NDVI, and it also had a general tendency to increase with NPP. The result was in line with recent findings in which vegetation coverage had a notable positive impact on NPP, and higher NDVI significantly improved CF [24,57]. Similar to that of BI, CF was negatively correlated with precipitation. Due to the complex spatial pattern of landscapes, CF had nonlinear relationships with the distance factors. CF exhibited a diminishing trend with the increased in population and GDP. Generally, in areas with a higher population and GDP, urbanization was usually higher. The demand for urban land had gradually increased, and a large amount of ecological land was occupied, which may be the major causes for CF loss [2].

Geomorphological factors were the dominant factors for SC, such as TPI (F = 1313.929), TNI and slope (Table S1, Figure S2). There was a negative and saturating relationship between SC and TPI. Specifically, SC decreased rapidly with increasing TPI when TPI was less than zero, whereas there was little change in SC when TPI was greater than zero. It means that the closer the valley, the higher the SC. According to Tsaaer et al. [82], the valleys were considered to be favorable sites for soil deposition, with large forest area and relatively lower soil loss. The areas with higher TNI and steep slope usually had higher vegetation coverage, the actual soil loss was considerably lower than potential soil loss, resulting in a higher SC [83]. SC displayed a tendency of increased first and then decreased with elevation. When the elevation was greater than 700 m, the higher the elevation, the higher the vegetation coverage, which can effectively intercept precipitation and reduce soil erosion, resulting in a decrease in sediment deposition from upstream, and a downward trend in SC [52]. Influenced by the mountain range trend, SC exhibited a double hump-shaped relation with aspect, peaking at about 90° and 270°. The CF had a positive linear correlation to NDVI and precipitation, and tended to nonlinearly increase with increasing NPP and soil organic content. The higher the vegetation coverage, the higher the NDVI and NPP, the thicker the litter layer, the more developed root systems, the higher the soil organic content, which can prevent soil erosion from rain splash [24], resulting in a decrease in actual soil loss and an upward trend in SC. The precipitation had the most direct positive effect on the potential soil loss, resulting in a higher SC. In general, the closer to the forestland and the further away from the construction land, the higher the SC.

Our findings suggest that WC was mainly controlled by precipitation (F = 472.458), which was in accordance with previous studies [70] (Table S1, Figure S3). Generally, in areas with a higher precipitation, water yield was usually higher, which led to higher WC. WC was also influenced by other factors, once the resistance of the controlling ability of precipitation for WC was broken, WC decreased significantly. Singh et al. [84] also corroborated a high rainfall did not guarantee a high WC in a region as it depends on variety of factors. WC increased linearly with elevation and nonlinearly with TNI. When the slope was less than 30°, WC increased rapidly with increasing slope, whereas there was a slightly decreased in WC when slope was greater than 30°. The higher the TNI with higher altitude and slope, the higher the vegetation coverage. The dense canopy had the ability to trap precipitation [24], the thick litter layer and developed root systems made more water infiltrated into the soil [26], and thus there will be an increase in WC. However, as soon as the slope exceeds 30°, water holding capacity strongly decreased, resulting in reduction of WC. Thus, there was lower WC on the ridges with steep slopes, while valleys had the highest WC. NDVI, NPP and soil organic content had extremely significant positive effects on WC, although slightly decreased or saturated at the highest values. Plants can not only absorb and store some precipitation, but also increase regional precipitation through plant transpiration [49]. The higher the soil organic content, the higher the water holding capacity. WC showed a fluctuating downward trend with the increasing distance to forestland while increased with the distance to waters, which was affected by the complex spatial pattern of mountains and rivers. Meanwhile, WC showed a general decreasing trend of decreased first and then increased or unchanged with the distance to grassland and cultivated land. WC had a linear positive correlation with the distance to construction land, whereas had an opposite correlation with population and GDP. With the increase of human activity intensity, a large amount of natural vegetation was converted into impervious surface [57], which directly led to the decrease of WC.

4.2. Scale Effect of ESs and Their Complex Interactions

Our study examined spatial-temporal heterogeneity and complex interactions of four key ESs, i.e., BI, CF, SC and WC. These ESs not only can reflect the key ecological problems facing the CUACM [85], but also have great contributions to regional welfare and sustainable development [24]. At grid scale, there was obvious spatial agglomeration of ESs in mountainous areas. Numerous studies have proved that multi-dimensional topographic factors in mountainous areas not only controlled the spatial distribution pattern of solar radiation and precipitation, but also affected the spatial heterogeneity of soil properties [86,87]. Characterized by diverse vegetation types and complex vegetation structures, mountainous areas not only provided high CF, but also provided various suitable habitats for species [88]. The mountainous areas, with high elevation, steep slope and high vegetation coverage, had a high potential soil erosion, low actual soil erosion and a low runoff that contributed greatly to SC and WC [61]. At the administrative districts scale, ESs were relatively higher in Beibei, Banan and Yubei, and drastically decreased in Jiangbei, which were closely related to regional vegetation coverage and urbanization development. Banan, Beibei and Yubei had relatively high average vegetation coverage, which were 0.86, 0.79 and 0.78, respectively. In 2008, according to CEFR, Banan, Beibei and Yubei had been designated as important ecological barriers of the study area. Affected by the development of Liangjiang New Area and the rapid expansion of urbanization, the ecological land in Jiangbei district was occupied by construction land, which directly led to the reduction of ESs. This was consistent with the research of Wang et al. [2], showing dramatically decreased in ESs with urban land expansion. It is worth noting that the spatial congruence of WC was relatively low from 2000 to 2018, which was mainly affected by the interannual variation of precipitation and potential evapotranspiration. Although the precipitation decreased first and then increased (992.89 mm, 980.42 mm and 1126.55 mm), the magnitude of its change did not exceed the increases in evapotranspiration (778.08 mm, 769.09 mm and 1364.77 mm). Growing studies has been discussing that the change of WC is not completely determined by land use change, but dominated by regional climate change and human activities [70,89]. At present, most studies have primarily focused on the spatiotemporal pattern analysis of ESs [90], but the scale-dependent effects have usually been ignored, such as land use and topographic heterogeneity. The spatiotemporal heterogeneity of ESs at different scales can be more clearly characterized by hotspots identification [2]. In our study, it can be clearly seen that multiple ES hotspots were mainly composed of forestland and distributed in high TNI with high altitude and steep slope, which proved the spatial-temporal patterns of ESs mentioned above.

ESs were dependent on dynamic processes that act and interact at different scales. Therefore, the interactions between ES pairs had obvious scale effects. Yang et al. [28] found that the degree of trade-offs between ESs tended to weaken gradually as the scale increased, and there were even signs that it would turn into synergies. Again, remarkably similar results were obtained in our study. In the entire CUACM, there was a positive relationship within ESs, indicating synergy or at least no conflicts. A possible explanation was that they shared common driving factors, such as land use, TNI and NDVI. Specifically, the forestland located in high TNI with high elevation, slope and NDVI not only provided suitable habitats for species, but also had the highest potential for CF [81]. The dense canopy, thick litter and developed root systems had an effective capability of rainfall interception, slope stability, and soil and water conservation [24]. However, there was distinct spatial heterogeneity in trade-offs and synergies between ES pairs at grid scale. This can be explained in two ways; on one hand, the natural and social driving factors had distinct geographical differences, on the other hand, each ES had different dominant driving factor. For example, SC exhibited a strong dependence on geomorphological factors, while WC was extremely sensitive to meteorological factors [70]. Bivariate LISA and ESBs had the potential to objectively divide study areas into different groups and make spatial explicit mapping, which can provide novel insights into the understanding and managing of multiple ESs.

4.3. Implications for ESs Management and Urban Planning

Detecting the complex relationships among multiple ESs was crucial for ecosystem management and urban planning. The spatial concordance and segregation among different ESs hotspots coincided with the formation of different types of ESBs, which can provide guidance for ecological function regionalization. Take 2018 as an example, ESB1 delivered 89.09% of class 4 ES hotspots and 50.95% of class 3 ES hotspots (Figure S4a). The multiple ES hotspots occupied 95.19% of the total area of ESB1 (Figure S4b). ESB2 covered 61.57%, 41.59% and 4.72% of class 2, class 3 and class 4 ES hotspots, respectively, with the proportion of multiple ES hotspots was up to 97.09%. This meant that both ESB1 and ESB2 could support multiple ESs simultaneously, which should be determined as the priority areas for ecological protection to ensure regional sustainable development [2]. ESB3 were mainly consisted of non-hotspots and class 1 ES hotspots, covering 60.72% and 29.77%, respectively. Considering that ESB3 was mainly located in the piedmont transition zone, with high proportion of cultivated land and high SC, land consolidation, reforestation and afforestation were urgently needed to enhance the water and soil conservation efficiencies [3]. ESB4 and ESB5 were primarily covered by non-hotspots, accounting for 82.37% and 92.97%. ESB4 was fringe areas of urban sprawl and provide an important buffering function in ecological environmental protection. Ecological protection and restoration projects using nature-based solutions should be implemented to address regional environmental challenges [91]. In terms of the lowest ESs in urban ecosystems that were in ESB5, the creation and protection of green spaces should be encouraged to relieve ecological pressures, whereas the disordered expansion of cities and the occupation of ecological land for construction should be restricted [33]. For example, as the Yangtze River and Jialing River crossed the urban center, establishing protected greenbelts along the rivers should be considered to reduce human disturbances. The ESBs were consistent with the classification of the major function-oriented zone planning [92], which can provide a theoretical reference for decision-making related to ecological protection and urban development. ESB1 and ESB2 may be considered as prohibited development areas, ESB3 as restricted development areas, ESB4 as prioritized development areas, and ESB5 as optimized development areas.

In general, the ESs synergies areas reflected the realizability of ESs multifunctional management and could be used as the pilot areas to solve trade-offs and encourage the formation of ESBs with more ESs hotspots. In synergy areas with the most ESs hotspots, the original ecosystem and vegetation should be strictly protected. The ESs trade-off areas illustrated the severity of trade-offs and should be the priority areas for ecological restoration. In ESs trade-off areas, different ecological restoration measures should be taken according to the different characteristics of trade-offs. For example, the intensity of development and utilization should be properly controlled to prevent construction land from occupying ecological land with high ESs. Implementing the project of returning farmland to forestland and grassland to increase vegetation coverage and maintain biodiversity are also effective means of regional ESs management.

4.4. Applications of Remote Sensing for ESs Evaluation

Remote sensing has the potential of performing synoptic, spatially continuous, and frequent observations at a wide range of locations [93], which can provide quantitative and spatially explicit thematic data regarding various biophysical characteristics that are often spatialized for ESs evaluation [94,95]. In this study, land use based on Landsat interpretation was used as a necessary input to the InVEST model to evaluate ESs. Land use has been widely used as a dominant proxy for ESs [96], such as forest types and cover are useful to estimate raw materials provisioning [97]. Galbraith et al. [98] stated that forest and other natural areas were positively related to the provision of pollination services. De Araujo Barbosa et al. [99] systematically reviewed the literatures on the evaluation of different ESs using remote sensing and pointed out that land cover was the most important proxy variable in all ESs assessments, followed by NDVI. NDVI was taken as an important parameter of cover-management factor in SC evaluation in this study. Extensive studies have demonstrated that vegetation indices derived from remote sensing inversion are an increasingly major data source for ESs evaluation [99]. For example, net primary production can be directly assessed by remote sensors such as MODIS [95]. It is proved that the vegetation indices based on remotely sensed spectral characteristics can also respond to different biodiversity conditions and serve as the basis for cultural services evaluation [18,100]. The geomorphological data obtained from ASTER GDEM V2 digital elevation model was also the key remote sensing information for BI, SC and WC evaluations in this study, as mapping topographic heterogeneity was helpful to reveal the possibility of providing biological refugia and the potential of soil and water loss [101]. Various remotely sensed information can also be used in combination as proxies for kinds of variables which in turn can be used as proxies for ESs [99,102]. Integrating land use, vegetation indices, meteorological, geomorphological, soil characteristics and other reference factors, we assessed BI, CF, SC and WC by InVEST model.

Generally speaking, remote sensing can be used to evaluate ESs in three different ways: direct and indirect measures and in combination with ecosystem models. Most provisioning services are tangible goods that can be quantified by the direct use of remotely sensed radiation information [103]. For example, the production capacity of agricultural and forestry ecosystems can be quantified based on satellite-derived chlorophyll measurements [96]. Since the regulating services are derived from complex ecological functions or processes, they are primarily assessed based on proxy variables reflecting various ecosystem conditions [104]. Taking climate regulation services as an example, it can be estimated by the net ecosystem exchange of CO₂ flux, which can be quantified using biomass as proxies [99,105], just as we evaluated CF. Cultural services arise from people’s perception of the natural phenomenon and environment, which are often evaluated by comprehensively considering biophysical characteristics of ecosystems and beneficiaries’ perception [106]. The combination of remotely sensed inversion data in a multi-agent modelling environment has the potential to provide necessary indicators for cultural services evaluation [99,107]. Supporting services are necessary for the maintenance of ecosystems and the production of all other ESs, which are usually estimated by related proxy variables [108]. Maintaining animal and plant diversity is one of the most important supporting services, which can be indirectly evaluated by remote sensing-based habitat mapping, individual spatial distribution mapping, spectral diversity, etc. [109]. In this study, habitat quality was assessed by the InVEST model, and the dimensionless relative value of 0–1 could be used as a surrogate to quantify biodiversity [61]. Although various methods and indicators have been developed to quantify and map ESs, there are still some limitations in the assessment of some ESs, such as cost-intensive and time-consuming [82]. Therefore, considering the main ecological problems faced by the study area and the relative importance of ESs to human well-being, four key ESs were evaluated in this study. Limited by the accessibility and availability of data, there was some inconsistency in the precision and resolution of the dataset, which had become a common problem in ESs evaluation [49].

5. Conclusions

This study demonstrated the effectiveness of the integrated method of InVEST model, the least-squares regression model and hotspots analysis for ESs assessment, spatial–temporal variation and hotspots identification. Based on spearman correlation coefficients, bivariate LISA and K-means clustering, the complex interactions among multiple ESs were detected. GAM was used to discuss the linear or nonlinear relationships between ESs and driving factors. The results showed that: (1) BI, SC and WC displayed decreased trends, while CF tended to remain basically unchanged. High values of ESs were concentrated in mountainous areas, especially in Jinyun, Zhongliang, Tongluo, Mingyue, Huaying mountains. In addition, the Yangtze River, Jialing River and other water bodies also had higher BI. At administrative districts scale, there were relatively higher ESs in Beibei, Banan and Yubei, and the dramatic decreases were observed in Jiangbei, which were closely related to regional vegetation coverage and urbanization development. (2) ES hotspots had distinct spatial heterogeneity in different land use types and TNI gradients. From the land use types perspective, forestland and grassland were major contributors to ESs, with the average area proportion of multiple ES hotspots accounting for 94.68% and 84.02%, respectively. Generally, the distribution index of multiple ES hotspots presented a tendency to increase with an increase in TNI gradients, which indicated that ESs hotspots were predominantly distributed in mountainous areas with high altitude and steep slope. (3) In the entire study area, all six pairs of ESs exhibited synergistic relationships. The correlation coefficients among BI, CF and SC showed a tendency to increase while the synergies between WC and other ESs were diminished. At grid scale, the synergies were mainly concentrated in the high-high and low-low clusters, i.e., mountainous areas and urban central areas. Five ESBs were determined according to the aggregation of multiple ESs on spatial and temporal scales, which can provide theory supports for land management decision-making and urban planning. (4) ESs were primarily driven by natural environment and socioeconomic factors. Land use types had profound impacts on BI and CF; in particular, forestland and grassland showed a clear positive affect. SC was mainly determined by geomorphological factors, while WC was more likely to be explained by precipitation. The unique and crisscrossing landscapes of mountains and rivers resulted in a strongly nonlinear relationship between ESs and the distance factors. Furthermore, policy factors were confirmed to have a certain positive effect on ESs. This study not only provides a scientific basis for the maximization of ESs benefits and multi-objective collaborative management, but is also beneficial to realize the coordinated and sustainable development of regional urban expansion and ecological environment.