Investigating Ecosystem Service Trade-Offs/Synergies and Their Influencing Factors in the Yangtze River Delta Region, China

A comprehensive understanding of the ecosystem services (ESs) trade-off/synergy relationships has become increasingly important for ecological management and sustainable development. This study employed the Yangtze River Delta (YRD) region in China as the study area and investigated the spatiotemporal changes in three ESs, namely, carbon storage (CS), water purification (WP), and habitat quality (HQ). A trade-off/synergy degree (TSD) indicator was developed that allowed for the quantification of the trade-off/synergy intensity, and the spatial pattern of the TSD between ESs in the YRD region to be analyzed. Furthermore, a geographically weighted regression (GWR) model was used to analyze the relationship between the influencing factors and trade-offs/synergies. The results revealed that CS, WP, and HQ decreased by 0.28%, 2.49%, and 3.38%, respectively, from 2005 to 2015. The TSD indicator showed that the trade-off/synergy relationships and their magnitudes were spatially heterogeneous throughout the YRD region. The coefficients of the natural and socioeconomic factors obtained from the GWR indicated that their impacts on the trade-offs/synergies vary spatiotemporally. The impact factors had both positive and negative effects on the trade-offs/synergies. The findings of this study could improve the understanding of the spatiotemporal dynamics of trade-offs/synergies and their spatially heterogeneous correlations with related factors.


Introduction
Ecosystem services (ESs) refer to the products and services obtained from ecosystems. ESs are the basis of human survival and are closely related to human wellbeing. As evidenced by the Millennium Ecosystem Assessment (MEA), the loss of ESs in recent decades directly threatens regional and global ecological security [1], and this worsening trend is likely to continue in the future due to the increasing intensity of human activities [2].
In the process of pursuing socioeconomic benefits, human beings have significant impact on the ESs interactions, which is manifested as the ES trade-off and synergy. The ES trade-off presents a win-lose situation that involves enhancing one ecosystem service at the expense of reducing another, and ES synergy presents a win-win situation that involves the enhancement of multiple ecosystems [3]. The dramatic socioeconomic development leads to an improvement in one ES at the expense of other ESs, thus threatening the stability and security of ecosystems and significantly influencing human wellbeing [4]. Ecosystem management relies on "our best understanding of the ecological interactions and processes necessary to sustain ecosystem composition, structure, and function" [5]. Ecosystem management should not only pursue single ecosystem service benefits but also consider the balance among multiple ESs to maximize the overall benefits and promote regional sustainable development. The importance and urgency of understanding the complex interactions among multiple ESs have been widely recognized by ecosystem managers and Land 2022, 11, 106 3 of 22 (RDA) [4], and multivariate regression trees [31] were applied to explain the global relationships among ES trade-offs/synergies and the related factors without considering spatial heterogeneity. It is recognized that the trade-off/synergy relationships among ESs are affected by spatial heterogeneity. For example, in South Africa, a trade-off between carbon sequestration and freshwater supply was found [34], while research on the Baiyangdian River Basin found a synergistic relationship between these two variables [35]. Therefore, it was assumed that the effects of influencing factors on trade-offs/synergies should be spatial heterogeneity, i.e., the variations in the effects over space. Global relationships reveal information only about the average conditions and might consequently overlook location-specific impacts [36,37]. If we ignore spatial heterogeneity, errors might arise in the results of analyses on the impacts of influencing factors on ES interactions, increasing uncertainty about ES management policies [29]. Despite the breadth of previous studies, a key question remains to be answered: how do different natural and socioeconomic factors impact the trade-offs/synergies among ESs with consideration of spatial heterogeneity? The geographical weighted regression (GWR) model is a local regression method that can explore the spatially varying relationship [38,39]. This method could provide more detailed spatial information on the complicated relationship between multiple variables and trade-offs/synergies among ESs.
The YRD region is an economically developed and densely populated area in China that is experiencing serious ecological degradation and environmental pollution. This paper aims at partly address this knowledge gap by using the YRD as the study area and investigating the spatiotemporal changes in the trade-offs/synergies among ESs and their spatial heterogeneous correlation with the related factors. Specifically, the objectives were to: (1) reveal the dynamics of the trade-offs/synergies that exist among ESs in the YRD region by using a new proposed indicator; (2) explore the spatially heterogeneous relationships that exist between the influencing factors and trade-offs/synergies by employing the GWR model; and (3) propose the implications of coordinating ESs to achieve sustainable development for ecosystem management.

Study Area
In this paper, the YRD region (115 • 47 -122 • 56 E, 27 • 55 -34 • 27 N) refers to associated urban agglomeration discussed in the "Yangtze River Delta Urban Agglomeration Development Plan" approved by the State Council in May 2016. The YRD region is located in a subtropical monsoon climate zone, with an annual precipitation of 800-2050 mm, presenting significant regional differences. The precipitation decreases from southeast to northwest, and the annual average temperature ranges from 9.3-17.3 • C. As shown in Figure 1, the YRD region includes 26 cities, including Shanghai, Nanjing, Hangzhou, Hefei, etc. The area of the whole region is 211,700 km 2 , accounting for 2.1% of the total land area of China.
The YRD region is one of the regions with the most rapid economic and urbanization development in China. Changes in land-use have significantly changed the structure and pattern of regional ecosystems, and many ecological environmental risks related to water regulation, air purification, and climate regulation are becoming increasingly severe. The contradiction between regional socioeconomic development and the protection of the ecological environment has become the main limiting factor that restricts the coordinated development of the YRD region and has a certain negative impact on human wellbeing. Due to socioeconomic development and the expansion of cities, the magnitude of ES trade-offs will continue to increase.
Therefore, the quantitative measurement of the spatiotemporal dynamics in the patterns of the trade-offs/synergies in the YRD and the identification of the spatially heterogeneous effects of influencing factors on trade-offs/synergies in the YRD region may not only enrich research on ESs but also provide a scientific basis for optimizing ecosystem management and regional planning. The YRD region is one of the regions with the most rapid economic and urbanization development in China. Changes in land-use have significantly changed the structure and pattern of regional ecosystems, and many ecological environmental risks related to water regulation, air purification, and climate regulation are becoming increasingly severe. The contradiction between regional socioeconomic development and the protection of the ecological environment has become the main limiting factor that restricts the coordinated development of the YRD region and has a certain negative impact on human wellbeing. Due to socioeconomic development and the expansion of cities, the magnitude of ES trade-offs will continue to increase.
Therefore, the quantitative measurement of the spatiotemporal dynamics in the patterns of the trade-offs/synergies in the YRD and the identification of the spatially heterogeneous effects of influencing factors on trade-offs/synergies in the YRD region may not only enrich research on ESs but also provide a scientific basis for optimizing ecosystem management and regional planning.

Data
The following data were used in this study: (1)

Data
The following data were used in this study: (1) land cover data on the YRD region for 2005 and 2015 with a resolution of 30 m, obtained from the Data Centre for Resources and Environmental Sciences, Chinese Academy of Sciences (RESDC) (http://www.resdc. cn/Default.aspx, accessed on 7 June 2021). Land cover was reclassified into six categories based on the data: cultivated land, forestland, grassland, water bodies, built-up land, and unused land. (2) Meteorological data, including information on precipitation, total solar radiation and the temperature in the YRD region from 2005 to 2015, derived from the China Meteorological Data Service Centre (http://data.cma.cn/, accessed on 13 July 2021).
(3) Soil property data, including soil depth and the content percentage of sand, clay, silt and organic matter at the scale of 1:1,000,000, obtained from the Harmonized World Soil Database version 1.1 (HWSD) (accessed on 13 July 2021). (4) Digital elevation model (DEM) data for the YRD region with a resolution of 30 m, obtained from the Geospatial Data Cloud (http://www.gscloud.cn/, accessed on 14 July 2021). (5) Socioeconomic data, including information on the gross domestic product (GDP) and population of the YRD region in 2005 and 2015, obtained from the statistical yearbook of 26 cities in the YRD region (accessed on 21 May 2021). In this study, spatial data were projected to the WGS_1984_Albers coordination system, and all raster data were resampled to a resolution of 100 m.
This study was divided into the following steps: (1) the estimation of ESs; (2) the measurement of the trade-offs/synergies among ESs during the study period; and (3) the estimation of the correlation between trade-off/synergy and multiple related factors. The method used in this study is presented in the following sections.

Estimation of ESs
In this study, we identified three key types of ESs with consideration of the following criteria: (1) the selected ESs should be strongly related to human wellbeing in the YRD region. (2) The data needed to calculate the selected ESs should be available. (3) The selected ESs should be significantly influenced by human activities and socioeconomic development. (4) The coordination of the selected ESs should be important for regional sustainable development. Based on these criteria, two key ESs (carbon storage, CS; water purification, WP) and an ES indicator (habitat quality, HQ) were selected for the YRD region. They are significantly related to human wellbeing and are important for the ecological environment. Their importance increases in the YRD region due to rapid socioeconomic development. Figure 2 shows the cascade model of ES generation and valuation [40,41]. The model demonstrates overlap between the ecosystem on the supply-side and human benefit and value on the demand side, with ecosystem service located at the intersection of the two. In this study, a biophysical approach was adopted to value the ecosystem service on the supply-side. method used in this study is presented in the following sections.

Estimation of ESs
In this study, we identified three key types of ESs with consideration of the following criteria: (1) the selected ESs should be strongly related to human wellbeing in the YRD region. (2) The data needed to calculate the selected ESs should be available. (3) The selected ESs should be significantly influenced by human activities and socioeconomic development. (4) The coordination of the selected ESs should be important for regional sustainable development. Based on these criteria, two key ESs (carbon storage, CS; water purification, WP) and an ES indicator (habitat quality, HQ) were selected for the YRD region. They are significantly related to human wellbeing and are important for the ecological environment. Their importance increases in the YRD region due to rapid socioeconomic development. Figure 2 shows the cascade model of ES generation and valuation [40,41]. The model demonstrates overlap between the ecosystem on the supply-side and human benefit and value on the demand side, with ecosystem service located at the intersection of the two. In this study, a biophysical approach was adopted to value the ecosystem service on the supply-side. The InVEST model is an open source software that has been widely used for calculating and mapping ecosystem services. It is appropriate for use in contexts where ecological processes are well understood. By using InVEST, the spatial patterns of CS, WP, and HQ were quantitatively estimated. The calculations used for the parameters required by the model are shown in Table 1. A carbon storage and sequestration model was used to estimate CS based on four types of carbon pools (aboveground biomass, belowground biomass, the soil and dead organic matter). WP was evaluated in terms of nitrogen export by adopting a nutrient delivery ratio model. A higher nitrogen export value indicates a lower level of WP. HQ was calculated using the habitat quality (HQ) model in InVEST, which considers land cover types and threats [42]. The InVEST model is an open source software that has been widely used for calculating and mapping ecosystem services. It is appropriate for use in contexts where ecological processes are well understood. By using InVEST, the spatial patterns of CS, WP, and HQ were quantitatively estimated. The calculations used for the parameters required by the model are shown in Table 1. A carbon storage and sequestration model was used to estimate CS based on four types of carbon pools (aboveground biomass, belowground biomass, the soil and dead organic matter). WP was evaluated in terms of nitrogen export by adopting a nutrient delivery ratio model. A higher nitrogen export value indicates a lower level of WP. HQ was calculated using the habitat quality (HQ) model in InVEST, which considers land cover types and threats [42].

Method Formula and Variable Description
Carbon Storage (CS) InVEST carbon storage and sequestration model CS = C above + C below + C soil + C dead C above : above ground biomass carbon stocks; C below : below ground biomass carbon stocks; C soil : soil carbon stocks; C dead : dead organic matter carbon stocks. The biomass carbon density for the different land-use types was derived from the results of a relevant study [43] (Table S1) Water Purification (WP) InVEST nutrient delivery model ALV x : adjusted nitrogen export value for pixel x; HSS x : hydrological sensitivity score for pixel x; pol x : export coefficient for pixel x. λ x : runoff index for grid cell x; λ w : mean runoff index in the watershed; ∑ U Y U : sum of the water yield of grid cell along the flow path above grid cell x. The biophysical parameters for water purification was presented in Table S2.
Habitat Quality (HQ) InVEST HQ model Q xj : HQ of land − use type j; H j : the habitat suitability of land − use type j; D xj : total threat level in a grid cell x with land − use type j; r y : intensity of the threat within the cell y; w r : relative impact of threat r; β x : level of accessibility in grid cell x; S jr : sensitivity of land − use type j to threat r; i rxy : impact of threat r from grid cell y on the habitat in grid cell x; d xy : linear distance between grid cells x and y; d rmax : maximum effective distance of the threat. The parameters used in the model were presented in Table S3. HQ is a measure of the quality of the ecological infrastructure.

Measurement of the Trade-Offs/Synergies among ESs
The trade-off describes the situation in which the increase of one ES leads to a decrease in another, and synergy describes the situation in which both ESs either increase or decrease simultaneously [44]. Spearman correlation analysis was employed to explore the relationships among ESs at the county scale based on changes in the ESs from 2000 to 2015 [45]. If the Spearman correlation coefficient for two services is higher than 0, then a synergistic relationship exists; otherwise, a trade-off relationship exists.
An innovative indicator, the trade-off/synergy degree (TSD), was proposed and applied to further explore the magnitudes and spatial patterns of trade-off/synergy relationships among the ESs. If the value of TSD is smaller (or greater) than 0, it indicates that there is a trade-off (or synergy) between the paired ESs. Its absolute value reflects the magnitude of the trade-off/synergy relationship.
The TSD indicator for ES i and ES j (TSD i−j ) over the periods of t2 and t1 was calculated with Equation (1):

Measurement of the Trade-Offs/Synergies among ESs
The trade-off describes the situation in which the increase of one ES leads to a decrease in another, and synergy describes the situation in which both ESs either increase or decrease simultaneously [44]. Spearman correlation analysis was employed to explore the relationships among Ess at the county scale based on changes in the Ess from 2000 to 2015 [45]. If the Spearman correlation coefficient for two services is higher than 0, then a synergistic relationship exists; otherwise, a trade-off relationship exists.
An innovative indicator, the trade-off/synergy degree (TSD), was proposed and applied to further explore the magnitudes and spatial patterns of trade-off/synergy relationships among the Ess. If the value of TSD is smaller (or greater) than 0, it indicates that there is a trade-off (or synergy) between the paired ESs. Its absolute value reflects the magnitude of the trade-off/synergy relationship.
The TSD indicator for and ( ) over the periods of 2 and 1 was calculated with Equation (1): where ∆ , and ∆ , are the relative changes in ESs i and j over the period t2-t1, respectively, and are calculated as follows: where , and , represent the value of at time points 2 and 1, respectively. As shown in Equation (1), if ∆ , × ∆ , = 0, then = 0, implying that no trade-off/synergy relationship exists between and . If ∆ , × ∆ , > 0, a synergy relationship exists between and , and the level of synergy can be measured by ((∆ , ) + (∆ , ) ) 2 ⁄ . If ∆ , × ∆ , < 0, then a trade-off relationship exists between the two ESs, and the level of the trade-off can be measured by − ((∆ , ) + (∆ , ) ) 2 ⁄ . Before conducting a spatial regression analysis, the global Moran's I index was first applied to detect the spatial dependence of the trade-offs/synergies among the ESs. At a given significance level, Moran's I value ranges from −1 to 1. A value closer to 1 implies that the trade-off/synergy relationship in a specific area has a trend similar to that of the surrounding areas. A value closer to −1 suggests that the trade-off/synergy relationship in a specific area has a trend that is dissimilar to that of the surrounding areas. If Moran's I is 0, then the trend of the trade-off/synergy relationship is randomly distributed [46].

Evaluation of the Impacts of Factors on ESs Trade-Offs/Synergies
With consideration of natural and socioeconomic influences, six potential influencing factors were selected ( (1) where ∆ES i,t2−t1 and ∆ES j,t2−t1 are the relative changes in ESs i and j over the period t2-t1, respectively, and are calculated as follows: where ES i,t2 and ES i,t1 represent the value of ES i at time points t2 and t1, respectively. As shown in Equation (1), if ∆ES i,t2−t1 × ∆ES j,t2−t1 = 0, then TSD i−j = 0, implying that no trade-off/synergy relationship exists between ES i and ES j . If ∆ES i,t2−t1 × ∆ES j,t2−t1 > 0, a synergy relationship exists between ES i and ES j , and the level of synergy can be measured relationship exists between the two ESs, and the level of the trade-off can be measured by Before conducting a spatial regression analysis, the global Moran's I index was first applied to detect the spatial dependence of the trade-offs/synergies among the ESs. At a given significance level, Moran's I value ranges from −1 to 1. A value closer to 1 implies that the trade-off/synergy relationship in a specific area has a trend similar to that of the surrounding areas. A value closer to −1 suggests that the trade-off/synergy relationship in a specific area has a trend that is dissimilar to that of the surrounding areas. If Moran's I is 0, then the trend of the trade-off/synergy relationship is randomly distributed [46].

Evaluation of the Impacts of Factors on ESs Trade-Offs/Synergies
With consideration of natural and socioeconomic influences, six potential influencing factors were selected ( Understanding how the influencing factors have contributed to the variations in the trade-offs/synergies among ESs is crucial for effective ecosystem management. Multiple linear regression (MLR) is one of the regression methods commonly used for exploring the relationships between multiple explanatory factors and independent variables [47,48]. Taking the study area as a whole, MLR can be used to explore the average regression coefficient by identifying a global relationship. The MLR used in this study is expressed in Equation (3): where y i is the estimated TSD value; β 0 and β k represent the intercept and regression coefficients of independent variable k, respectively; x ik is the value of the independent variable k at sample i; and ε is the random error term for sample i. Understanding how the influencing factors have contributed to the variations in the trade-offs/synergies among ESs is crucial for effective ecosystem management. Multiple linear regression (MLR) is one of the regression methods commonly used for exploring the relationships between multiple explanatory factors and independent variables [47,48]. Taking the study area as a whole, MLR can be used to explore the average regression coefficient by identifying a global relationship. The MLR used in this study is expressed in Equation (3): where is the estimated TSD value; and represent the intercept and regression coefficients of independent variable , respectively; is the value of the independent variable at sample ; and ε is the random error term for sample . A geographically weighted regression (GWR) model can estimate a location-specific regression coefficient for each unit in the study area by considering spatially varying re- A geographically weighted regression (GWR) model can estimate a location-specific regression coefficient for each unit in the study area by considering spatially varying relationships that exist between the dependent and independent variables. In this study, a GWR model was adopted to investigate the relationships between the related factors and the trade-offs/synergies while considering spatial heterogeneity. The formula used for the GWR is similar to that used for the MLR, except the former takes spatial heterogeneity into consideration. The GWR can be expressed as follows [49]: where y i represents the TSD value for sample i; (u i , v i ) represents the coordinates of sample i; β 0 (u i , v i ) is the constant term estimated for sample i; and β k (u i , v i ) is the estimated regression coefficient of factor k at sample i. x ik represents the value of the independent variable k that could explain TSD at sample i; ε i is the random error term for sample i. The Gaussian kernel was selected to calculate the spatial weighting [50]: where w ij is the weight of sample j for sample i. d ij represents the Euclidean distance between sample i and sample j. h is the kernel bandwidth, which determines the range of spatial dependency. In this study, we selected the optimal bandwidth identified by the Akaike's information criterion (AIC) method for the purposes of minimizing the AIC value: the estimated results with lower AIC values better reflect reality. The performance of the regression models (MLR and GWR) can be measured and compared by considering R 2 and AIC. R 2 is the square of the correlation coefficient, indicating the degree of agreement between the estimated value and the observed value. AIC is not an absolute goodness of fit, but it can be used to effectively compare several different models. A difference between the AIC values of two models higher than 3 suggests that there is a significant difference between the two models. A more suitable model, which is assessed by its capability to better explain the variations in the dependent variables, tends to have a higher R 2 and lower AIC.

Spatial Patterns of ESs
As shown in Figure 4, forestland is mainly located in the southern part of the YRD region. High-intensity built-up land is concentrated along the Yangtze River. Cultivated land is mainly distributed in the northern part of the YRD region. As shown in Table 2, cultivated land accounted for 48.04% of the total area, the largest proportion of the YRD region, in 2015. Built-up land accounted for 11.90% of the total area in 2015.
independent variable that could explain TSD at sample ; is the random error term for sample . The regression coefficients of each sample are estimated by weighting all units around sample . The units closer to sample have a stronger impact on the estimated coefficient. The Gaussian kernel was selected to calculate the spatial weighting [50]: where is the weight of sample for sample . represents the Euclidean distance between sample and sample . ℎ is the kernel bandwidth, which determines the range of spatial dependency. In this study, we selected the optimal bandwidth identified by the Akaike's information criterion (AIC) method for the purposes of minimizing the AIC value: the estimated results with lower AIC values better reflect reality.
The performance of the regression models (MLR and GWR) can be measured and compared by considering R 2 and AIC. R 2 is the square of the correlation coefficient, indicating the degree of agreement between the estimated value and the observed value. AIC is not an absolute goodness of fit, but it can be used to effectively compare several different models. A difference between the AIC values of two models higher than 3 suggests that there is a significant difference between the two models. A more suitable model, which is assessed by its capability to better explain the variations in the dependent variables, tends to have a higher R 2 and lower AIC.

Spatial Patterns of ESs
As shown in Figure 4, forestland is mainly located in the southern part of the YRD region. High-intensity built-up land is concentrated along the Yangtze River. Cultivated land is mainly distributed in the northern part of the YRD region. As shown in Table 2, cultivated land accounted for 48.04% of the total area, the largest proportion of the YRD region, in 2015. Built-up land accounted for 11.90% of the total area in 2015.   Along with land cover changes in the YRD region, the selected ESs also varied over the study period. As shown in Table 3, the overall magnitude of CS, WP (measured by nitrogen export), and HQ in the YRD region decreased from 2005 to 2015. CS decreased marginally from 5.028 × 10 8 to 5.014 × 10 8 Mg. The total amount of WP decreased from 3.891 × 10 7 kg to 3.794 × 10 7 kg. The average value for HQ shows a declining trend, with the value decreasing from 0.385 in 2005 to 0.372 in 2015. There is significant spatial heterogeneity in CS, WP, and HQ in 2005 and 2015. As shown in Figure 5, the high-provision areas for CS were mainly distributed in the southern mountainous area of the YRD, while the lower CS values were located in urban areas, where a large amount of non-built-up land had changed to built-up land. The high-provision areas for WP were mainly distributed in the northern area. WP in the YRD shows a gradual decreasing trend from southeast to northwest, which is consistent with the spatial patterns of the proportion of cultivated land. High values of HQ were found in the south areas, which are at higher elevations, while low value areas were mainly distributed in the north area. This result can be explained by the fact that areas with high elevations are mountainous with dense vegetation and rich biodiversity. In these areas, the natural environment, such as mountains, face fewer threats.
The majority of counties (157 out of 160 counties) in the YRD region show a decreasing trend in CS from 2005 to 2015 ( Figure 6). Among these cities, Kunshan (Suzhou city) and Jiading (Shanghai city) experienced the most significant decline in average CS, decreasing by 17.831 Mg/ha and 15.236 Mg/ha, respectively. Dafeng, Xiaoshan, and Chongming showed a slight increase in average CS over the study period.
Overall, WP exhibited a downward trend across the YRD region. At the county level, over one-half (89) of the counties in the YRD region showed that WP is increasing, mainly in the eastern part of the YRD, accounting for 52.8% of the total area. Cultivated land was converted to BUL. The counties exhibiting a significant decrease in WP are in the southeastern part of the YRD, including Wenling, Yuhuan, and Huangyan.
Most of the counties (156) presented a declining trend for HQ, especially the Suzhou city core and Xiaoshan (Hangzhou city), decreasing by 33.3% and 18.3%, respectively. Only 4 counties experienced an increase from 2005-2015, mainly centralized in Nantong, Chuzhou, Yangzhou, and Taaizhou cities, accounting for 4.43% of the total area. The reduction in HQ was mainly due to significant urban expansion.  Overall, WP exhibited a downward trend across the YRD region. At the county level, over one-half (89) of the counties in the YRD region showed that WP is increasing, mainly in the eastern part of the YRD, accounting for 52.8% of the total area. Cultivated land was converted to BUL. The counties exhibiting a significant decrease in WP are in the southeastern part of the YRD, including Wenling, Yuhuan, and Huangyan.   Overall, WP exhibited a downward trend across the YRD region. At the county level, over one-half (89) of the counties in the YRD region showed that WP is increasing, mainly in the eastern part of the YRD, accounting for 52.8% of the total area. Cultivated land was converted to BUL. The counties exhibiting a significant decrease in WP are in the southeastern part of the YRD, including Wenling, Yuhuan, and Huangyan.

Spatiotemporal Dynamics of ESs Trade-Offs/Synergies
Spearman correlation analysis was employed to analyze the trade-off/synergy relationships between ES pairs. Both trade-offs and synergies were observed ( Figure S1). We found that all three pairs of ESs (CS-WP, CS-HQ, WP-HQ) showed significant positive or negative correlations (p < 0.01). CS was significantly positively correlated with HQ. This result was supported by our findings that CS increased from 2005 to 2015, and HQ decreased over this period. The trade-offs between CS/HQ and WP/HQ suggested that a decrease or increase in one ES could lead to an increase or decrease in another ES in a specific county. CS had a strong and significant synergy with HQ, with a Spearman correlation coefficient larger than 0.7.
The TSD indicator was further applied to reveal the spatial distribution and magnitude of the trade-off/synergy relationships among the ESs. The TSD indicator is an extension of the traditional correlation analysis because it can detect local rather than global relationships between paired ESs. The relationships between paired ESs varied over space (Figure 7).
The TSD indicator was further applied to reveal the spatial distribution and magnitude of the trade-off/synergy relationships among the ESs. The TSD indicator is an extension of the traditional correlation analysis because it can detect local rather than global relationships between paired ESs. The relationships between paired ESs varied over space (Figure 7). The negative TSD value for CS-WP were mainly distributed in the eastern and southwestern areas, which indicates that a decrease in CS could lead to an increase in WP in these regions. A larger amount of ecological land (forestland and grassland), the main type of land that supports CS, was occupied by built-up land due to the rapid development in this region. Significant urban expansion resulted in a decrease in CS. The loss in ecological land led to an increase in nitrogen export. Furthermore, the TSD values for CS-WP are lower in eastern and central regions and higher in southern and northwestern regions. The relationship between CS and WP is synergetic (TSD > 0) in the southeastern and northwestern areas, and this result contrasts with the implications of the Spearman The negative TSD value for CS-WP were mainly distributed in the eastern and southwestern areas, which indicates that a decrease in CS could lead to an increase in WP in these regions. A larger amount of ecological land (forestland and grassland), the main type of land that supports CS, was occupied by built-up land due to the rapid development in this region. Significant urban expansion resulted in a decrease in CS. The loss in ecological land led to an increase in nitrogen export. Furthermore, the TSD values for CS-WP are lower in eastern and central regions and higher in southern and northwestern regions. The relationship between CS and WP is synergetic (TSD > 0) in the southeastern and northwestern areas, and this result contrasts with the implications of the Spearman correlation coefficient for CS-WP. This result suggests that a decrease in CS might decrease nitrogen export. This decrease may have occurred because cultivated land was converted to built-up land, and the increase in built-up land area might decrease CS. Additionally, built-up land produced less nitrogen loads than cultivated land. CS and HQ mainly exhibited a synergistic relationship throughout the YRD region; this relationship can be attributed to high consistency in the land cover types of CS and HQ. For example, ecological land normally has higher values of CS and HQ. However, built-up land has a lower capacity for CS and HQ. The high TSD values for CS-HQ were mainly distributed in the eastern part of the YRD region, and the values gradually decreased from the eastern area to the western area of the YRD region. The Shanghai and Suzhou-Wuxi-Changzhou regions in the eastern area of the YRD region are the most developed areas and exhibit a significant increase in built-up land to meet increasing demand for socioeconomic development. Therefore, the trade-off intensity in this area is higher than in other areas.
Similar to the spatial distribution of the TSD value for CS-WP, a large proportion of counties exhibited trade-off relationships between WP and HQ in the eastern and southern areas of the YRD region, which means that the increase in nitrogen export might decrease HQ. In addition, the trade-off intensity gradually decreased from the eastern areas to the southern areas. This decrease may have occurred because newly cultivated land occupies a large amount of ecological land, which leads to a decrease in HQ and an increase in WP. Positive TSD values are observed in the northwestern and southeastern areas, which indicates that the relationship between WP and HQ is synergistic. Cultivated land was converted to built-up land due to the rapid socioeconomic development. Cultivated land performed worse in water purification because of its high proportion of agricultural The global Moran's I for the TSD of each paired ES is presented in Table 4. The values of 0.8309, 0.3979, and 0.7216 for CS-WP, CS-HQ, and WP-HQ, respectively, suggest that the patterns of the trade-off/synergy relationships between the paired ESs were spatially autocorrelated.

The Impacts of Related Factors on the Trade-Offs/Synergies between ESs
We applied the MLR and GWR models to investigate the spatial relationships between the TSD and six natural and socioeconomic factors (ELE, PRE, TEM, GDP, POP, and BUL). In the study, VIF (variance inflation factor) was used to verify whether there is multicollinearity in the explanatory variables. If VIF exceeds 7.5, it means that the variable may be redundant. The results show that VIF values are all less than 7.5, which implies that only slight collinearity exists in the explanatory variables. The range of standardized residual values of the GWR model is −5.07-3.26, of which approximately 98.1% is −2.58-2.58. Therefore, the standardized residual values of the GWR model are randomly distributed with 95% confidence. As shown in Table 5, the R 2 values of the GWR results are larger than those of the MLR model. The AIC values of the GWR are smaller than those of the MLR model. In contrast to the results of the MLR model, which is a global model for estimating the relationship between two variables, the GWR results provide a specific coefficient for each sample. The effects of the spatial relationship between the explanatory and dependent variables on the trade-offs/synergies between the ESs is explored.   Previous studies have highlighted elevation as one of the significant influencing factors affecting ecosystem services [39,51]. Elevation was positively correlated with the tradeoff/synergy relationships between CS and WP in the eastern part of the YRD region. In contrast, a negative effect of elevation was mainly observed in the western area. It is found that the areas where a positive coefficient between elevation and TSD was present experienced a significant urban expansion, for example Shanghai, Suzhou, Nanjing. For a large area of the YRD region, there was a negative association between elevation and the CS-HQ relationship.
Vegetation is often influenced by climatic factors that affect their growth [52], which means that ESs are also influenced by the climatic factors. The direction and extent of the response of ESs to climatic factors normally present regional characteristics [53]. Precipitation is an important variable of the trade-off/synergy relationships from 2005-2015 ( Figure 8). Changes in precipitation had a negative effect on CS-WP and WP-HQ. This result implies that an increase in precipitation may either increase the trade-off intensity or reduce the synergy intensity. The negative correlation between precipitation and the trade-off/synergy intensity for CS-HQ in the southern area was related to the relatively slow changes in CS and HQ.
The changes in TEM also had significant effects on the trade-off/synergy relationships between the paired ESs. Both positive and negative correlations between TEM and TSD were found across the YRD region. As presented in Figure 8, the spatial distributions of the coefficients for CS-WP and WP-HQ are similar. The change in TEM in the eastern and western areas exhibited a negative correlation with the trade-off/synergy intensity of CS-WP and WP-HQ. TEM had a positive effect on the TSD in the central and northern areas of the YRD region. As evidenced by Schuur (2003), the effects of climate factors on ESs are spatially heterogeneous [53].
As illustrated in Figure 8, GDP and POP have positive and negative correlations with the trade-offs/synergies between ESs. TSD values for CS-WP, CS-HQ, and WP-HQ have significant positive relationships with GDP and POP in the western area of the YRD region, where the intensity of human activity and the socioeconomic level are relatively weak compared with those in the eastern region. In the eastern area of the YRD region, however, an increase in GDP or POP will increase the probability that the relationships will convert to a trade-off.
The changes in built-up land are negatively correlated with CS-WP and WP-HQ relationships in the eastern area of the YRD region, where built-up land is concentrated and increased significantly over the study period ( Figure 8). In addition, BUL exhibited a positive correlation with the synergy intensity of CS-HQ across the YRD region, and its effects gradually increased from southern areas to northern areas. This result implies that the synergy between CS and HQ is more sensitive to changes in BUL. CS and HQ are affected by a higher level of land-use consistency, and the increase in BUL could decrease the area of land cover types that support CS and HQ.

Spatially Varying Trade-Offs/Synergies among the ESs and the Underlying Factors
By understanding the spatial patterns of the ES relationships' responses to influencing factors, decision makers can regulate these factors to improve the synergy intensity and reduce unwanted trade-offs [44]. Our findings show that the trade-offs and synergies among multiple ESs (CS, WP, and HQ) were widespread in the YRD region, which is in accordance with the findings of previous studies [54,55]. However, studies concerning the spatial distribution of the intensity of trade-offs/synergies are scarce. Pairwise correlation coefficients are commonly used to reveal the global relationships between paired ESs [31,56]. This type of analysis reflects only the magnitude of a relationship and does not provide insight into the spatial patterns of a relationship. The hot spots of various ESs were compared to investigate the spatial relationship between two ESs. However, the results generated by hot spot comparisons do not indicate where trade-offs and synergies exist.
The RMSE indicator was widely adopted in previous studies to measure the intensity of trade-offs for one temporal point [10,25], but changes in ESs over time are not considered in this indicator. Zhang et al. (2020) used a binary value (0, 1) to measure the relationship by comparing the changes in ESs between two points in time [29]. Although this type of analysis can identify the ES relationships, the magnitude of the relationships is not measured. Therefore, it is difficult to fully reveal the impact mechanisms of the tradeoffs/synergies. In this study, we proposed a new indicator (TSD) to quantify the intensity of the trade-off/synergy relationships among ESs and their spatial heterogeneity rather than their spatial homogeneity. Compared with the findings of previous studies, our findings clarify the locations where trade-offs/synergies existed and their intensity.
Our study supports the argument that the relationships between ESs can be affected by natural and socioeconomic factors [57,58]. For appropriate ES management, decision makers should consider not only global correlations but also the spatial heterogeneity of the impacts to mitigate unwanted ES trade-offs. In contrast to the findings obtained by using correlation analysis and a global regression model [14,33], our results obtained by the GWR model can reveal the spatial relationship between trade-offs/synergies and the related factors. The intensity of trade-off/synergy relationships was significantly affected by natural factors and socioeconomic factors. A larger number of studies have confirmed that ESs are correlated with urbanization [59]. As indicated by the coefficient of determination (R 2 ) from the GWR model, BUL had the largest effects on the relationships between paired ESs among these factors in the study. Su et al. (2014) argued that the effects of urbanization on Ess are spatially heterogeneous, with the magnitude of the effects varying spatially, which is similar to our findings [37]. However, previous studies considering topics related to the effects of influencing factors on the complex relationships that exist among multiple Ess were limited by the use of a global regression model and correlation analysis without considering spatial non-stationarity. In this study, goodness of fit tests employed to compare the MLR and GWR models were assessed using R 2 and AIC.
The results indicate that the GWR model has a strong explanatory capability for estimating the relationships between the influencing factors and the trade-off/synergy intensity of ESs. The coefficients estimated by the GWR model correspond to a specific geographic area, and the spatially varying effects of various factors can be explored through the creation of a spatial distribution map of the GWR coefficients.
A decrease in elevation could lead to an increase in synergy intensity of CS-HQ. This result may be related to the distribution of ecological land across different elevations. People prefer to live in lower terrain gradients, leading to a low ecosystem service values, and vice versa. When there is a decrease in elevation, the proportion of forest decreases and that of cultivated land and built-up land increases. The urban expansion rate in this area is relatively higher than that in areas in higher elevations. The decreases in CS and HQ are relatively significant; thus, TSD value was enhanced. In the YRD region, the vegetation coverage area increased as precipitation increased, CS and HQ were enhanced, and WP decreased; thus, trade-off intensities of CS-WP and WP-HQ increased. Additionally, the effects of precipitation in the southern area of the YRD region are greater than those in other areas due to its higher proportion of ecological land and significantly higher levels of precipitation. With increases in GDP and population, socio-economic dependent land continuously occupied ecological land, and waste gas and water put added pressure on regulating and supporting services of the ecosystem. The newly expanded built-up land occupied a larger amount of cultivated and ecological land, which led to a significant decline in CS and HQ because of the loss of land cover classes that support CS and HQ. In contrast, urbanization was positively correlated with nitrogen export due to heavy pollutant loads on impervious surfaces and the weak purification ability of the limited amount of vegetation coverage.

Policy Implications
Integrating information about the trade-offs/synergies between ESs into ecological planning and management is a prerequisite for coordinating the relationship between ecological protection and socioeconomic development [60]. This study quantitatively explores the influencing mechanisms of natural and socioeconomic processes on the tradeoff/synergy relationships among ESs, which is helpful for developing scientific and effective ecosystem service management practices and regional planning. Recommendations for ES regulation and promotion that correspond to our results are as follows: First, in the context of the collaborative development of the YRD region, it is crucial to enhance the implementation of precisely targeted ecological projects and refine ecological management policies based on the analysis of the trade-offs/synergies among ESs. In recent years, a range of ecological projects have been carried out, such as the Grain for Green Project, the Natural Forest Protection Project, and Reclaiming Lake from Farmland [61]. These ecological projects play a key role in enhancing the ecological functions of regional water conservation and soil and water conservation and accelerating the construction of a green ecological corridor in the YRD region. However, ES degradation was still observed in some ecologically sensitive and rapidly developed areas. The YRD region still faces a great challenge in achieving the 'win-win' goal of socioeconomic development and ecological protection. Reaching this goal may be difficult because the trade-off/synergy relationships and their spatial heterogeneity, and impact factors have been neglected in ecological projects. Therefore, it is important to integrate the complex interactions among ESs into the design and implementation of various ecological projects to promote the efficiency of ecological protection and restoration.
Second, our findings reveal the specific regions where trade-off/synergy relationships exist, which is helpful information for decision making regarding trade-offs in specific locations and implementing ES management more effectively. The YRD region can be separated into three areas based on the intensity of land development, population concentration and changes in ESs: optimized development area, key development area and restricted development area. Different development strategies should be developed for various areas. The optimized development area refers to the area where the carrying capacity of resources and the environment is saturated. This area is mainly distributed in Shanghai, Hangzhou, and southern Jiangsu. We should take effective measures to strictly control the scale and intensity of new construction land and appropriately expand agricultural and ecological spaces. The key development area refers to the area with great potential for the carrying capacity of resources and the environment. This area is mainly distributed in central Jiangsu, central Zhejiang, central Anhui and the coastal area. It is necessary to strengthen the capacity of industrial and population agglomeration, appropriately expand industrial and urban spaces, optimize the rural living space, and strictly protect the green ecological space. The restricted development area refers to the area with significant ecological sensitivity and a low carrying capacity of resources and the environment. This area is mainly distributed in northern Jiangsu, western Anhui and western Zhejiang. It is necessary to strictly control the scale of new construction land; promote the concentrated development of cities; strengthen the protection of water resources, ecological restoration and construction; and maintain the stability of the ecosystem structure and its functioning.
Third, the land development intensity in the YRD region was 17.1% in 2013, which was 15% higher than that due to Japan's Pacific coastal urban agglomeration, and the potential space available for subsequent construction is insufficient. Shanghai's development intensity is as high as 36%, far more than that of Paris at 21% and that of London at 24%. Due to extensive and unrestrained development, new development zones and industrial parks occupy too much land, which seriously affects the overall ecosystem structure and utilization efficiency of regional land space [62]. Therefore, the "three red lines" (Ecological Protection Red Line, Cultivated land Protection Red Line, and Urban Development Boundary Red Line) policy for spatial planning was proposed to improve land-use efficiency [63]. Despite the common acknowledgements of the importance of the three red lines in sustainable

Limitations and Further Study
Although this study obtained some interesting findings, some limitations exist. Firstly, the TSD indicator was developed and applied to analyze the trade-offs/synergies among ESs in the YRD region and effectively quantify the complex relationships that exist between ESs. However, a temporal scale was not included in this study. Some scholars have pointed out that the relationship between ESs varies temporally [19,21]. In addition, the reliability of the results regarding the trade-offs/synergies can be improved by using long-term time series data. Therefore, future studies need to identify the relationships between ESs using long-term series data to mitigate uncertainty about the trade-offs/synergies that exist among ESs and analyze the impacts at the temporal scale. Secondly, the TSD indicator did not consider the spatial dependence of ESs. As some studies suggested that the supply of ESs in any given location depends upon juxtaposition with other ESs and on the overall mosaic of ESs across a region [64,65]. It is therefore, crucial to consider spatial dependence of ESs in future studies. Thirdly, limitations exist in terms of the analysis of the influence mechanisms of the trade-offs/synergies among ESs. This study investigated the effects of natural and socioeconomic factors on the trade-off/synergy relationships among ESs by considering six factors (ELE, PRE, TEM, GDP, POP, BUL). The findings suggest that these factors have significant effects on ES relationships. However, the complex relationships between ESs may be affected by other factors (e.g., soil properties, slope, solar radiation) [66,67]. It is of great importance to consider more potential factors to improve the understanding of the influence mechanisms of the trade-offs/synergies that exist among ESs. Finally, it is widely acknowledged that a better understanding of the relationships between ES supply and demand is important for sustainable ES management and the improvement of human wellbeing. However, the study analyzed the trade-off/synergy relationships from the perspective of ES supply. Therefore, the supply-demand relationships should be involved in future studies. The results of such an analysis will provide a more scientific support for sustainable development.

Conclusions
The results of the trade-off/synergy analysis enable decision makers to develop effective regional planning and ecological management policies that promote sustainable development. In this study, we investigated spatiotemporal changes in three ESs, namely, CS, WP, and HQ, in the YRD region from 2005 to 2015. The results revealed that CS, WP, and HQ decreased by 0.28%, 2.49%, and 3.38%, respectively. These variations in three key ESs in the YRD region imply that ESs have degraded mainly due to intensive human disturbance in recent years. In addition to measuring the trade-off/synergy relationships between multiple ESs, this study measured their magnitude using an innovative indicator (TSD). The TSD indicator was proven to be efficient for identifying the types of relationships that exist between ESs and quantifying the intensity of the trade-offs/synergies. The relationships between the factors and trade-offs/synergies were revealed by using MLR and GWR models. This study demonstrated that the trade-off/synergy relationships presented significant autocorrelations. The GWR model performed better than the MLR model in explaining variations in the trade-off/synergy intensity, as the GWR model generated higher R 2 values and lower AIC values. The impacts of natural and socioeconomic factors on the trade-offs/synergies were spatially heterogeneous rather than spatially consistent. Our findings show that the GWR model can be used to obtain precise information on the various roles of the related factors at different sites in the study area rather than producing a global coefficient for the entire area. The trade-off/synergy intensity is significantly correlated with meteorological, urbanization, and terrain factors. The findings improve the understanding of the trade-off/synergy relationships that exist among ESs and their influencing mechanisms. The trade-off/synergy relationships among ESs in the YRD region are not only affected by environmental factors but also significantly related to socioeconomic development. Furthermore, the spatial heterogeneity of the trade-offs/synergies is well explained by the influencing factors. The GWR model revealed that these relationships are spatially heterogeneous rather than spatially consistent relationships.
Supplementary Materials: The following are available online at https://www.mdpi.com/article/ 10.3390/land11010106/s1, Figure S1: Coefficients of the Spearman correlations between paired ESs (** p < 0.01), Table S1: Carbon density per unit area of different land cover types in the YRD region, China (Unit: Mg/ha), Table S2: Biophysical table in water purification module for the YRD region, China, Table S3: Parameters for habitat quality in the YRD region, China.