Analysis of interrelated characteristics between ecosystem services and ecosystem health in the Guangdong–Hong Kong–Macao Greater Bay Area

Ecosystem health (EH) underpins the capacity of vegetation ecosystems to provide essential ecosystem services (ESs), which together are fundamental to regional sustainability. In regions undergoing rapid urbanization, the interrelationships between EH and ESs become increasingly complex, yet they remain largely unexplored in previous studies. This study integrates the VOR and InVEST models to quantify EH and four key ESs in the Guangdong–Hong Kong–Macao Greater Bay Area (GBA) from 2000 to 2020 and further analyzes their interrelationships using a bivariate spatial autocorrelation model and the XGBoost-SHAP approach. The results indicate that: (1) From 2000 to 2020, low-value areas of most ESs and EH expanded, regions of EH deterioration accounted for 71.75% of the study area, indicating the profound impact of rapid urbanization. (2) EH showed strong positive global spatial correlations with CS and NPP, but weak negative spatial correlations with FP and WY. (3) Interrelationships between ESs and EH can be divided into stable synergy type and dynamic trade-off type based on their differing ecological processes; climate factors can significantly impact the interrelationships primarily by affecting the dynamic trade-off type. This study integrates spatial analysis and machine learning approaches to examine the relationships between EH and ESs, thereby advancing the understanding of ecosystem states and functions and providing a theoretical basis for formulating ecological restoration targets.

1 Introduction

Ecosystems sustain human societies by delivering a wide range of benefits, which are commonly understood through the concept of ecosystem services (ES). Ecosystem services refer to the ecological attributes, processes, and functions that underpin human well-being, highlighting the ways in which functioning ecosystems contribute to people’s lives (Costanza et al., 1997). Urbanization profoundly alters natural ecosystems (Grimm et al., 2008), causing declines in ESs and significantly affecting the well-being of humans (Millennium Ecosystem Assessment, 2005). In this context, enhancing ESs to sustain the well-being of humans has become a significant issue guiding environmental conservation. However, ecosystems are complex systems influenced by multiple factors. To ensure the sustainable provision of ESs, ecosystem management must focus not only on ESs but also on their mechanisms of generation (Bennett et al., 2009), which corresponds to the condition of ecosystems. Therefore, the ultimate goal of ecosystem management is to ensure the provision of diverse ESs while maintaining optimal ecosystem health (EH) (Rapport et al., 1999). EH refers to the structural and functional integrity of an ecosystem, encompassing its ability to remain active, organized, and self-sustaining, while also demonstrating resilience to stress and disturbances (Costanza, 1992; Rapport et al., 1998). In addition to the provision of ESs, EH is also crucial, as it represents the concept of an ideal ecosystem. Therefore, a thorough investigation of the EH-ESs relationship in a specific region is crucial, which will serve as a scientific foundation for effectively promoting their integrated enhancement and coordinated improvement.

Regarding the interrelationship between ESs and EH, widespread consensus indicates that healthy ecosystems can provide ESs sustainably and steadily (Costanza, 2012; Hernández-Blanco et al., 2022). However, quantitative research reveals that relationships between the two may not always reflect strictly synergistic characteristics (Zhu and Cai, 2024). In regions undergoing rapid urbanization, where dramatic changes are occurring in the structure, function, and state of ecosystems, the complexity of interrelationship between EH and ESs may increase. Therefore, it is imperative to test the hypothesis that a healthy ecosystem produces high levels of ESs, especially in rapidly urbanizing areas. However, research on the integration of ESs and EH remains limited. Most related studies tend to examine EH and ESs separately. In studies on EH, scholars have primarily employed evaluation models such as PSR and VOR to quantitatively assess ecosystem health (Das et al., 2021; Li et al., 2024). Subsequent analyses often include investigations of driving factors and corresponding ecological zoning (Han et al., 2024; Li et al., 2025; Shao et al., 2025). Research on ecosystem services, which has become a central focus in contemporary ecology, mainly encompasses trade-off and synergy analyses based on the quantification of ecosystem services (Feng et al., 2021; Xue et al., 2023), cluster analysis (Dou et al., 2020; Xia et al., 2023; Su et al., 2024), and mechanism analysis of influencing factors (Liu et al., 2019; Jianying et al., 2020; Qi et al., 2020).

Existing studies that combine EH and ESs mainly link ESs to the assessment of EH (Peng et al., 2015; Pan et al., 2021; Yu et al., 2022; Huang et al., 2024), grounded in the theory that healthy ecosystems can provide abundant service functions. Some studies also take them as indicators to quantitatively measure ecological indices. For example, some scholars analyzed regional ecological risk by taking the square root of the product of ecosystem health and service function indicators in the Fen River basin (Wang et al., 2023a) and the Beijing-Tianjin-Hebei urban agglomeration (Kang et al., 2018). Other scholars have integrated EH and ESs with additional indicators and conducted comprehensive analysis of these indicators to measure ecological security in the Huaihe River Basin (Zhu and Cai, 2023). However, the above studies focus primarily on the integrated application of these two indicators rather than investigating the underlying relationship mechanisms between them. Given the fundamental differences between EH and ESs in conceptual frameworks and quantification methods, they are likely to exhibit distinct spatiotemporal patterns and response mechanisms. Therefore, it is necessary to understand their intrinsic relationships rather than simply integrating these indicators. In existing research on the interrelationship between EH and ESs, Liu et al. (2015), using correlation coefficient calculation methods, found a generally synergistic relationship in the forest areas of Northeast China. Similarly, Zhu et al. (2021), also employing correlation coefficient calculation methods, discovered that in the Pingjiang watershed, higher regional ecological quality is associated with lower food supply services, with the trade-off between NDVI and food supply being the most significant. In reality, ESs and EH are the results of interactions among numerous ecological variables, the relationship between them may exhibit complex non-linear characteristics. Correlation coefficients overlook the uncertainties introduced by these characteristics and fails to reflect the complex relationship. Therefore, further systematic analysis and research are needed to reveal the complex characteristics and mechanisms of the interrelationship between ESs and EH amid rapid urbanization.

As a prominent economic hub in China, the Guangdong–Hong Kong–Macao Greater Bay Area (GBA) is one of the four leading bay areas worldwide. By 2023, the total population of GBA had exceeded 86 million, and its economic output had surpassed 14 trillion yuan. Taking up a fraction of one percent of China’s land area, the GBA contributes to 1/9th of the national economic total, playing a pivotal role in the quality-driven growth of China (Wang et al., 2023b). The GBA boasts favorable natural geographic conditions, but rapid urbanization has stressed the natural ecosystems. Previous research has shown that the GBA is facing increasing environmental and resource pressure (Wu et al., 2021), and urban carbon emissions grew rapidly by 7% from 2000 to 2011 (Zhou et al., 2018). Therefore, the GBA is a typical example of changes in EH and ESs driven by rapid urbanization. Guided by the Green and Beautiful Guangdong Campaign, the GBA plays an essential role as a model for the development of bay areas in China and is positioned to be a pioneer for the Beautiful China Initiative, which makes ecological construction in the GBA extremely important. Consequently, enriching research on ecosystems in the GBA is important for providing the theoretical foundation for more scientifically informed ecological management practices. Therefore, this study integrates the InVEST model, VOR model, bivariate Moran’s I method and machine learning techniques to explore the interrelated characteristics between EH and ESs. This research aims to: (1) quantify EH and ESs, exploring their spatiotemporal distribution patterns; (2) examine the spatial interrelationships between EH and ESs; and (3) investigate the variable interrelated characteristics between EH and ESs. The findings will provide deeper insights into the interactions between EH and ESs, which will offer a valuable foundation for designing ecological restoration strategies in regions sharing comparable geographical characteristics.

2 Materials and methods

2.1 Study sites

The GBA spans an area of 5.6×10⁴ km², covering coordinates from 21°25’ to 24°30’N and 111°12’ to 115°35’E (Figure 1A). It is characterized by an overlap of rainfall and temperature seasonality under a subtropical monsoon climate. It experiences warm and humid conditions year-round, which create favorable conditions for both heat and water. The mountains are concentrated in the northern part, while plains are predominantly located in the central region and along the coast. The dominant vegetation type is subtropical evergreen broadleaf forest. As a key economic hub of China, the GBA has undergone rapid urbanization in the past few decades, significantly transforming its surface landscape pattern (Figures 1B-D).

Figure 1

Figure 1. Geographical location (a) and land use types (b-d) of the Guangdong–Hong Kong–Macao Greater Bay Area.

2.2 Data sources and methodology

Multivariate data were collected to quantify EH and ESs, and the data used are summarized in Table 1, where the first column lists the types of data employed, and the second and third columns provide the source of each dataset along with additional details such as spatial resolution and temporal coverage (Table 1).

Table 1

Table 1. Data and their sources.

To enhance the accuracy of spatial interpolation, we incorporated meteorological station data from both the GBA and its adjacent regions into calculation. In ArcGIS Pro 3.0.2 software, we obtained meteorological raster data by applying the Kriging spatial interpolation technique. We standardized all spatial datasets to the WGS 1984 UTM Zone 49N coordinate system based on the study area’s location.

2.3 Quantification of ecosystem health

We quantified EH using the VOR framework, which is considered to characterize ecosystem structure and function completely (Costanza, 2012; Das et al., 2021).

The VOR framework includes the assessment of ecosystem vitality, organization and resilience. We measured ecosystem vitality using the Normalized Difference Vegetation Index (NDVI), which is regarded as an effective way to evaluate the vitality (Box et al., 1989; Peng et al., 2015); Using the moving window method by Fragstats 4.2 software, we quantified ecosystem organization thorough three landscape metrics: Shannon’s Diversity Index (SHDI), Contagion Index (CONTAG), and Patch Cohesion Index (COHESION) (Frondoni et al., 2011; Peng et al., 2017). Using InVEST v3.14.1 software, we assessed ecosystem resilience through the Habitat Quality sub-module (Ye et al., 2013; Sharp et al., 2018; Chen et al., 2022). We used Equations 1–4 for the calculation:

$\begin{array}{ll}{i}_{rxy}=1-\left(\frac{d_{xy}}{d_{rmax}}\right),\text{if linear}& (1)\end{array}$

$\begin{array}{ll}{i}_{rxy}=\text{exp} (-(\frac{2.99}{d_{rmax}}\left)d_{xy}\right),\text{if exponential}& (2)\end{array}$

$\begin{array}{ll}{D}_{xj}={\displaystyle {\sum }_{r=1}^R{\displaystyle {\sum }_{y=1}^{Y_r}(\frac{W_r}{{\sum }_{r=1}^R{W}_r}}})r_y{i}_{rxy}{\beta }_x{S}_{jr}& (3)\end{array}$

$\begin{array}{ll}{Q}_{xj}=H_j(1-\frac{D_{xj}^2}{D_{xj}^2+{\text{k}}^2})& (4)\end{array}$

where $d_{xy}$ represents the linear distance between raster cells $x$ and $y$; $d_{rmax}$ denotes the maximum effective distance of threat $r$’s reach across space; $i_{rxy}$ is the impact of the threat factor $\text{r}$ originating from the cell $x$ on the cell $\text{y}$; $D_{xj}$ represents the habitat degradation risk index of the cell $x$ in the habitat type $j$; $y$ is the cell under the threat factor $r$; $r_y$ is used to determine if the cell $y$ provides a source of threat factor $r$; $W_r$ is the threat weight of the factor $r$; ${\beta }_x$ denotes the accessibility of the threat factor of the cell $x$; $S_{jr}$ represents the sensitivity coefficient of the land use type $j$ to the factor $r$; $Q_{xj}$ is the quality of habitat of the raster $\text{x}$ in habitat type $j$; and $k$ is the half-saturation constant, often be set to 0.05.

The TOPSIS method is a widely used evaluation technique that identifies the best choice by calculating the distances between the data set and both the ideal and worst solutions (Hwang and Yoon, 1981). To reduce the subjectivity of weighting decisions, we combined the Entropy Weight method with the TOPSIS method for the quantification of EH (Lv et al., 2023; Wu et al., 2021).

2.4 Quantification of ecosystem services

We selected four ecosystem service indicators as representatives based on the context of the GBA, including carbon sequestration (CS), net primary productivity (NPP), food production (FP), and water yield (WY). (1) The dense population and advanced economic development in the GBA make the ecosystem’s provision of material and energy services crucial, and NPP reflects the material production capacity of the ecosystem (Nemani et al., 2003). (2) With the reduction in per capita farmland area over the past few decades (Zhang et al., 2021), research on FP is crucial for ensuring food security and understanding the effects of urbanization on FP. (3) As economic and population growth increases demand for clean water and hydroelectric resources, assessing WY helps optimize water resource management for sustainable development (Wang et al., 2021). (4) Carbon emissions are high in the GBA, while the superior vegetation conditions provide a basis for CS. Studying CS is important for achieving carbon neutrality (Zhou et al., 2018).

2.4.1 Carbon sequestration

CS in a given area is determined by the size of four carbon pool: above-ground biomass, below-ground biomass, soil carbon, and dead organic matter. We used the Carbon Storage and Sequestration sub-module in InVEST v3.14.1 (Delphin et al., 2013; Sharp et al., 2018; Chu et al., 2019) to assess CS. We calculated the total carbon storage for each land use type using Equation 5:

$\begin{array}{ll}{C}_{totali}=(C_{abovei}+C_{belowi}+C_{soili}+C_{deadi})\times A_i& (5)\end{array}$

where $A_i$ represents the area of land use type $i$; $C_{abovei}$ is the above-ground carbon storage per unit area of land use type $i$; $C_{belowi}$ is the below-ground carbon storage per unit area of land use type $i$; $C_{soili}$ is soil carbon storage per unit area of land use type $i$; $C_{deadi}$ is dead organic carbon storage per unit area of land use type $i$; and $C_{totali}$ represents the total carbon storage of land use type $i$.

2.4.2 Net primary productivity

NPP is the total amount of organic matter accumulated by vegetation per unit area and per unit time, reflecting the material supply capacity of the ecosystem (Nemani et al., 2003; Imhoff et al., 2004). We downloaded the MOD17A3HGF006 vegetation NPP dataset from the Analytical Insight of Earth Platform to assess NPP.

2.4.3 Water yield

The Annual Water Yield sub-module in InVEST v3.14.1 is an estimation method using the water balance method (Redhead et al., 2016; Sharp et al., 2018), calculated WY by Equation 6:

$\begin{array}{ll}Y(x)=\left(1-\frac{AET(x)}{P(x)}\right)·P(x)& (6)\end{array}$

where $Y(x)$ represents the annual water yield (mm) of the raster cell $x$; $AET(x)$ is the annual actual evapotranspiration (mm) of the cell $x$; $P(x)$ is the annual precipitation (mm) of the cell $x$. We calculated evapotranspiration using the Penman-Monteith equation. The maximum root burial depth data required by the model was substituted with the absolute depth to bedrock data, while the plant available water content was calculated from available soil water capacity until wilting point data, which is recommended in the user’s guide in InVEST software (Sharp et al., 2018). Zhang’s coefficient is a constant that represents the seasonal characteristics of precipitation, we adjusted it for different years to ensure the model simulation results align with government water resource reports.

2.4.4 Food production

Due to the significant linear correlation between grain production and NDVI (Chen et al., 2020), we combined the NDVI values with grain production statistics collected from the Statistical Yearbook to better quantify the value of food production (Huang et al., 2023). We computed FP using Equation 7:

$\begin{array}{ll}F{P}_i=G_{\text{sum }}\times \frac{{\overline{NDVI}}_i}{{\overline{NDVI}}_{\text{sum }}}& (7)\end{array}$

where $F{P}_i$ represents the food production value of raster cell $i$; $G_{\text{sum }}$ represents the total regional production of crops; ${\overline{NDVI}}_i$ denotes the annual mean NDVI value of the cell whose land use type is farmland; and ${\overline{NDVI}}_{\text{sum }}$ represents the sum of NDVI values across all raster in the GBA. To achieve more realistic results, we performed calculations separately for each city based on the respective statistical data.

2.5 Analysis of the interrelated characteristics between EH and ESs

ESs and EH are two comprehensive indicators of ecosystems, and the relationship between the two is likely to be complex, as they result from the interaction of many environmental and ecological variables. Due to spatial interactions and diffusion, spatial autocorrelations can cause biases in the relationship between these indicators when using methods like ordinary least squares or geographically weighted regression model (Zhang et al., 2018). In contrast, the bivariate Moran’s I method effectively reflects the interrelationships of two variables spatially (Anselin, 1995; Xu et al., 2019), making it suitable for revealing bivariate spatial correlations. Thus, we use this method to reveal spatial interrelated characteristics between the indicators.

Besides research on spatial correlation features, combining methods that intuitively reveal the interrelated characteristics between variables may further enhance the understanding of their relationship. In recent years, machine learning models have gained popularity in geographic and ecological analysis due to their ability to effectively process multidimensional data and reveal nonlinear relationships between variables without predefined norms. The Extreme Gradient Boosting (XGBoost) model, which adds regularization rules to reduce the risk of overfitting, has significantly improved algorithm efficiency and accuracy, making it highly effective in ecological studies (Yang et al., 2021; Yuan et al., 2024). Additionally, SHAP (Shapley Additive Explanations), based on game theory and local interpretability, is a classic method used to explain the results of machine learning. The partial dependence plots generated by SHAP help in understanding the complex nonlinear interrelationships between explanatory variables and dependent variables (Lundberg et al., 2018). Based on this knowledge, we employed the XGBoost-SHAP model to further analyze the relationship between ESs and EH.

2.5.1 Spatial interrelated characteristics

We examined the spatial interrelationships between EH and ESs via the GeoDa 1.18.0 software, employing both global and local spatial autocorrelation methods. We applied global Moran’s I to identify the overall spatial correlation between EH and ESs, with values ranging from -1 to 1. We used Equation 8 to calculate the global Moran's I.

$\begin{array}{ll}{I}_{global}^{ab}=\frac{n}{S_0}\frac{{\sum }_{i=1}^n{\sum }_{j=1}^n{w}_{ij}(X_{ai}-{\overline{X}}_a)(X_{bj}-{\overline{X}}_b)}{{\sum }_{i=1}^n{(X_{ai}-{\overline{X}}_a)}^2}& (8)\end{array}$

where $n$ is the total number of spatial units; $X_{ai}$ and $X_{bj}$ are the observed values of variables $a$ and $b$ for spatial units $i$ and $j$, respectively; ${\overline{X}}_a$​ and ${\overline{X}}_b$​ are their mean values; $w_{ij}$​ is the spatial weight describing the spatial relation between units $i$ and $j$; $S_0$ is the sum of spatial weights. $I_{global}^{ab}$ is the global bivariate Moran’s I index. A positive value suggests a positive spatial correlation, while a negative value suggests a negative spatial correlation. When $I_{global}^{ab}$ = 0, no spatial autocorrelation exists between the variables.

At the local level, we used the local Moran’s I (LISA) method to identify spatial agglomeration patterns between EH and ESs. We calculated local Moran's I using Equation 9:

$\begin{array}{ll}{I}_i^{ab}=z_{a,i}{\displaystyle {\sum }_{j=1}^n{w}_{ij}{z}_{b,j}}& (9)\end{array}$

where $z_{a,i}=(X_{ai}-{\overline{X}}_a)/{\sigma }_a$ and $z_{b,j}=(X_{bj}-{\overline{X}}_b)/{\sigma }_b$ are the standardized values of variables $a$ and $b$; ${\sigma }_a$ and ${\sigma }_b$ are their standard deviations; other symbols are as defined above. $I_i^{ab}$ is the local bivariate Moran’s I index. Based on the results of the bivariate Moran’s I analysis, we categorized the study area into five types: high-high type, indicating regions with high ESs and EH; low-low type, where both ESs and EH are low; high-low and low-high type, representing spatial heterogeneity; and non-significant type, indicating no obvious spatial agglomeration. In this study, we assessed the statistical significance of the results using a Z-test (P< 0.05).

2.5.2 Variable interrelated characteristics

We implemented XGBoost-SHAP analysis using Python’s scikit-learn, XGBoost, and SHAP packages for modeling and visualization. Using ArcGIS Pro 3.0.2 software, we generated 1.0×10⁵ random points across the study area to extract ESs and EH data. Before modeling, we standardized all features using the StandardScaler to ensure comparability.

For constructing the XGBoost regression model, we utilized continuous EH values as the dependent variable, with multiple ecosystem service indicators serving as independent variables. The XGBoost algorithm builds predictive models through iterative addition of decision trees, and its core objective function is formulated as shown in Equation 10:

$\begin{array}{ll}ℒ={\displaystyle {\sum }_{i=1}^{n}l(y_i,{\widehat{y}}_i)}+{\displaystyle {\sum }_{k=1}^K\text{Ω}(f_k)}& (10)\end{array}$

where $n$ represents the total training samples, $y_i$​ and ${\widehat{y}}_i$​ denote the observed and predicted values for sample $i$, respectively; $l(y_i,{\widehat{y}}_i)$ is the loss function; $K$ indicates the total number of trees, $f_k$​ represents the $k$ -th tree, and $\text{Ω}(f_k)$ is the regularization component that prevents overfitting.

We divided the dataset into training (80%) and validation (20%) sets. A grid search cross-validation (GridSearchCV) was performed to optimize the XGBoost hyperparameters, including maximum tree depth (3–6 nodes), learning rate (0.01-0.1), number of n estimators (100-300), min child weight (1, 3, 5 nodes) and subsample (0.8-1.0). We evaluated model performance by the coefficient of determination (R²). We employed Five-fold cross-validation to assess the model’s predictive accuracy and stability, ensuring robust performance across different data subsets.

For interpretability analysis, we employed SHAP methodology to quantify individual feature contributions to EH predictions. For a given observation $i$ and feature $j$, the SHAP value is mathematically defined as shown in Equation 11:

$\begin{array}{ll}{\varphi }_j(i)={\displaystyle {\sum }_{S\subseteq F\setminus \left\{j\right\}}\frac{\left|S\right|!(\left|F\right|-\left|S\right|-1)!}{\left|F\right|!}}[f(S\cup \left\{j\right\})-f(S)]& (11)\end{array}$

where $F$ represents the complete feature space, $S$ denotes any feature subset excluding $j$, and $[f(S\cup \left\{j\right\})-f(S)]$ quantifies the marginal contribution of feature $j$ when added to subset $S$. The SHAP value ${\varphi }_j(i)$ measures the contribution of feature $j$ to the prediction for observation $i$. Feature importance ranking was determined by computing mean absolute SHAP values across all observations, providing insights into the relative influence of different ecosystem services on environmental health outcomes.

3 Results

3.1 Spatial-temporal variations of ESs and EH in the GBA from 2000 to 2020

Both ESs and EH in the GBA revealed prominent spatial heterogeneity, with the distribution of various indicators exhibiting spatial autocorrelation and significant differences in their spatial patterns (Figure 2). CS and NPP showed similar spatial distribution patterns to EH (Figures 2A-I), generally decreasing from the periphery to the center and from north to south. The highest values were concentrated in the northern parts of Zhaoqing, Guangzhou, and Huizhou, whereas the areas with lower values were predominantly concentrated in the core area of the GBA. Apart from areas such as Kowloon, Yau Tsim Mong, and Kwun Tong, where the values of EH, CS and NPP were lower, Hong Kong generally demonstrated exceptional EH and ESs conditions, emerging as concentrated zone of high ecological value within the central regions of the GBA. The high-value areas of FP demonstrated a more balanced spatial arrangement throughout the landscape with lower spatial autocorrelation (Figures 2J-L). Regions with high FP mostly concentrated in suburban areas, while regions with high WY formed a northeast-southwest corridor through the central of the GBA (Figures 2M-O).

Figure 2

Figure 2. Quantification results and annual variations of EH (A-C) and ecosystem CS (D-F), NPP (G-I), FP (J-L), and WY (M-O) service functions from 2000 to 2020.

Over the 20-year period, the median EHI showed a gradual increasing trend. The median EHI value rose from 0.29 in 2000 to 0.32 in 2020, while the dispersion of the data gradually increased over the 20-year period (Figure 3A). Spatially, low-value zones expanded, with 71.75% of the area experiencing EH degradation. However, values in most original high-value zones increased, resulting in an overall polarizing trend. The overall change trends of CS and NPP values were not pronounced (Figures 3B, D), but their spatial distribution patterns evolved similarly to EHI, with both indicators showing increased diffusion of low-value areas over time while surrounding areas exhibited relatively minor changes. FP decreased most significantly, with the median value per raster cell declining from 2.19t to 0.50t (Figure 3C). This sharp decrease in FP occurred primarily in the core areas of the GBA. Meanwhile, the data dispersion of FP also decreased significantly. WY exhibited the greatest magnitude of fluctuation, varying irregularly over time, with the median WY value changing by up to 197.77 mm (Figure 3E).

Figure 3

Figure 3. Boxplots of the values of EH (A) and ecosystem CS (B), FP (C), NPP (D), and WY (E) service functions from 2000 to 2020.

3.2 Spatial interrelated characteristics between ESs and EH in the GBA

Overall, the positive spatial correlations between ESs and EH predominated (Figure 4), with distinct spatial variations across different ES-EH combinations (Figure 5).

Figure 4

Figure 4. The bivariate Moran’s I scatter plots of EH and CS (A, B), NPP (C, D), FP (E, F), and WY (G, H) from 2000 to 2020. The x-axis represents the standardized values of ESs, while the y-axis represents the spatially lagged EH values. Each graph includes a line of best fit and the corresponding Moran’s I value.

Figure 5

Figure 5. The bivariate LISA map of EH and CS (A, B), NPP (C, D), FP (E, F), and WY (G, H) from 2000 to 2020. The pie chart represents the proportion of different clustering types in relation to the total area of the GBA.

Global Moran’s I values exceeded 0.65 between EH and either CS or NPP, suggesting notable positive spatial correlation (P<0.05) (Figures 4A-D). The spatial clustering patterns between the two types of combination showed a high degree of overlap (Figures 5A-D), primarily characterized by the extensive distribution of high-high and low-low clusters. Over 30% of the area was occupied by high-high clusters, primarily located in the forested zones in the outskirts of the GBA. The low-low clusters covered over 20%, concentrated in the center of Guangzhou, Foshan, Shenzhen, and Dongguan.

Global Moran’s I values were negative for both EH-FP and EH-WY combinations (Figures 4E-H). The spatial clustering patterns for FP and EH were dominated by low-high clusters, which occupied more than 33% of the area and were located in the outskirts of the GBA (Figures 5E, F). The low-low clusters were predominantly located in the central regions of the GBA, surrounded by the low-high clusters. The spatial clustering pattern for WY and EH showed a more even distribution across types (Figures 4G, H), with the western part mainly exhibiting low-high and low-low clusters, and the eastern part showing high-low and high-high clusters. Temporally, the global Moran’s I between EH and either FP or WY decreased from 2000 to 2020. Low-low cluster areas expanded most notably in the EH-FP correlation pattern, with an increase of more than 8% in the GBA. The spatial interrelated characteristics between WY and EH varied dramatically with less discernible patterns.

3.3 Variable interrelated characteristics between ESs and EH in the GBA

The XGBoost model performance varied over time and among different ES and EH combinations (Table 2). The model of EH and CS achieved the highest accuracy with an average R² of 0.68, while the model of EH and NPP had an average R² of 0.53. The EH-WY model showed lower prediction accuracy with substantial temporal variability.

Table 2

Table 2. Model fitting accuracy values (R²).

Results of the XGBoost-SHAP model indicated that CS and EH showed a positive relationship close to linear (Figure 6A). The relationship between NPP and EH showed a clear threshold effect, with SHAP values stabilizing and even slightly declining when NPP reaches around 10,000 kg/m² (Figure 6B). The relationship between FP and EH initially exhibited a positive trend before shifting to a negative trend at lower values. When FP exceeded 10 tons, the EH prediction curve flattened. From 2000 to 2020, the curve’s peak point shifted rightward (Figure 6C). The relationship between WY and EH demonstrated a negative correlation because the prediction curves showed notable variation across years (Figure 6D).

Figure 6

Figure 6. Interrelated characteristics between EH and CS (A), NPP (B), FP (C), and WY (D) indicators from 2000 to 2020.

4 Discussion

Scholarly consensus suggests that healthy ecosystems can provide diverse ESs (Costanza, 2012). However, despite their conceptual relationship, the empirical evidence for specific relationships between EH and ESs within defined spatial-temporal contexts remains limited. To address this gap, our study leverages a combination of bivariate spatial autocorrelation modeling and the XGBoost-SHAP approach to quantitatively examine the interrelated characteristics between EH and ESs in the GBA, while the further understanding of the mechanisms behind these relationships and the comprehensive interaction between EH and ESs requires analysis of ecosystem processes.

4.1 Stable synergistic interrelated pattern between ESs and EH

EH showed similar spatial aggregation patterns with CS and NPP, characterized by significant positive correlations. Core areas with intense human activity experienced rapid urbanization over the past two decades, leading to declines in both EH and ESs and resulting in low-value and deteriorated zones (Gordon et al., 2009). In contrast, the outskirts of the GBA, less affected by urbanization, exhibited strong ecosystem vitality, organization, and resilience, supporting more stable material cycling and energy flow, and allowing rapid recovery from disturbances, thereby promoting CS and NPP (Chen et al., 2018; Reader et al., 2023). The gradual increase of the EH in surrounding areas can be attributed to ecological restoration projects in nearby mountains (Feng et al., 2021), plant growth promoted by climate change (Friend et al., 2013), or the self-recovery capacity of the ecosystem (Liu et al., 2023); all of these factors may enhance EH, CS, and NPP and reinforce their synergistic relationship.

At the process level, healthy vegetation ecosystems maintain high NPP and CS through photosynthesis, during which carbon is fixed and transformed into organic matter (Johnson, 2016; Walther et al., 2019; Smith et al., 2024). This establishes the coupling between CS and NPP and contributes to the stability of the ecosystem. Similar synergistic patterns have been observed elsewhere; for instance, Zhang et al. (2022) reported a correlation coefficient of 0.81 between habitat quality and carbon storage in the Chengdu–Chongqing Urban Agglomeration, and Huang et al. (2023) found a strong synergy between habitat quality and NPP in the Wujiang Basin. Urban land encroachment disrupts these ecological processes, reducing associated service functions, which has been widely observed in previous studies (Wu et al., 2024; Yushanjiang et al., 2024). Methodologically, previous studies often relied on correlation coefficients to assess overall relationships between ecosystem indicators, while more detailed analyses revealed spatial patterns of trade-offs and synergies (Xue et al., 2023). In this study, SHAP analysis further elucidated specific associations between EH and ESs. Notably, we identified a threshold effect in the NPP–EH relationship: areas with high NPP often corresponded to plantations around the GBA. Despite rapid material production, these plantations exhibit simplified ecosystem structures and are less capable of maintaining ecosystem health and multiple services (Dislich et al., 2017; Guillaume et al., 2018).

4.2 Dynamic trade-off interrelated pattern between ESs and EH

FP and WY exhibited more complex non-linear relationships with EH. EH increased initially with FP, but after a certain point, it decreased and then stabilized. FP is mostly enhanced by the use of pesticides, herbicides, and monoculture (Reader et al., 2023); these agricultural activities have resulted in simpler ecosystem structures and lower biodiversity, despite crop growth directly corresponding to increased ecosystem vitality on cropland. This dual effect explains the observed non-linear relationship between FP and EH. Temporally, from 2000 to 2020, the peak of the SHAP curve describing the relationship between FP and EH shifted to the right. The GBA used to be a major grain-producing region in Guangdong Province; however, rapid urbanization has directly consumed farmland, leading to decreased grain production (Wang et al., 2019). Consequently, regions with low EH and high FP shifted to low-low clusters, weakening the trade-off intensity between them, possibly contributing to the shift of the peak of the curve. Previous studies have indicated negative correlations between FP and other services; for example, Xia et al. (2023) revealed this pattern while also noting an expansion of synergy areas over time using a geographically weighted regression model. In contrast, our analysis found relatively few FP–EH trade-off zones, likely because the highly urbanized GBA contains extensive areas with simultaneously low EH and FP, weakening the observed trade-offs.

The interrelationship between EH and WY was not dominated by any single clustering type, and WY increased overall with noticeable fluctuations as EH decreased. The irregular clustering patterns and sharp fluctuations observed in the curves can largely be attributed to climatic variability, particularly precipitation. As an external driver, precipitation exhibits high uncertainty and weak correlation with surface conditions (Runting et al., 2017). Since WY is strongly controlled by precipitation, it exhibited greater variability than the relatively stable EH, leading to unstable interrelationship patterns. High–low clustering of WY and EH was mainly concentrated on the impervious surfaces at the urban core of the GBA, where rising WY due to blocked infiltration was accompanied by low EH. By contrast, forested areas showed higher EH but lower WY because of intense transpiration (Sun et al., 2006). Previous studies have reported higher WY in mountainous forests and lower WY in densely populated areas (Darvishi et al., 2022; Mo et al., 2023). The WY distribution in the GBA differs from these patterns, likely due to its monsoon-influenced coastal climate and extensive impervious surfaces typical of a large urban agglomeration. Regarding trade-offs between WY and other ecological indicators, prior studies have documented similar phenomena. For example, Liu et al. (2020) quantified runoff coefficients in Southwest China, revealing trade-offs between WY and carbon sequestration, while Zhang et al. (2022) found that in the Chengdu–Chongqing urban agglomeration, WY and habitat quality were largely dominated by trade-offs, with negative synergies exceeding positive ones. In contrast, our results show that in the GBA, high–high EH–WY clusters are more extensive than low–low clusters, likely reflecting abundant precipitation supporting higher WY. Methodologically, Zhang et al. focused on temporal trends, capturing WY declines during forest restoration, whereas our study emphasizes spatial patterns, revealing the coexistence of different ecological indicators across the GBA. Notably, SHAP analysis further identified the non-linear characteristics of the EH–WY relationship, highlighting its complexity under the combined influence of climate variability and human activities.

4.3 Hierarchical interrelated patterns between ESs and EH

The above analysis reveals that the interrelated patterns between ESs and EH vary depending on the type of ESs. This difference arises from the varying dependencies of ESs on ecosystem structures and processes.

CS and NPP are strongly correlated with ecosystem photosynthetic processes and serve as crucial indicators reflecting fundamental ecosystem metabolic processes. Their stable synergistic relationship with EH fully demonstrates the intrinsic resilience characteristics of ecosystems (Holling, 1973). In contrast, FP and WY are driven mainly by external factors such as human activities and climate change, reflecting ecosystems’ sensitive responses to external disturbances. Their dynamic trade-off relationship with EH profoundly reflects the adaptive adjustment processes of ecosystems within disturbance-recovery cycles (White and Pickett, 1985). Notably, the intensity of ecosystem responses to external disturbances tends to gradually diminish as ecosystem health status and fundamental supporting service functions improve. For instance, photosynthetic processes enhance the ecosystems’ buffering capacity against climate change by efficiently absorbing carbon dioxide and promoting long-term carbon storage. Under sustained moderate disturbance conditions, ecosystems can continuously evolve and enhance their stability (Holling, 1973; Walker et al., 2004). However, rapid urbanization has disrupted the inherent self-regulation thresholds of ecosystems (Muradian, 2001), making the differentiation mechanism of stability-dynamics interrelationship patterns more pronounced. This differentiation primarily stems from profound changes in land surface environmental characteristics induced by urbanization. These modifications not only amplify the fluctuation magnitude of key dynamic factors—such as temperature and precipitation—but also significantly intensify both the impact intensity and influence depth of these factors on corresponding ecosystem services.

The identification of hierarchical interrelated patterns represents an advancement beyond the traditional binary cognition of synergy-trade-off relationships between ecosystem attributes. By revealing the intrinsic hierarchical nature of ecosystem service-health relationships under rapid urbanization, this study provides empirical evidence. It offers novel insights into the classical theory that healthy ecosystems can sustainably provide abundant ESs. The stable synergistic characteristics confirm that the healthy ecosystems can maintain fundamental ecological processes. Meanwhile, the dynamic trade-off characteristics indicate that the supporting role of EH in supporting ESs is not a simple linear facilitation, but rather exhibits hierarchical differences based on the characteristics of ecological processes.

4.4 Limitations and future perspectives

This study investigates the interrelated characteristics between ESs and EH, advancing the understanding of their interactions. However, there remain certain limitations. Firstly, methodological constraints, including parameter dependencies in the InVEST model and subjective evaluations in the VOR model, may introduce quantification biases that potentially compromise regression fitting accuracy. Secondly, the types of ESs examined in this study are limited, necessitating future comprehensive investigations across a broader spectrum of ecosystem service types. Finally, while this study reveals the interrelated characteristics among indicators and provides reasonable interpretations under specific conditions, the underlying causal mechanisms require further investigation. This could be achieved through multi-scale and spatial-temporal quantitative analyses to provide more robust evidence for understanding the mechanistic interactions between ecosystem health and ecosystem services.

5 Conclusions

This study employs a comprehensive approach combining the InVEST model, VOR model, bivariate spatial autocorrelation analysis, and XGBoost-SHAP model to systematically investigate the interrelationships between ESs and EH in the Guangdong–Hong Kong–Macao Greater Bay Area. The findings reveal three key insights:

1. Over the past 20 years, values of EH and FP have decreased in the urban expansion areas of the GBA, while EH has improved in the peripheral zones. WY achieved high values in built-up areas, and the overall distribution pattern shows a spatial configuration of high values in the east and low values in the west. Rapid urbanization has significantly influenced the spatiotemporal evolution of EH and ESs in the GBA.

2. CS and NPP, as indicators closely linked to fundamental ecosystem metabolism, exhibit a stable synergistic relationship with EH. In the peripheral regions of the GBA, EH and corresponding service functions may be synergistically enhanced, requiring attention to the protection of natural forest lands.

3. FP and WY respond dramatically to human activities and climate change factors, showing a dynamic trade-off relationship with EH. Rapid urbanization has made ecosystems more sensitive to external disturbances, thereby intensifying the trade-off effects between EH and ESs.

By quantifying the relationship characteristics between ecosystem health and service functions, this research, based on the analysis of ecosystem structural and process characteristics, explores how ecosystems respond to external disturbances such as climate change in urban contexts, and further reveals the hierarchical differences in the relationship characteristics between EH and ESs as well as the impact of rapid urbanization on their relationship. This validates existing theories while providing new insights into their real-time manifestations during rapid urbanization. Furthermore, the results offer a theoretical basis for ecosystem management and planning in regions worldwide with similar natural geographic characteristics and urbanization development traits.

Data availability statement

The datasets presented in this article are not readily available because Dataset available on request from the authors. The raw data supporting the conclusions of this article will be made available by the authors on request. Requests to access the datasets should be directed to Y2h1YW5mdXphbmdAMTYzLmNvbQ==.

Author contributions

XW: Conceptualization, Data curation, Investigation, Methodology, Validation, Writing – original draft, Writing – review & editing. YW: Data curation, Investigation, Methodology, Validation, Writing – review & editing. LS: Investigation, Methodology, Software, Writing – review & editing. SD: Funding acquisition, Project administration, Writing – review & editing, Software, Supervision. CZ: Conceptualization, Data curation, Funding acquisition, Investigation, Methodology, Project administration, Resources, Validation, Writing – review & editing.

Funding

The author(s) declare financial support was received for the research and/or publication of this article. This study was supported by the Guangzhou Collaborative Innovation Center on Science-tech of Ecology and Landscape (202206010058); and the National Natural Science Foundation of China (42471026). The funders provided basic data, contributed to improving the manuscript, and offered financial support.

Acknowledgments

We express our gratitude to the Field Scientific Observation and Research Station of Han River in Guangdong province, as well as our colleagues, for their invaluable input and recommendations that have significantly offered helps to improve this manuscript. We would also like to express our sincere gratitude to Xuan Dong and Qianwen Zhan for their assistance with data computation.

Conflict of interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Generative AI statement

The author(s) declare that no Generative AI was used in the creation of this manuscript.

Any alternative text (alt text) provided alongside figures in this article has been generated by Frontiers with the support of artificial intelligence and reasonable efforts have been made to ensure accuracy, including review by the authors wherever possible. If you identify any issues, please contact us.

Publisher’s note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

References

Anselin, L. (1995). Local indicators of spatial association—LISA. Geogr. Anal. 27, 93–115. doi: 10.1111/j.1538-4632.1995.tb00338.x

Crossref Full Text | Google Scholar

Bennett, E. M., Peterson, G. D., and Gordon, L. J. (2009). Understanding relationships among multiple ecosystem services. Ecol. Lett. 12, 1394–1404. doi: 10.1111/j.1461-0248.2009.01387.x

PubMed Abstract | Crossref Full Text | Google Scholar

Box, E. O., Holben, B. N., and Kalb, V. (1989). Accuracy of the AVHRR vegetation index as a predictor of biomass, primary productivity and net CO 2 flux. Vegetatio 80, 71–89. doi: 10.1007/BF00048034

Crossref Full Text | Google Scholar

Chen, T., Feng, Z., Zhao, H., and Wu, K. (2020). Identification of ecosystem service bundles and driving factors in Beijing and its surrounding areas. Sci. Total Environ. 711, 134687. doi: 10.1016/j.scitotenv.2019.134687

PubMed Abstract | Crossref Full Text | Google Scholar

Chen, W., Gu, T., and Zeng, J. (2022). Urbanisation and ecosystem health in the Middle Reaches of the Yangtze River urban agglomerations, China: A U-curve relationship. J. Environ. Manag. 318, 115565. doi: 10.1016/j.jenvman.2022.115565

PubMed Abstract | Crossref Full Text | Google Scholar

Chen, S., Wang, W., Xu, W., Wang, Y., Wan, H., Chen, D., et al. (2018). Plant diversity enhances productivity and soil carbon storage. Proc. Natl. Acad. Sci. 115, 4027–4032. doi: 10.1073/pnas.1700298114

PubMed Abstract | Crossref Full Text | Google Scholar

Chu, X., Zhan, J., Li, Z., Zhang, F., and Qi, W. (2019). Assessment on forest carbon sequestration in the Three-North Shelterbelt Program region, China. J. Clean. Prod. 215, 382–389. doi: 10.1016/j.jclepro.2018.12.296

Crossref Full Text | Google Scholar

Costanza, R. (1992). Toward an operational definition of ecosystem health. Ecosyst. Health: New Goals Environ. Manag. 239, 269. doi: 10.4337/9781035303427.00013

Crossref Full Text | Google Scholar

Costanza, R. (2012). Ecosystem health and ecological engineering. Ecol. Eng. 45, 24–29. doi: 10.1016/j.ecoleng.2012.03.023

Crossref Full Text | Google Scholar

Costanza, R., d’Arge, R., de Groot, R., Farber, S., Grasso, M., Hannon, B., et al. (1997). The value of the world’s ecosystem services and natural capital. Nature 387, 253–260. doi: 10.1038/387253a0

Crossref Full Text | Google Scholar

Darvishi, A., Yousefi, M., Dinan, N. M., and Angelstam, P. (2022). Assessing levels, trade-offs and synergies of landscape services in the Iranian province of Qazvin: towards sustainable landscapes. Landscape Ecol. 37, 305–327. doi: 10.1007/s10980-021-01337-0

Crossref Full Text | Google Scholar

Das, M., Das, A., Pereira, P., and Mandal, A. (2021). Exploring the spatio-temporal dynamics of ecosystem health: A study on a rapidly urbanizing metropolitan area of Lower Gangetic Plain, India. Ecol. Indic. 125, 107584. doi: 10.1016/j.ecolind.2021.107584

Crossref Full Text | Google Scholar

Delphin, S., Escobedo, F., Abd-Elrahman, A., and Cropper, W., Jr. (2013). Mapping potential carbon and timber losses from hurricanes using a decision tree and ecosystem services driver model. J. Environ. Manag. 129, 599–607. doi: 10.1016/j.jenvman.2013.08.029

PubMed Abstract | Crossref Full Text | Google Scholar

Dislich, C., Keyel, A. C., Salecker, J., Kisel, Y., Meyer, K., Auliya, M., et al. (2017). A review of the ecosystem functions in oil palm plantations, using forests as a reference system. Biol. Rev. 92, 1539–1569. doi: 10.1111/brv.12295

PubMed Abstract | Crossref Full Text | Google Scholar

Dou, H., Li, X., Li, S., Dang, D., Li, X., Lyu, X., et al. (2020). Mapping ecosystem services bundles for analyzing spatial trade-offs in inner Mongolia, China. J. Clean. Prod. 256, 120444. doi: 10.1016/j.jclepro.2020.120444

Crossref Full Text | Google Scholar

Feng, Z., Jin, X., Chen, T., and Wu, J. (2021). Understanding trade-offs and synergies of ecosystem services to support the decision-making in the Beijing–Tianjin–Hebei region. Land Use Policy 106, 105446. doi: 10.1016/J.LANDUSEPOL.2021.105446

Crossref Full Text | Google Scholar

Friend, A. D., Lucht, W., Rademacher, T. T., Keribin, R., Betts, R., Cadule, P., et al. (2013). Carbon residence time dominates uncertainty in terrestrial vegetation responses to future climate and atmospheric CO2. Proc. Natl. Acad. Sci. 111, 3280–3285. doi: 10.1073/pnas.1222477110

PubMed Abstract | Crossref Full Text | Google Scholar

Frondoni, R., Mollo, B., and Capotorti, G. (2011). A landscape analysis of land cover change in the Municipality of Rome (Italy): Spatio-temporal characteristics and ecological implications of land cover transitions from 1954 to 2001. Landscape Urban Plan. 100, 117–128. doi: 10.1016/j.landurbplan.2010.12.002

Crossref Full Text | Google Scholar

Gordon, A., Simondson, D., White, M., Moilanen, A., and Bekessy, S. A. (2009). Integrating conservation planning and landuse planning in urban landscapes. Landscape Urban Plan. 91, 183–194. doi: 10.1016/j.landurbplan.2008.12.011

Crossref Full Text | Google Scholar

Grimm, N. B., Faeth, S. H., Golubiewski, N. E., Redman, C. L., Wu, J., Bai, X., et al. (2008). Global change and the ecology of cities. Science 319, 756–760. doi: 10.1126/science.1150195

PubMed Abstract | Crossref Full Text | Google Scholar

Guillaume, T., Kotowska, M. M., Hertel, D., Knohl, A., Krashevska, V., Murtilaksono, K., et al. (2018). Carbon costs and benefits of Indonesian rainforest conversion to plantations. Nat. Commun. 9, 2388. doi: 10.1038/s41467-018-04755-y

PubMed Abstract | Crossref Full Text | Google Scholar

Han, P., Hu, H., Zhou, J., Wang, M., and Zhou, Z. (2024). Integrating key ecosystem services to study the spatio-temporal dynamics and determinants of ecosystem health in Wuhan’s central urban area. Ecol. Indic. 166, 112352. doi: 10.1016/j.ecolind.2024.112352

Crossref Full Text | Google Scholar

Hernández-Blanco, M., Costanza, R., Chen, H., DeGroot, D., Jarvis, D., Kubiszewski, I., et al. (2022). Ecosystem health, ecosystem services, and the well-being of humans and the rest of nature. Global Change Biol. 28, 5027–5040. doi: 10.1111/gcb.16281

PubMed Abstract | Crossref Full Text | Google Scholar

Holling, C. S. (1973). Holling CS. Resilience and stability of ecological systems. Ann Rev Ecol Syst 4: 1-23. Annu. Rev. Ecol. Syst. 4, 1–23. doi: 10.1146/annurev.es.04.110173.000245

Crossref Full Text | Google Scholar

Huang, Y., Gan, X., Feng, Y., Li, J., Niu, S., and Zhou, B. (2024). A new framework for assessing ecosystem health with consideration of the sustainable supply of ecosystem services. Landsc Ecol. 39, 37. doi: 10.1007/s10980-024-01834-y

Crossref Full Text | Google Scholar

Huang, F., Zuo, L., Gao, J., Jiang, Y., Du, F., and Zhang, Y. (2023). Exploring the driving factors of trade-offs and synergies among ecological functional zones based on ecosystem service bundles. Ecol. Indicat. 146, 109827. doi: 10.1016/j.ecolind.2022.109827

Crossref Full Text | Google Scholar

Hwang, C. L. and Yoon, K. P. (1981). Multiple Attribute Decision Making, Lecture Notes in Economics and Mathematical Systems (Berlin, Heidelberg: Springer Berlin Heidelberg). doi: 10.1007/978-3-642-48318-9

Crossref Full Text | Google Scholar

Imhoff, M. L., Bounoua, L., Ricketts, T., Loucks, C., Harriss, R., and Lawrence, W. T. (2004). Global patterns in human consumption of net primary production. Nature 429, 870–873. doi: 10.1038/nature02619

PubMed Abstract | Crossref Full Text | Google Scholar

Jianying, X., Jixing, C., and Yanxu, L. (2020). Partitioned responses of ecosystem services and their tradeoffs to human activities in the Belt and Road region. J. Clean. Prod. 276, 123205. doi: 10.1016/j.jclepro.2020.123205

Crossref Full Text | Google Scholar

Johnson, M. P. (2016). Photosynthesis. Essays Biochem. 60, 255–273. doi: 10.1007/0-306-48137-5

Crossref Full Text | Google Scholar

Kang, P., Chen, W., Hou, Y., and Li, Y. (2018). Linking ecosystem services and ecosystem health to ecological risk assessment: A case study of the Beijing-Tianjin-Hebei urban agglomeration. Sci. Total Environ. 636, 1442–1454. doi: 10.1016/j.scitotenv.2018.04.427

PubMed Abstract | Crossref Full Text | Google Scholar

Li, J., Pei, W., Li, Y., Liu, S., Chen, Y., Wang, B., et al. (2024). Evaluating and diagnosing ecosystem health of the “three-lake” watershed in Yuxi, Yunnan, China from 2010 to 2020 by PSR-KDE. Environ. Res. 258, 119406. doi: 10.1016/j.envres.2024.119406

PubMed Abstract | Crossref Full Text | Google Scholar

Li, X., Yu, K., Xu, G., Li, P., Li, Z., Shi, P., et al. (2025). Quantifying thresholds of key drivers for ecosystem health in large-scale river basins: A case study of the upper and middle Yellow River. J. Environ. Manage. 383, 125480. doi: 10.1016/j.jenvman.2025.125480

PubMed Abstract | Crossref Full Text | Google Scholar

Liu, N., Sun, P., Caldwell, P., Harper, R., Liu, S., and Sun, G. (2020). Trade-off between watershed water yield and ecosystem productivity along elevation gradients on a complex terrain in southwestern China. J. Hydrol. 590, 125449. doi: 10.1016/j.jhydrol.2020.125449

Crossref Full Text | Google Scholar

Liu, Z., Wang, S., and Fang, C. (2023). Spatiotemporal evolution and influencing mechanism of ecosystem service value in the Guangdong-Hong Kong-Macao Greater Bay Area. J. Geogr. Sci. 33, 1226–1244. doi: 10.1007/s11442-023-2127-5

Crossref Full Text | Google Scholar

Liu, Y., Xu, G., Jiang, H., Hu, X., Peng, J., Song, Z., et al. (2015). Synergy between ecosystem services and ecosystem health in the forest area of Northeast China. Prog. Geogr. 34, 761–771. doi: 10.18306/dlkxjz.2015.06.011

Crossref Full Text | Google Scholar

Liu, W., Zhan, J., Zhao, F., Yan, H., Zhang, F., and Wei, X. (2019). Impacts of urbanization-induced land-use changes on ecosystem services: A case study of the Pearl River Delta Metropolitan Region, China. Ecol. Indic. 98, 228–238. doi: 10.1016/j.ecolind.2018.10.054

Crossref Full Text | Google Scholar

Lundberg, S. M., Nair, B., Vavilala, M. S., Horibe, M., Eisses, M. J., Adams, T., et al. (2018). Explainable machine-learning predictions for the prevention of hypoxaemia during surgery. Nat. Biomed. Eng. 2, 749–760. doi: 10.1038/s41551-018-0304-0

PubMed Abstract | Crossref Full Text | Google Scholar

Lv, T., Zeng, C., Lin, C., Liu, W., Cheng, Y., and Li, Y. (2023). Towards an integrated approach for land spatial ecological restoration zoning based on ecosystem health assessment. Ecol. Indicat. 147, 110016. doi: 10.1016/j.ecolind.2023.110016

Crossref Full Text | Google Scholar

Millennium Ecosystem Assessment (2005). Ecosystems and human well-being (Washington, DC: Island press).

Google Scholar

Mo, W., Zhao, Y., Yang, N., and Xu, Z. (2023). Ecological function zoning based on ecosystem service bundles and trade-offs: a study of Dongjiang Lake Basin, China. Environ. Sci. pollut. Res. 30, 40388–40404. doi: 10.1007/s11356-022-24782-z

PubMed Abstract | Crossref Full Text | Google Scholar

Muradian, R. (2001). Ecological thresholds: a survey. Ecol. Econ 38, 7–24. doi: 10.1016/S0921-8009(01)00146-X

Crossref Full Text | Google Scholar

Nemani, R. R., Keeling, C. D., Hashimoto, H., Jolly, W. M., Piper, S. C., Tucker, C. J., et al. (2003). Climate-driven increases in global terrestrial net primary production from 1982 to 1999. Science 300, 1560–1563. doi: 10.1126/science.1082750

PubMed Abstract | Crossref Full Text | Google Scholar

Pan, Z., He, J., Liu, D., Wang, J., and Guo, X. (2021). Ecosystem health assessment based on ecological integrity and ecosystem services demand in the Middle Reaches of the Yangtze River Economic Belt, China. Sci. Total Environ. 774, 144837. doi: 10.1016/j.scitotenv.2020.144837

Crossref Full Text | Google Scholar

Peng, J., Liu, Y., Li, T., and Wu, J. (2017). Regional ecosystem health response to rural land use change: A case study in Lijiang City, China. Ecological Indicators 72, 399–410. doi: 10.1016/j.ecolind.2016.08.024

Crossref Full Text | Google Scholar

Peng, J., Liu, Y., Wu, J., Lv, H., and Hu, X. (2015). Linking ecosystem services and landscape patterns to assess urban ecosystem health: A case study in Shenzhen City, China. Landscape Urban Plan. 143, 56–68. doi: 10.1016/j.landurbplan.2015.06.007

Crossref Full Text | Google Scholar

Qi, Y., Lian, X., Wang, H., Zhang, J., and Yang, R. (2020). Dynamic mechanism between human activities and ecosystem services: A case study of Qinghai lake watershed, China. Ecol. Indic. 117, 106528. doi: 10.1016/j.ecolind.2020.106528

Crossref Full Text | Google Scholar

Rapport, D. J., Böhm, G., Buckingham, D., Cairns, J., Jr., Costanza, R., Karr, J., et al. (1999). Ecosystem health: the concept, the ISEH, and the important tasks ahead. Ecosyst. Health 5, 82–90. doi: 10.1046/j.1526-0992.1999.09913.x

Crossref Full Text | Google Scholar

Rapport, D. J., Costanza, R., and McMichael, A. J. (1998). Assessing ecosystem health. Trends Ecol. Evol. 13, 397–402. doi: 10.1016/S0169-5347(98)01449-9

PubMed Abstract | Crossref Full Text | Google Scholar

Reader, M. O., Eppinga, M. B., de Boer, H. J., Damm, A., Petchey, O. L., and Santos, M. J. (2023). Biodiversity mediates relationships between anthropogenic drivers and ecosystem services across global mountain, island and delta systems. Global Environ. Change 78, 102612. doi: 10.1016/j.gloenvcha.2022.102612

Crossref Full Text | Google Scholar

Redhead, J., Stratford, C., Sharps, K., Jones, L., Ziv, G., Clarke, D., et al. (2016). Empirical validation of the InVEST water yield ecosystem service model at a national scale. Sci. Total Environ. 569, 1418–1426. doi: 10.1016/j.scitotenv.2016.06.227

PubMed Abstract | Crossref Full Text | Google Scholar

Runting, R. K., Bryan, B. A., Dee, L. E., Maseyk, F. J., Mandle, L., Hamel, P., et al. (2017). Incorporating climate change into ecosystem service assessments and decisions: a review. Global Change Biol. 23, 28–41. doi: 10.1111/gcb.13457

PubMed Abstract | Crossref Full Text | Google Scholar

Shao, S., Pan, Y., and Huang, A. (2025). Delineating ecological improvement zones based on impact thresholds of ecosystem health changes in Northwest Hubei, China. Ecol. Indic. 178, 113898. doi: 10.1016/j.ecolind.2025.113898

Crossref Full Text | Google Scholar

Sharp, R., Chaplin-Kramer, R., Wood, S., Guerry, A., Tallis, H., Ricketts, T., et al. (2018). InVEST User’s Guide. Stanford, CA, USA: Natural Capital Project.

Google Scholar

Smith, N. G., Zhu, Q., Keenan, T., and Riley, W. (2024). Acclimation of photosynthesis to CO2 increases ecosystem carbon storage due to leaf nitrogen savings. Global Change Biol. 30, e17558. doi: 10.1111/gcb.17558

PubMed Abstract | Crossref Full Text | Google Scholar

Su, H., Du, M., Liu, Q., Kang, X., Zhao, L., Zheng, W., et al. (2024). Assessment of regional Ecosystem Service Bundles coupling climate and land use changes. Ecol. Indic. 169, 112844. doi: 10.1016/j.ecolind.2024.112844

Crossref Full Text | Google Scholar

Sun, G., Zhou, G., Zhang, Z., Wei, X., McNulty, S. G., and Vose, J. M. (2006). Potential water yield reduction due to forestation across China. J. Hydrol. 328, 548–558. doi: 10.1016/j.jhydrol.2005.12.013

Crossref Full Text | Google Scholar

Walker, B., Holling, C. S., Carpenter, S. R., and Kinzig, A. P. (2004). Resilience, adaptability and transformability in social-ecological systems. Ecol. Soc. 9, 5. doi: 10.5751/es-00650-090205

Crossref Full Text | Google Scholar

Walther, S., Duveiller, G., Jung, M., Guanter, L., Cescatti, A., and Camps-Valls, G. (2019). Satellite observations of the contrasting response of trees and grasses to variations in water availability. Geophys. Res. Lett. 46, 1429–1440. doi: 10.1029/2018GL080535

Crossref Full Text | Google Scholar

Wang, X., Liu, L., and Zhang, S. (2021). Integrated model framework for the evaluation and prediction of the water environmental carrying capacity in the Guangdong–Hong Kong–Macao Greater Bay Area. Ecol. Indicat. 130, 108083. doi: 10.1016/j.ecolind.2021.108083

Crossref Full Text | Google Scholar

Wang, J., Peng, J., and Fu, B. (2023b). Integrated ecological protection and restoration in the Guangdong–Hong Kong–Macao Greater bay area: thoughts and suggestions. Bull. Chin. Acad. Sci. (Chinese Version) 38, 288–293. doi: 10.16418/j.issn.1000-3045.20221221002

Crossref Full Text | Google Scholar

Wang, Y., Yang, Z., Yu, M., Lin, R., Zhu, L., and Bai, F. (2023a). Integrating ecosystem health and services for assessing ecological risk and its response to typical land-use patterns in the eco-fragile region, North China. Environ. Manage. 71, 867–884. doi: 10.1007/s00267-022-01742-4

PubMed Abstract | Crossref Full Text | Google Scholar

Wang, J., Zhou, W., Pickett, S. T., Yu, W., and Li, W. (2019). A multiscale analysis of urbanization effects on ecosystem services supply in an urban megaregion. Sci. Total Environ. 662, 824–833. doi: 10.1016/j.scitotenv.2019.01.260

PubMed Abstract | Crossref Full Text | Google Scholar

White, P. S. and Pickett, S. T. A. (1985). “Natural disturbance and patch dynamics: an introduction,” in The Ecology of Natural Disturbance and Patch Dynamics (New York, America: Elsevier), 3–13. doi: 10.1016/B978-0-08-050495-7.50006-5

Crossref Full Text | Google Scholar

Wu, M., Wu, J., and Zang, C. (2021). A comprehensive evaluation of the eco-carrying capacity and green economy in the Guangdong–Hong Kong–Macao Greater Bay Area, China. J. Clean. Prod. 281, 124945. doi: 10.1016/j.jclepro.2020.124945

Crossref Full Text | Google Scholar

Wu, W., Zeng, H., Guo, C., You, W., Xu, H., Hu, Y., et al. (2024). Spatial heterogeneity and management challenges of ecosystem service trade-offs: a case study in Guangdong Province, China. Environ. Manage. 73, 378–394. doi: 10.1007/s00267-023-01851-8

PubMed Abstract | Crossref Full Text | Google Scholar

Xia, H., Yuan, S., and Prishchepov, A. V. (2023). Spatial-temporal heterogeneity of ecosystem service interactions and their social-ecological drivers: Implications for spatial planning and management. Resour. Conserv. Recycl. 189, 106767. doi: 10.1016/j.resconrec.2022.106767

Crossref Full Text | Google Scholar

Xu, D., Huang, Z., and Huang, R. (2019). The spatial effects of haze on tourism flows of Chinese cities: Empirical research based on the spatial panel econometric model. Acta Geographica Sinica 74, 04000814. doi:  doi: 10.11821/dlxb201904014

Crossref Full Text | Google Scholar

Xue, C., Chen, X., Xue, L., Zhang, H., Chen, J., and Li, D. (2023). Modeling the spatially heterogeneous relationships between tradeoffs and synergies among ecosystem services and potential drivers considering geographic scale in Bairin Left Banner, China. Sci. Total Environ. 855, 158834. doi: 10.1016/j.scitotenv.2022.158834

PubMed Abstract | Crossref Full Text | Google Scholar

Yang, C., Chen, M., and Yuan, Q. (2021). The application of XGBoost and SHAP to examining the factors in freight truck-related crashes: An exploratory analysis. Accident Anal. Prev. 158, 106153. doi: 10.1016/j.aap.2021.106153

PubMed Abstract | Crossref Full Text | Google Scholar

Yang, J. and Huang, X. (2021). The 30 m annual land cover dataset and its dynamics in China from 1990 to 2019 [dataset. Earth Syst. Sci. Data 13, 3907–3925. doi: 10.5194/essd-13-3907-2021

Crossref Full Text | Google Scholar

Ye, X., Skidmore, A. K., and Wang, T. (2013). Within-patch habitat quality determines the resilience of specialist species in fragmented landscapes. Landscape Ecol. 28, 135–147. doi: 10.1007/s10980-012-9826-0

Crossref Full Text | Google Scholar

Yu, W., Zhang, D., Liao, J., Ma, L., Zhu, X., Zhang, W., et al. (2022). Linking ecosystem services to a coastal bay ecosystem health assessment: A comparative case study between Jiaozhou Bay and Daya Bay, China. Ecol. Indic. 135, 108530. doi: 10.1016/j.ecolind.2021.108530

Crossref Full Text | Google Scholar

Yuan, Y., Guo, W., Tang, S., and Zhang, J. (2024). Effects of patterns of urban green-blue landscape on carbon sequestration using XGBoost-SHAP model. J. Clean. Prod. 476, 143640. doi: 10.1016/j.jclepro.2024.143640

Crossref Full Text | Google Scholar

Yushanjiang, A., Zhou, W., Wang, J., and Wang, J. (2024). Impact of urbanization on regional ecosystem services— a case study in Guangdong-Hong Kong-Macao Greater Bay Area. Ecol. Indic. 159, 111633. doi: 10.1016/j.ecolind.2024.111633

Crossref Full Text | Google Scholar

Zhang, H., Deng, W., Zhang, S., Peng, L., and Liu, Y. (2022). Impacts of urbanization on ecosystem services in the Chengdu-Chongqing Urban Agglomeration: Changes and trade-offs. Ecol. Indic. 139, 108920. doi: 10.1016/j.ecolind.2022.108920

Crossref Full Text | Google Scholar

Zhang, Y., Liu, Y., Zhang, Y., Liu, Y., Zhang, G., and Chen, Y. (2018). On the spatial relationship between ecosystem services and urbanization: A case study in Wuhan, China. Sci. Total Environ. 637, 780–790. doi: 10.1016/j.scitotenv.2018.04.396

PubMed Abstract | Crossref Full Text | Google Scholar

Zhang, P., Zhou, Y., Xie, Y., Wang, Y., Li, B., Li, B., et al. (2021). Assessment of the water-energy-food nexus under spatial and social complexities: a case study of Guangdong-Hong Kong-Macao. J. Environ. Manag. 299, 113664. doi: 10.1016/j.jenvman.2021.113664

PubMed Abstract | Crossref Full Text | Google Scholar

Zhou, Y., Shan, Y., Liu, G., and Guan, D. (2018). Emissions and low-carbon development in Guangdong–Hong Kong–Macao Greater Bay Area cities and their surroundings. Appl. Energy 228, 1683–1692. doi: 10.1016/j.apenergy.2018.07.038

Crossref Full Text | Google Scholar

Zhu, Q. and Cai, Y. (2023). Integrating ecological risk, ecosystem health, and ecosystem services for assessing regional ecological security and its driving factors: Insights from a large river basin in China. Ecol. Indic. 155, 110954. doi: 10.1016/j.ecolind.2023.110954

Crossref Full Text | Google Scholar

Zhu, Q. and Cai, Y. L. (2024). Impact of ecological risk and ecosystem health on ecosystem services. Acta Geogr. Sin. 79, 1303–1317. doi: 10.11821/dlxb202405013

Crossref Full Text | Google Scholar

Zhu, Q., Guo, J., Guo, X., Chen, L., Han, Y., and Liu, S. (2021). Relationship between ecological quality and ecosystem services in a red soil hilly watershed in southern China. Ecol. Indicat. 121, 107119. doi: 10.1016/j.ecolind.2020.107119

Crossref Full Text | Google Scholar