ORIGINAL RESEARCH article

Front. Environ. Sci., 28 September 2026

Sec. Ecosystem Restoration

Volume 14 - 2026 | https://doi.org/10.3389/fenvs.2026.1912598

Remotely sensed vegetation resilience in the four major estuarine ecosystems of eastern China

Abstract

Introduction:

Estuarine wetlands are among the most productive ecosystems on Earth, but they have been experiencing severe degradation induced by human activities and climate change worldwide. Traditional assessments of degradation mostly rely on static indicators, which are insufficient to reflect dynamical ecosystem properties such as systemic resilience.

Methods:

Here, we used remote sensing data and critical slowing down indicators to evaluate ecological resilience of vegetation in four major estuaries in eastern China: the Liaohe Estuary, the Yangtze Estuary, the Minjiang Estuary, and the Pearl Estuary. We identified hotspots of resilience loss during the past two decades, and then assessed relative ecological risk using the Pressure–State–Response (PSR) framework. Using random forest modeling and SHAP (SHapley Additive exPlanations) analysis, we further identified key factors associated with the classification of resilience-loss hotspots.

Results:

Our results showed that all four estuaries have undergone resilience loss, with resilience loss areas accounting for ∼40% of the total study area. Moreover, hotspots of resilience loss were found to be mostly/strongly/largely associated with faster sea level rise and macroclimate conditions, but showed a spatial mismatch with relatively high ecological risk areas. Defined as resilience-loss hotspots with relatively low ecological risk, potential risk areas accounted for ∼35% of the total area. Despite still exhibiting relatively favorable environmental conditions, these areas might degrade quickly under ongoing unprecedented environmental change and thus need further monitoring.

Discussion:

Our study combined risk assessment with resilience evaluation to provide an integrated framework detecting unrecognized ecosystem states. The framework contributes to better understanding resilience dynamics of estuarine wetlands, and has useful implications for adaptive management of these important ecosystems in a changing world.

1 Introduction

Estuarine ecosystems are among the most intensively utilized yet highly vulnerable ecosystems worldwide (; Pein and Staneva, 2026). Situated at the interface among terrestrial, marine, and riverine ecosystems, estuarine wetlands are shaped by multiple interacting processes, including tidal dynamics, sediment transport, and fresh-salt water exchanges (; ). These complex ecological and geomorphological processes create highly heterogeneous landscapes that support a wide range of ecosystem services (Rodrigues-Filho et al., 2023; ). However, due to climate change and anthropogenic stresses such as reclamation and species invasion, estuarine ecosystems are undergoing a severe and irreversible process of degradation (; van Maren et al., 2025). A recent study estimated that the global wetlands, including estuarine wetlands, has lost approximately 22% of its area since 1970 (). The estuarine wetlands have exhibited not only landscape type fluctuation, but also fragmentation and biogeochemical cycle change (; ), resulting in substantial declines in biodiversity, reduced carbon storage capacity, weakened coastal protection, and diminished ecosystem-service provision (Zhang et al., 2025; ).

China’s estuarine wetlands are facing serious challenges as induced by reclamation, urbanization, and biological invasion. A national scale study indicated that from 1984 to 2018, the area of coastal wetlands in China generally showed a downward trend, and reclamation was one of the main driving factors causing the loss of natural wetlands (Wang et al., 2014; Wang et al., 2021). In recent years, the implementation of protection and restoration policies has led to partial recovery of wetland area in some regions, but the restoration of ecosystem structure and functions still faces challenges, such as incomplete recovery of biodiversity and ecosystem services (), delayed recovery of soil carbon and biogeochemical functions (Wang B. et al., 2026), and ecological trade-offs associated with invasive species (Zhao et al., 2026). These developments underscore the urgency of monitoring degradation in coastal and estuarine wetlands before irreversible transitions occur.

Specifically, along the eastern coast of China, Liaohe Estuary, Yangtze Estuary, Minjiang Estuary, and the Pearl River Estuary are the four most representative estuaries that have experienced substantial ecological degradation. In Liaohe Estuary, agriculture activities, aquaculture expansion and land-use change have led to salt marsh degradation and habitat fragmentation. The expansion of Spartina alterniflora and the conversion of wetlands to artificial ponds were identified as major drivers of wetland degradation (). The Yangtze Estuary, located in one of the most densely populated and economically developed regions of China, has been increasingly affected by intensive reclamation, port construction, and urban expansion. Consequently, the extent of natural wetlands has continuously declined, leading to substantial changes in ecosystem structure and functions (). The Minjiang Estuary has been influenced by reclamation, aquaculture activities, and pollutant inputs, resulting in mangrove degradation and reduced habitat connectivity. The Pearl River Estuary, as a typical highly urbanized area, has also experienced increasing ecological pressures associated with industrialization and urban expansion. Mangroves and intertidal wetlands have been continuously threatened, manifested by habitat loss, ecosystem fragmentation, and the reduction of natural shorelines (Wei et al., 2021).

In order to shift from a static assessment of the ecosystem’s state to a dynamic evaluation that assesses resilience (i.e., the capacity of maintaining stable in the face of disturbances), a series of methods have been developed (Scheffer et al., 2001; Scheffer et al., 2012). Theoretical and empirical studies suggest that environmental change and external disturbances can drive ecosystems from one stable state to another, for example, from a high-biomass, densely vegetated state to a low-biomass, sparsely vegetated state, or vice versa (). Such regime shifts are often characterized by hysteresis and limited reversibility, making recovery difficult even after the original environmental conditions are restored. A key prediction of resilience theory is that ecosystems approaching a critical transition exhibit a characteristic dynamical behavior known as critical slowing down (CSD, Scheffer et al., 2009). As approaching tipping points, ecosystems recover more slowly from perturbations, resulting in increased temporal autocorrelation and variance in ecosystem state variables, which could be treated as a proxy for resilience change. This is of great significance for understanding the state of the ecosystem and making reasonable predictions about the impending tipping point.

Although resilience assessments provide important information on intrinsic stability of ecosystems, they capture only one dimension of ecosystem vulnerability and need to be complemented by information on external pressures and current ecological conditions. Conventional ecological risk assessments, including the Pressure-State-Response (PSR) model, provide such complementary information by evaluating the pressures imposed on ecosystems, their current states, and their responses to these pressures. However, risk assessments do not characterize the intrinsic stability of ecosystems, meaning that ecosystems exposed to similar levels of pressure may differ substantially in their probability of future degradation. On the other hand, declining resilience does not necessarily indicate that an ecosystem has already undergone degradation, particularly when external pressures remain relatively low, ecological conditions are favorable, and response capacity remains strong. Therefore, integrating ecological risk assessment with resilience measurement could provide a more comprehensive characterization of ecosystem vulnerability by jointly considering both current risk conditions and changes in intrinsic stability. Importantly, this integration helps to identify areas that currently exhibit relatively low ecological risk but have already experienced resilience loss. Such areas may represent potential risk areas where reduced resilience would decrease the resistance to future disturbances and increase the likelihood of catastrophic shifts to alternative ecosystem states (Pichon et al., 2026).

A large number of studies have shown that climate change is already driving widespread resilience loss across terrestrial ecosystems. For example, tropical and northern forests, which are crucial components of the Earth system, are experiencing substantial resilience loss and have shifted from a nearly closed state to an open one, thereby triggering profound social-ecological impacts (Verbesselt et al., 2016; ; Rotbarth et al., 2023). Once a dryland ecosystem crosses the drought threshold, it will rapidly transform into a state without vegetation cover, leading to a further expansion of the arid area (Wang R. et al., 2026). Despite growing interest in ecosystem resilience, studies focusing on estuarine wetlands remain extremely limited. Although increasing evidence suggests that climate change, sea level rise, and accelerating urbanization are progressively eroding the resilience of estuarine ecosystems, quantitative assessments of the spatial extent, geographical hotspots, and primary drivers of resilience loss remain scarce, hindering our ability to identify vulnerable regions and anticipate future ecological tipping points. Meanwhile, it is challenging to assess resilience change and delineate main factors driving resilience change in estuarine wetlands (Moffett et al., 2015). Different from terrestrial or aquatic ecosystems, estuarine wetlands are situated at the land-sea interface and are shaped by multiple processes originating from both land and ocean. This imparts greater spatiotemporal heterogeneity to estuarine wetlands, introducing higher levels of spatiotemporal noise signals (Wang C. et al., 2025). Consequently, special attention must be paid to filtering time-series noise and considering spatial scales when evaluating ecological risks and assessing resilience dynamics. Furthermore, driven by multiple stressors, like climate change, fresh-saline water exchange dynamics, tidal dynamics, sedimentary processes, and human activities, the underlying mechanisms driving changes in estuarine wetland resilience are exceptionally complex, posing significant challenges for ecosystem regulation and management (Mitchell et al., 2015).

Here, we assessed ecological resilience in four major estuarine wetlands along the eastern coast of China using satellite-based vegetation index time series based on critical slowing down, and quantified spatial patterns of resilience change and hotspots of resilience loss. Then, we developed a PSR (Pressure-State-Response)-based ecological risk assessment framework and compared ecological risk patterns with resilience-loss hotspots to identify potential risk areas where signals of resilience decline were detected despite seemingly low ecological risk. Finally, to provide scientific support for future conservation interventions and adaptive management, we used random forest modeling and SHAP (SHapley Additive exPlanations) analysis to examine environmental and anthropogenic factors associated with resilience-loss hotspot classification.

2 Materials and methods

2.1 Datasets and preprocessing

The definition of the study area was derived from the National Coastal Wetland Map of China (Zhang et al., 2026). We extracted consistent wetland areas from 2000 to 2024 and clipped them to the boundaries of four major estuaries along the eastern coast of China - the Liaohe Estuary, Yangtze Estuary, Minjiang Estuary, and Pearl River Estuary (Figure 1). The original dataset distinguishes among salt marshes, mangroves, and tidal flats. Although there exist specific differences in dynamics and ecological processes among the three types, we conducted analysis based on the assumption that estuarine wetlands exhibit alternative stable states and that transitions among tidal flats, salt marshes, and mangroves can occur naturally under changing environmental conditions. Therefore, all three wetland types were treated as a unified wetland category in this study. To test the robustness of the corresponding results, we also extracted invariable salt marsh, mangrove, tidal flat and conducted similar analysis (Supplementary Figure S8).

FIGURE 1

Vegetation dynamics were characterized using the Normalized Difference Vegetation Index (NDVI) derived from the MODIS/061/MOD13Q1 product (16-day composite, 250m spatial resolution, ). This product is generated from a single sensor, thereby minimizing potential data inconsistencies and temporal drift associated with sensor transitions. To ensure data quality, observations flagged in the SummaryQA layer as snow- or ice-covered pixels, cloudy pixels, or marginal-quality data were excluded from the analysis. Considering that estuarine wetlands are frequently inundated, we evaluated the robustness of our results under different screening thresholds for the original NDVI values. Specifically, we compared three scenarios: retaining all valid observations, excluding NDVI values below 0, and excluding NDVI values below 0.2 (a commonly used threshold for vegetated surfaces, Tucker, 1979). The results were largely consistent across these scenarios, indicating that the findings were robust to the choice of threshold (Supplementary Figure S2C). Consequently, to maximize the length of the available time series and minimize data gaps that could affect subsequent analyses, no low-value filtering was applied to the original NDVI observations.

Reliable estimation of resilience indicators requires stationary vegetation time series. Therefore, long-term trends and seasonal signals were removed prior to resilience analysis. Following previous studies, we tried three commonly used methods for time series decomposition, and adopted a combined method in which long-term trends were removed using a rolling-window mean and seasonal components were modeled and removed using a third-order harmonic regression based on the robustness test on the synthetic time series (Smith and Boers, 2023). Although the STL method appeared to capture the residual trend more closely (Supplementary Figure S1A), the method adopted in this study showed better performance in recovering the underlying lag-1 autocorrelation trend relevant to resilience assessment (Supplementary Figure S1B,C). We also did a robust test based on the acquired NDVI time series and it showed that the results were roughly robust under different decomposition methods (Supplementary Figure S2A). All subsequent analyses were conducted on the residual NDVI time series after the removal of both long-term trends and seasonal variability.

2.2 Ecological resilience based on critical slowing down indicators

We assessed ecological resilience of each estuary ecosystem by computing one of the most robust resilience loss signals, the lag-1 autocorrelation (AC1):where:

The metric was calculated separately using 1-year, 3-year, 5-year moving windows applied to residual NDVI time series to evaluate the robustness of the signals to moving-window size. The resulting AC1 trends were highly consistent across window sizes, indicating that the choice of moving-window length did not substantially alter the results (Supplementary Figure S2B). Considering the high level of environmental noise in estuarine wetlands caused by multiple interacting stressors, the 5-year moving window was adopted for all subsequent analyses because it provides more stable estimates of temporal dynamics. To statistically evaluate the temporal trends in these windowed sequences, we employed a robust Mann-Kendall test and calculated Kendall’s tau rank correlation coefficient (). Significant increasing trends in AC1 were interpreted as ecological resilience losing indicators based on the theory of CSD (Scheffer et al., 2009). Although the temporal trend in variance has also been proposed as an indicator of ecological resilience loss, it was not considered in this study because variance is not expected to increase universally before critical transitions. Under certain conditions, such as when environmental forcing is predominantly stochastic or when critical slowing down weakens the system’s ability to respond to high-frequency environmental fluctuations, variance may even decline prior to a transition ().

In addition, we conducted a robust test to analysis whether abrupt changes in NDVI residual would influence the signals. The abrupt changes were detected following a previous study (Smith and Boers, 2023). However, we did not require the mean values of the residual series before and after abrupt changes to show no significant difference. Instead, significant changes in either direction, including both increases and decreases, were considered and detected. The differences in AC1 trends between pixels with and without abrupt changes were evaluated using the Mann-Whitney U test (Supplementary Figure S3), with p-values adjusted using the Benjamini–Hochberg (BH) procedure to control the false discovery rate associated with multiple comparisons ().

To identify spatial clusters of resilience loss, pixels with a positive Kendall’s tau rank correlation coefficient and a statistically significant trend (p < 0.05) were first classified as ecological resilience-loss pixels. This binary variable was then used to calculate the local Getis–Ord hotspot statistic:

Here, denotes the binary resilience-loss indicator of grid j (1 = resilience-loss pixel; 0 = otherwise), is the spatial weight between grid i and j. n is the total number of grids, and are the mean and standard deviation of the binary resilience-loss variable across all grids, respectively. To minimize the uncertainty associated with selecting a single neighborhood scale, the statistic was calculated under three neighborhood radii (r = 1,2, and 3 grid cells). Based on the two-tailed significance test of the statistic, all grids were classified as hotspots or coldspots (at the 90%, 95%, and 99% confidence levels when p < 0.1, 0.05, or 0.01) or as non-significant areas. The resulting p-values were adjusted using the Benjamini–Hochberg false discovery rate (FDR) procedure to account for multiple comparisons (). For a set of valid pixels, the p values were sorted as:

The adjusted p value for the -th ordered test was calculated as:

Pixels were classified as FDR-corrected hotspots when and , corresponding to 90%, 95%, and 99% confidence levels, respectively. Pixels were classified as FDR-corrected coldspots when and passed the same thresholds. This correction was performed separately for each estuary, and neighborhood radius.

Since ecological resilience-loss pixels had already been identified through the initial Mann–Kendall significance test, the objective of the hotspot analysis was not to reassess whether resilience loss was statistically significant, but to identify whether these significant resilience-loss pixels formed spatial clusters. Accordingly, a grid was classified as a resilience-loss hotspot/coldspot if it was detected as a significant hotspot/coldspot after FDR correction under at least one of the three neighborhood radii. The confidence level assigned to each hotspot/coldspot corresponded to the highest significance level achieved across the tested neighborhood radii. This strategy inevitably increases the likelihood of including isolated or scale-sensitive clusters that arise from random variation or local noise. However, it substantially reduces the omission of small but potentially meaningful resilience-loss hotspots. We considered this trade-off preferable for practical monitoring and ecosystem management because overlooking genuinely vulnerable areas may delay the detection of impending regime shifts and ultimately lead to substantially higher restoration costs than falsely identifying a limited number of noisy hotspots.

To evaluate the sensitivity of hotspot detection to neighborhood scale and multiple-testing correction, we compared the proportions of hotspot and coldspot areas obtained using the original significance test and the FDR-corrected results under different neighborhood radii. We further quantified the spatial consistency of hotspot and coldspot patterns by calculating the Jaccard similarity coefficient between the results from each neighborhood radius and those obtained using the smallest neighborhood radius (r = 1) as the reference (Supplementary Figure S5; ). Finally, a spatial permutation test was performed by randomly permuting the resilience-loss pixels while preserving the total number of positive pixels and recalculating the statistic repeatedly (n = 1,000). Empirical p-values were then derived by comparing the observed numbers of hotspots and coldspots with those expected under random spatial configurations (Supplementary Figure S6).

2.3 Environmental covariates and risk classification

A pressure-state-response (PSR) framework was used to construct a covariate-based ecological risk index for estuarine wetlands (). The PSR framework organizes environmental indicators into three subsystems: external pressures, current ecological state, and ecological response or buffering capacity. All the covariates were prepared using Google Earth Engine (GEE) for each estuarine region. Covariate values were extracted at the centroid of each 250-m sampling grid. Because several climate and oceanographic datasets have coarser spatial resolutions than the analysis grid, neighboring sampling grids may receive identical covariate values. We consider this appropriate because these variables characterize regional environmental forcing that is inherently spatially smooth over the extent of individual estuaries. Our objective was to assess ecological risk for each 250-m vegetation grid while incorporating environmental drivers at their native spatial resolutions, rather than artificially resampling coarse-resolution datasets to create spurious fine-scale variability. The covariates were also used to analyze the influence of external factors on the resilience changes of vegetation in estuaries and were grouped into six categories (Table 1; Supplementary Figure S11).

TABLE 1

VariableDatasetRaw temporal resolutionRaw spatial resolutionPSR typeExplanatory type
NightlightVisible Infrared Imaging Radiometer Suite (VIIRS) Day/Night Band product ()Monthly; 2012–2024500 mPressureHuman activity
Temperature meanThe Fifth generation of ECMWF atmospheric reanalyses for climate (ERA5) Land monthly aggregated dataset ()Monthly; 2000–20240.1°PressureClimate
Temperature varianceThe Fifth generation of ECMWF atmospheric reanalyses for climate (ERA5) Land monthly aggregated dataset ()Monthly; 2000–20240.1°PressureClimate
Precipitation meanThe Fifth generation of ECMWF atmospheric reanalyses for climate (ERA5) Land monthly aggregated dataset ()Monthly; 2000–20240.1°PressureClimate
Precipitation varianceThe Fifth generation of ECMWF atmospheric reanalyses for climate (ERA5) Land monthly aggregated dataset ()Monthly; 2000–20240.1°PressureClimate
Sea level meanCopernicus Marine Environment Monitoring Service ()Monthly; 2000–20241/12°PressureSea level
Sea level varianceCopernicus Marine Environment Monitoring Service ()Monthly; 2000–20241/12°PressureSea level
Inundation durationMODIS/061/MOD09A1 (Vermote, 2021)8-day composite; 2000–2024500 mPressureHydrology
NDVI meanMODIS/061/MOD13Q1 NDVI mean ()16-day composite; 2000–2024250 mStateVegetation
MNDWI P10MODIS/061/MOD09A1 (Vermote, 2021)8-day composite; 2000–2024500 mStateHydrology
ElevationShuttle Radar Topography Mission (SRTM) digital elevation model (Reuter et al., 2007)Static30 mStateTerrain
SlopeShuttle Radar Topography Mission (SRTM) digital elevation model (Reuter et al., 2007)Static30 mStateTerrain
NDVI residual Theil-Sen slopeMODIS/061/MOD13Q1 NDVI 90th percentile ()16-day composite; 2000–2024250 mResponse/
NDVI residual varianceMODIS/061/MOD13Q1 NDVI 90th percentile ()16-day composite; 2000–2024250 mResponse/

Variables used in risk assessment and explanatory model.

Predictor variables were derived from long-term remote-sensing and reanalysis records. For time-varying variables, mean denotes the long-term average over the study period, variance denotes the long-term temporal variance, P10 denotes the 10th percentile. MNDWI, was calculated from MODIS, surface reflectance as (green-SWIR1)/(green+SWIR1), and inundation duration was calculated as the fraction of clear MODIS, observations with MNDWI>0. NDVI, residuals were obtained from the MODIS NDVI, time series after removing the fitted seasonal/trend component. The NDVI, residual Theil-Sen slope represents the robust temporal trend of the residual series, and NDVI, residual variance represents the temporal variance of the residual series.

All indicators were transformed so that larger values represented higher ecological risk. For positive risk indicators, such as nighttime light, temperature, and sea level, higher raw values were interpreted as higher risk. For negative risk indicators, such as NDVI, MNDWI, and elevation, lower raw values were interpreted as higher risk. Within each estuary, each indicator was rescaled to the range 0–1 using min-max normalization. For positive indicators, the normalized risk value was calculated as:

For negative indicators, the normalized value was reversed:where is the normalized risk value of grid cell for indicator , and is the corresponding raw covariate value.

Indicator weights were estimated using the entropy weight method separately for each estuary. For each indicator, the proportional contribution of grid cell was calculated as:

The entropy value of indicator was then computed as:

And the information utility value was:

The final entropy weight of each indicator was:

Within each PSR subsystem, normalized indicators were combined using their entropy weights to obtain pressure, state, and response risk scores. The three subsystem scores were then weighted again using the entropy weight method to calculate the composite PSR ecological risk score:

Here, is the subsystem risk score of grid cell for subsystem , and is the entropy-derived weight of the corresponding subsystem. The resulting PSR ecological risk score is a continuous index ranging from low to high relative ecological risk within each estuary.

For categorical mapping and comparison, the continuous PSR ecological risk scores were further analyzed using the Getis-Ord Gi* hotspot method same as above. Unlike the resilience-loss hotspot analysis, a grid was classified as an ecological risk hotspot/coldspot only when it showed an FDR-corrected significant result in the same direction under at least two of the three neighborhood radii. This stricter criterion was adopted because the PSR score is a continuous risk index, and the statistic identifies spatial clustering directly from these continuous values. Requiring agreement across multiple neighborhood scales not only enabled us to identify more stable spatial clusters of relatively high and low ecological risk, but also reduced the dependence of the classification on any single neighborhood radius.

To comprehensively evaluate the robustness of the PSR framework, we first examined the sensitivity of the PSR index to different construction methods by comparing the Spearman rank correlation coefficients between the baseline approach used in this study and alternative combinations of normalization and weighting methods (Supplementary Figure S12A). Second, a leave-one-out analysis was conducted by sequentially removing each of the three PSR subsystems (Pressure, State, and Response) as well as each individual indicator, followed by recalculation of the PSR scores and their Spearman rank correlation coefficients with the baseline results (Supplementary Figure S12C). Finally, the spatial robustness of the PSR-based ecological risk hotspot/coldspot patterns was evaluated by comparing the proportions of hotspots/coldspots obtained using the original significance test and the FDR-corrected results under different neighborhood radii. We further quantified the spatial consistency of hotspot/coldspot patterns by calculating the Jaccard similarity coefficient between the results from each neighborhood radius and those obtained using the smallest neighborhood radius (r = 1) as the reference (Supplementary Figure S13). In addition, a spatial permutation test was performed to determine whether the observed numbers of hotspots and coldspots deviated significantly from those expected under random spatial configurations (Supplementary Figure S12B).

Therefore, potential risk areas were identified by overlaying resilience-loss hotspots with areas where the PSR ecological risk assessment indicated relatively low risk. Conceptually, these areas represent locations where the PSR-based risk framework did not detect clear ecosystem degradation or high current ecological risk, but where the detected ecological resilience change indicated significant critical-slowing-down signals. Such areas may still appear relatively stable according to conventional pressure–state–response indicators, yet their declining ecological resilience suggests that they may be approaching a critical transition or regime shift before substantial changes become evident in ecosystem state variables (). Because such regime shifts can be discontinuous, nonlinear, and characterized by time-lagged responses, and post-shift recovery may not simply return the ecosystem to its previous desirable state (Scheffer et al., 2001; Suding et al., 2004), we treated these overlap areas as potential risk zones requiring enhanced monitoring and management attention.

2.4 Explainable machine learning analysis

To identify the environmental and anthropogenic factors associated with resilience-loss hotspots, we fitted a binary random forest classification model using the resilience-loss hotspot product as the response variable. To reduce the influence of spatial autocorrelation during model evaluation, pixels were assigned to five spatial folds based on their geographic coordinates within each estuary, and model performance was evaluated using spatial-block cross-validation. Model performance was assessed using overall accuracy, hotspot sensitivity, other-class specificity, balanced accuracy, hotspot precision, and macro F1. Overall accuracy denotes the proportion of correctly classified pixels, hotspot sensitivity and other-class specificity quantify the proportions of correctly identified hotspot and non-hotspot pixels, respectively, balanced accuracy is the average of hotspot sensitivity and other-class specificity, hotspot precision represents the proportion of predicted hotspots that were truly hotspots, and macro F1 summarizes the balance between precision and recall across both classes.

We further calculated the Shapley additive explanations (SHAP) values via 50 Monte Carlo simulations with the adjustment option activated based on the trained random forest classification model to clarify the relative contributions of different environmental processes and anthropogenic activity factors to the classification of resilience-loss hotspots. SHAP values represent the marginal contribution of each predictor to the model prediction for every grid cell (; Nohara et al., 2019). In this study, positive SHAP values indicate that a predictor increases the probability of a grid cell being classified as a resilience-loss hotspot, whereas negative SHAP values indicate that it decreases this probability. To quantify the overall explanatory importance of individual predictors, we first calculated the mean absolute SHAP value of each predictor across all samples, which reflects the average magnitude of its contribution to the classification regardless of direction. The relative contribution of each predictor was then calculated by dividing its mean absolute SHAP value by the sum of the mean absolute SHAP values of all predictors. To further evaluate the contributions of different categories of environmental drivers, the mean absolute SHAP values of all predictors within the same driver category were summed and divided by the total mean absolute SHAP value across all predictors. This yielded the relative contributions of the major driver categories, including climate, hydrology, vegetation, topography, sea level, and human activities.

To further characterize the response patterns of key predictive variables, partial dependence plots (PDPs) were generated for the top-ranked predictors based on SHAP importance rankings (). For each selected predictor in the PDP analysis, individual conditional expectation curves were first estimated, followed by aggregation to obtain the average partial dependence response across all samples. The 99% confidence interval was derived from the variability in individual conditional expectation predictions. To facilitate interpretation, predictors exhibiting highly skewed distributions were visualized after excluding the outer 5% of their observed value range. This trimming was applied solely for visualization purposes to reduce the influence of sparsely sampled tails and to avoid overinterpreting response patterns supported by very few observations, and the underlying model predictions were not modified. In addition, because partial dependence plots may involve predictions for combinations of predictor values that are weakly represented or absent in the training data, portions of the curves may reflect model extrapolation rather than empirical observations (). Therefore, our interpretation focused primarily on the overall response direction and broad nonlinear trends within well-sampled regions of the predictor space, rather than on local fluctuations or responses at the distribution tails.

3 Results

3.1 Ecological resilience based on critical slowing down indicators

We calculated the trend of lag-1 autocorrelation (AC1) based on NDVI time series from 2000 to 2024 for each 250 m grid cell (Figure 2; Table 2), and defined areas with significant positive trend of AC1 as ecological resilience loss based on CSD. Our results showed that approximately 42.26% of the estuarine regions were losing resilience. Among them, the area of the Liaohe Estuary that were losing resilience was the largest, accounting for about 46.66%, while the area of the Pearl Estuary that were losing resilience was the smallest, accounting for about 38.99%. Spatially, the areas of resilience loss in the Liaohe Estuary were mainly concentrated on its southwestern side, whereas those in the Yangtze, Minjiang, and Pearl River estuaries showed more fragmented distributions.

FIGURE 2

TABLE 2

EstuaryNumber of valid pixelsTrend of AC1
Sig. PositiveNon-sig. TrendSig. Negative
Liaohe Estuary17,92046.66%11.22%42.12%
Yangtze Estuary5,73941.14%14.13%44.73%
Minjiang Estuary1,52942.25%15.63%42.12%
Pearl River Estuary2,86538.99%15.25%45.76%

Summary of AC1 trend categories across the four estuaries. Sig. represents significant.

We further used the Getis-Ord Gi* statistic based on the ecological resilience to identify the hotspots of resilience loss (Figure 3; Table 3). The results showed pronounced spatial heterogeneity among the four estuaries. The Liaohe Estuary exhibited the strongest signal of resilience loss, with 32.98% of valid pixels classified as hotspots at the 99% confidence level, indicating a highly significant spatial clustering of resilience loss. By contrast, most pixels in the Yangtze, Minjiang, and Pearl River estuaries were classified as not significant, accounting for 87.80%, 96.60%, and 87.96% of valid pixels, respectively. Spatially, resilience-loss hotspots in the Liaohe Estuary were mainly concentrated in the southwestern part of the wetland area, whereas those in the Yangtze and Minjiang estuaries were primarily distributed within the estuarine interior. In the Pearl River Estuary, hotspot pixels were relatively sparse and mainly located on the seaward side of the estuarine wetland.

FIGURE 3

TABLE 3

EstuaryColdspot 99%Coldspot 95%Coldspot 90%Not sigHotspot 90%Hotspot 95%Hotspot 99%
Liaohe Estuary32.98%5.08%3.21%18.83%2.33%5.01%32.56%
Yangtze Estuary0.82%2.27%2.49%87.80%3.00%3.00%0.63%
Minjiang Estuary0.00%0.46%0.72%96.60%0.78%0.72%0.72%
Pearl River Estuary0.10%1.47%3.18%87.96%3.73%2.23%1.33%

Percentage distribution of hotspot and coldspot classes of resilience loss at different confidence levels across the four estuaries. Sig. represents significant.

3.2 Relative ecological risk based on Press-State-Response indicators

We also developed an indicator system based on the Press-State-Response (PSR) framework to assess the relative ecological risk across the four estuarine wetlands. Due to the significant differences in the range of each estuary area, the sample distribution was largely uneven. Therefore, we conducted risk classification comparisons within each estuary. A hotspot analysis of ecological risk revealed spatial patterns that were broadly similar to those of resilience-loss hotspots (Figure 4; Table 4), although notable spatial mismatch existed between the two assessments. About 20.27% of the regions were classified as high ecological risk hotspots. In the Liaohe Estuary, areas with relatively high ecological risk were primarily concentrated in the seaward. In both the Yangtze River Estuary and the Minjiang Estuary, risk levels were generally higher in the south than in the north. In the Pearl River Estuary, high-risk areas were mainly distributed near the sea. When the classification criterion was tightened to require significance across all three neighborhood radii, the spatial extent of hotspots and coldspots showed no substantial difference (Supplementary Figure S14).

FIGURE 4

TABLE 4

EstuaryMean PSR scoreColdspot 99%Coldspot 95%Coldspot 90%Not sigHotspot 90%Hotspot 95%Hotspot 99%
Liaohe Estuary0.2739.41%16.17%4.56%17.07%0.41%1.12%21.32%
Yangtze Estuary0.269.21%4.96%4.25%65.17%3.18%4.16%9.07%
Minjiang Estuary0.359.83%6.59%4.76%56.29%3.97%3.79%14.77%
Pearl River Estuary0.341.98%7.73%5.64%65.36%2.02%4.65%12.62%

Percentage distribution of hotspot and coldspot classes of relative ecological risk assessed by indicators derived from PSR at different confidence levels across the four estuaries. Sig. represents significant.

Notably, some areas identified as hotspots of resilience loss were not classified as high-risk areas in the PSR-based ecological risk assessment (Figure 5; Table 5). We defined these regions as potential risk areas by overlaying the resilience-loss hotspots with the PSR ecological-risk coldspots, and pixels classified as resilience-loss hotspots and PSR coldspots at the corresponding confidence levels were retained. These areas therefore represent locations where the PSR framework still indicated relatively low ecological risk, but vegetation temporal dynamics had already shown increased AC1 and reduced ecological resilience, reflecting a possible pre-degradation state before clear degradation signals emerged in the PSR-based assessment. Our results showed that there remained pronounced spatial heterogeneity among the four estuaries. Approximately 60.13% of all valid pixels were classified as potential risk areas, accounting for 61.3% of the total resilience-loss hotspots in the Liaohe Estuary, while the proportions were relatively smaller in other estuaries. Spatially, these areas were mainly clustered in the central-southwestern part of the Liaohe Estuary, while they appeared as scattered patches in the other three estuaries.

FIGURE 5

TABLE 5

EstuaryRatio of resilience-loss hotspotsRatio of potential risk areasRatio of potential risk area in resilience-loss hotspots
Liaohe Estuary39.89%60.13%61.30%
Yangtze Estuary6.62%18.42%15.00%
Minjiang Estuary2.22%21.18%17.65%
Pearl River Estuary7.29%15.36%16.27%

Proportions of resilience-loss hotspots, ecological risk areas, and ecological risk areas within resilience-loss hotspots across the four estuaries.

3.3 Impacts of environmental factors on resilience change

We fitted a random forest classification model to identify the factors associated with resilience-loss hotspots. Under spatial-block cross-validation, the model achieved an overall accuracy of 72.41% (Supplementary Figure S15). SHAP analysis was then used to identify the dominant environmental and anthropogenic drivers contributing to hotspot classification (Figure 6). Among all the predictor variables, the sea level variance contributed the most (21%), followed by the sea level mean (19%) and precipitation mean (13%). To analyze the key processes related to resilience-loss hotspots, we classified the explanatory variables. Among them, climate variables overall explained the largest proportion of the variation (44.22%), followed by sea level (39.75%), and human activity (6.82%).

FIGURE 6

Partial dependence analysis summarized how the fitted random forest model classified resilience-loss hotspots in relation to the main predictors (Figure 7). The five strongest predictors showed clear but mixed associations with hotspot probability. In general, sea level variance and mean sea level were negatively associated with predicted hotspot probability, although the strongest positive SHAP contributions were concentrated at the lower end of their observed ranges (Figures 7A,B). Mean precipitation showed a similar pattern, with higher hotspot probability under relatively low precipitation conditions (Figure 7C). Temperature variance was positively associated with hotspot probability within the main range of the data distribution, suggesting that areas with higher temperature variability were more likely to be classified as resilience-loss hotspots (Figure 7D). Mean temperature showed approximately negative responses, with positive contributions mainly at low values and weaker or negative contributions at higher values (Figure 7E). Overall, these patterns in model -response suggest that sea-level conditions and climatic variability were key factors distinguishing resilience-loss hotspots from other areas.

FIGURE 7

4 Discussion

Based on resilience theory and the hypothesis of alternative stable-state transitions, we detected CSD signals of resilience loss across tidal flats, salt marshes, and mangroves in four major estuarine wetlands of eastern China. The three habitat types were treated as components of a single estuarine wetland system here, which was based on the assumption that estuarine wetlands can exhibit alternative stable states and that transitions among the three habitats may occur naturally under changing environmental conditions (Wang and Temmerman, 2013; Saintilan et al., 2014; Moffett et al., 2015). From an ecosystem-level perspective, such transitions represent state shifts within the estuarine wetland system rather than conversions between independent ecosystem types (Moffett et al., 2015). Previous studies have shown that resilience loss signals are robust to different ecosystem types, and many studies assessing resilience loss combined different ecosystem types within biomes (e.g., forests, drylands, and wetlands, Smith et al., 2022; ; ). We also conducted the hotspot analysis separately for each habitat type. The results showed that most of the pixels remain stable during the study period (Supplementary Figure S4), and no pixels exhibited opposite hotspot/coldspot classifications between the two approaches. Observed inconsistencies mainly resulted from differences in statistical significance (Supplementary Figure S8). Therefore, we evaluated resilience at the ecosystem level instead of separately for each habitat type. It should be noted that this simplification may mask habitat-specific resilience dynamics and responses to environmental drivers. In this case, our results should be interpreted as reflecting the overall resilience of estuarine wetland ecosystems rather than the resilience of individual wetland types.

In addition, it should be noted that in some cases, CSD signals like trend of AC1 might fail to capture the approaching shift (). Recent studies have further emphasized that CSD indicators should not be regarded as universal predictors of regime shifts, because their performance can be affected by non-stationary environmental forcing, changing noise characteristics, limited time-series length, and scale mismatches between observations and underlying ecological processes (; ). Accordingly, recent approaches increasingly recommend combining CSD indicators with complementary evidence, such as direct recovery metrics, spatial resilience patterns, and process-based information, particularly in spatially heterogeneous coastal and deltaic systems (; Wang Y. et al., 2025). However, our study focuses mainly on the resilience change detected by historical vegetation dynamics, rather than regime shifts had already occurred. Our robustness test showed that there were no significant differences in trend of AC1 between pixels with and without regime shifts (Supplementary Figure S3). Consequently, our results should help identify ecosystems that gradually approach tipping points and should be interpreted primarily as indicators of resilience loss rather than deterministic predictors of regime shifts.

Furthermore, it has been evidenced that estuarine wetlands can exhibit alternative stable states and tipping points arising from coupled feedbacks among vegetation, hydrological processes, and geomorphological dynamics (Moffett et al., 2015; ). Empirical observations and transplantation experiments in tidal marshes have demonstrated progressively slower vegetation recovery with increasing inundation stress, providing direct evidence of critical slowing down in real-world wetlands (van Belzen et al., 2017). However, it is always difficult to reliably identify perturbation events across a large scale purely from remote sensing data, and there may be some inconsistency between the change of AC1 trend and recovery rate (). Therefore, we chose to focus on the ability of the ecosystem to maintain its stable state in the face of incoming disturbance, which could be characterized by the AC1 trend used here. Moreover, theoretical analyses of salt-marsh dynamics have shown that temporal autocorrelation would increase as the system approaches a critical transition (Neijnens et al., 2021).

Importantly, the results in line with our expectation persisted across multiple scenarios of sensitive analysis, which were largely consistent with previous studies (Figure 2). For instance, in the western part of the Liaohe Estuary, some wetlands have been recorded as having been converted into farmland or impervious surfaces (). The particularly high proportion of resilience-loss hotspots in the Liaohe Estuary may be partly explained by its relatively large and continuous wetland distribution, which could facilitate the formation of spatially clustered patterns compared with more fragmented wetland landscapes in other estuaries. The Yangtze Estuary wetland experienced the dual disturbance of reclamation and biological invasion, which led to the decline of blue carbon storage (Wang T. et al., 2026). Areas losing resilience in the Minjiang Estuary showed biological invasion (), while the Pearl River Estuary faces a sudden decrease in sediment supply and habitat fragmentation (Zhao et al., 2010; Yuan et al., 2025).

We further overlaid the resilience-loss hotspot/coldspot with protected areas and compared the environmental characteristics of hotspot and coldspot areas. Our results showed that coldspots were not consistently concentrated within protected areas (Supplementary Figure S7; ). In addition, some protected-area pixels were still identified as resilience-loss hotspots, indicating that vegetation resilience within protected areas was not necessarily better than that outside protected areas. This pattern suggests that formal protection boundaries alone do not always correspond to areas with enhanced resilience, and that continued monitoring of vegetation resilience within protected areas remains necessary (). In addition, there exist significant differences between hotspot and coldspot areas in several environmental factors (Supplementary Figure S10), which might suggest association between certain environmental conditions and spatial variance in vegetation resilience. However, the direction of these differences varied among estuaries, indicating that vegetation resilience in estuarine wetlands emerges from complex interactions among hydrological dynamics, sediment processes, climatic variability, and human disturbances (; Moffett et al., 2015).

Merely considering the dynamic characteristics of the ecosystem resilience is a one-sided view. It is also important to assess the current ecological risk to comprehensively assess the condition of estuarine wetlands. Therefore, we introduced the PSR framework to evaluate the current risk of four major estuarine wetlands. In order to use a simple, easily obtained indicator system, we carefully considered the pressure, state, and response factors to be included in the framework, and selected as much remote sensing data as possible. For example, to comprehensively reflect the environmental pressure over the past 24 years, we not only included the mean value of temperature, precipitation, and sea level from 2000 to 2024, but also calculated their variances to reflect the fluctuation features, which might represent the intensity of disturbance (). We introduced the Modified Normalized Difference Water Index (MNDWI) as a proxy for hydrological conditions. The 10th percentile of MNDWI (P10) was used as a state indicator to characterize the baseline hydrological condition of the wetland. Inundation duration was calculated based on the periods when MNDWI >0 and was included in the PSR pressure subsystem because prolonged inundation represents hydrological stress on wetland vegetation. We also conducted a comparison to test if the MNDWI>0 criterion was suitable for all four estuaries to capture flooding (Supplementary Figure S9). Similar threshold-based approaches have been widely applied in large-scale wetland and surface water mapping studies because they provide a consistent method for characterizing hydrological conditions across broad regions (Xu, 2006; Pekel et al., 2016). Similarly, mean NDVI was considered in the PSR state subsystem to characterize the vegetation state, while robust slope and variance calculated with NDVI residual were considered as response indicators to reflect the ability of wetland vegetation to maintain stability under environmental stress. It also should be noted that we normalized the indicators with a linear method. Although the invasion of S. alterniflora, which has been one of the most important factors influencing estuarine wetlands, would increase NDVI, it might also protect the shoreline and promote carbon storage, thus mitigate ecological risks (Wan et al., 2009; Zhang et al., 2023). Potential redundancy among these indicators was evaluated through leave-one-out sensitivity test (Supplementary Figure S12C), and the results showed that it would be unlikely to introduce substantial differences for choosing indicators from the same source.

A considerable proportion of estuarine wetland regions has been identified as potential risk areas (Figure 5; Table 5). The interpretation of these areas should be considered within the complementary framework of resilience assessment and PSR-based ecological risk evaluation. The PSR framework characterizes current ecological risk by integrating information on environmental pressures, ecosystem states, and responses, whereas resilience indicators provide insights into changes in ecosystem stability. Therefore, PSR-based ecological risk coldspots do not necessarily indicate the absence of ecological stress or future vulnerability. Rather, they represent areas where the current combination of pressure, state, and response indicators has not resulted in relatively high ecological risk. The overlap between resilience-loss hotspots and ecological risk coldspots thus highlights a mismatch between current ecological risk status and ecosystem stability dynamics. These areas may maintain relatively favorable observable conditions while simultaneously experiencing resilience loss, indicating a potential decline in their capacity to absorb future disturbances. Such divergence is ecologically important because ecosystem transitions may involve nonlinear processes and time lags, and changes in ecosystem stability can occur before substantial alterations become evident in conventional state indicators (Scheffer et al., 2001; ; Törnqvist et al., 2020).

Importantly, potential risk areas should not be interpreted as confirmed degradation zones or deterministic predictions of future regime shifts. Instead, they represent priority areas for further monitoring and adaptive management. The robustness of this classification was evaluated through sensitivity analyses involving alternative index construction methods, normalization procedures, weighting schemes, and spatial permutation tests (Supplementary Figure S11, S12). Despite uncertainties associated with the complexity of estuarine ecosystems and the availability of selectable indicators, the identified patterns remained generally consistent across different analytical settings. Therefore, potential risk areas provide a useful screening framework for identifying locations where conventional risk indicators may not yet indicate pronounced degradation, while vegetation dynamics suggest declining ecosystem resilience. Early intervention in such areas may provide greater opportunities for effective conservation, and reduce financial and ecological costs associated with restoration after severe degradation has occurred ().

We focused on resilience-loss hotspots rather than continuous AC1 trends because identifying spatially clustered regions of resilience decline provides more meaningful information for ecological resilience assessment and management. Although pixel-level AC1 trends capture temporal changes in ecological resilience, they can be sensitive to local time-series quality, disturbances, and short-term environmental variability, resulting in scattered significant pixels that may have limited operational value for restoration. By incorporating spatial context and statistical significance, hotspot analysis identifies coherent clusters of increasing AC1, highlighting areas where resilience loss is concentrated. These spatially connected regions represent more meaningful management units, where early intervention, restoration efforts, and conservation resources can be more effectively prioritized. We also tested alternative response variables, including continuous AC1 trend and significantly increasing AC1 class (Supplementary Figure S16, S17). However, these models showed weak predictive performance. This may be because pixel-level AC1 trend is influenced by multiple sources of variability that are not fully represented by the available environmental covariates, including local disturbance events, unmeasured management actions, vegetation community composition, and residual noise in the NDVI time series (; Pettorelli et al., 2014). In addition, the relationship between environmental drivers and AC1 trend may be spatially nonstationary and highly nonlinear, making continuous pixel-level predictions difficult (Wiens, 1989; ). The significantly increasing AC1 class response may further amplify uncertainty because it combines trend magnitude and significance thresholds at the pixel level without accounting for spatial clustering. Therefore, the hotspot-based response provided a more interpretable and spatially robust target for linking resilience-loss patterns to environmental and anthropogenic drivers.

The variability and mean value of sea level and precipitation, together with temperature variability, were identified as the most important predictors for classifying resilience-loss hotspots (Figure 6). The observed response patterns suggest that resilience-loss hotspots are primarily associated with variations in climate and sea-level conditions. Nevertheless, the Partial Dependence Plots (PDPs) are intended to illustrate the model’s learned associations rather than provide direct estimates of ecological responses (; ). Owing to the relatively small spatial extent of the estuarine wetlands and the limited distribution of the remaining natural coastal wetlands in China, several explanatory variables exhibit inherently restricted variation, resulting in substantial differences in sample density across predictor ranges. Consequently, some local oscillations and non-monotonic changes in the PDP curves may reflect the underlying data structure and model fitting rather than robust ecological responses. Therefore, our interpretation focused on the overall response patterns within well-sampled regions rather than localized fluctuations or extrapolated portions of the curves. Despite these limitations, the consistently nonlinear response patterns observed across multiple predictors suggest that resilience responses to environmental change are unlikely to be governed by simple linear relationships, but instead are shaped by ecological thresholds, feedback mechanisms, and interactions among multiple environmental drivers (; Scheffer et al., 2001).

Estuarine wetlands occupy the land-sea interface and are simultaneously influenced by both marine and freshwater processes (). In our model, mean sea level and sea-level variability were among the strongest predictors associated with resilience-loss hotspot classification, highlighting that broad-scale hydrogeomorphic settings provide important information for distinguishing areas vulnerable to resilience loss (). By contrast, inundation duration showed a relatively weak contribution despite representing local hydrological conditions. One possible explanation is that estuarian wetland vegetation has adapted to regular tidal inundation through long-term physiological and morphological adjustments, such that routine flooding duration alone provides limited discriminatory information within the observed environmental gradients (Silvestri et al., 2005; Munns and Tester, 2008). Given that mean sea level may better characterize the background hydrological stress experienced by vegetation, whereas and that sea-level variability may capture fluctuations associated with episodic flooding or other transient hydrodynamic conditions (Morris et al., 2002), these two variables are more strongly associated with the spatial distribution of resilience-loss hotspots than average inundation duration.

Climate-related variables also showed strong associations with resilience-loss hotspot classification. Mean temperature and precipitation characterize the long-term climatic setting under which estuarine vegetation develops, whereas their temporal variability reflects the frequency and magnitude of climatic fluctuations (Woodward, 1987; Smith, 2011). Compared with long-term averages, climatic variability may better distinguish locations exposed to more dynamic environmental conditions, and have been widely recognized as an important characteristic of estuarine ecosystems (). Importantly, at the land-sea interface, climatic variability may influence estuarine vegetation resilience through both direct and indirect pathways. In addition to directly affecting vegetation physiological processes, precipitation variability can alter freshwater discharge and subsequently modify estuarine salinity regimes, sediment delivery, and inundation conditions, while temperature variability may interact with these hydrological processes to further influence vegetation growth and recovery (). Such coupled climate-hydrology-sediment interactions can amplify or buffer climatic effects on vegetation, depending on local environmental conditions (). Therefore, the high importance of precipitation and temperature variability in the classification model suggests that resilience-loss hotspots are statistically associated with areas experiencing greater hydroclimatic variability.

These interpretations should nevertheless be cautioned with a few limitations of the modeling framework in mind. Environmental predictors in coastal wetlands are inherently spatially autocorrelated and often correlated with one another because they are governed by common climatic and hydrological gradients (; ). In addition, the relatively limited spatial extent of the remaining natural coastal wetlands in China results in uneven predictor distributions and restricted value ranges for several variables. Consequently, the variable importance and partial dependence patterns identified by the random forest model should be interpreted as model-derived associations that summarize the combined effects of correlated environmental gradients, rather than independent causal relationships. Despite these limitations, the consistent dominance of sea-level and hydroclimatic variables across both SHAP importance and PDP analyses suggests that these environmental gradients provide robust explanatory information for distinguishing resilience-loss hotspots.

5 Conclusion

Based on the lag-1 autocorrelation trend of the NDVI time series, we evaluated the spatial patterns of ecological resilience in four typical estuaries in eastern China. The results showed that approximately 42.26% of the estuarine areas were losing resilience, among which the Liaohe Estuary had the largest area of resilience loss (46.66%), and the Pearl Estuary had the lowest (38.99%). Hotspot analysis showed substantial spatial inconsistency in resilience-loss hotspots among these four estuaries. The PSR-based ecological risk assessment framework indicated that approximately 23.46% of the study area were classified as relatively high ecological risk hotspots. Potential risk areas, defined as the overlap between resilience-loss hotspots and ecological-risk coldspots, also showed pronounced differences among the four estuaries. The random forest classification model showed moderate predictive performance under spatial-block cross-validation, with an overall accuracy of 72.05%. SHAP analysis indicated that sea-level variability, mean sea level, precipitation conditions, and temperature variability were the main factors associated with resilience-loss hotspot classification. These predictors showed complex model-response patterns, overall suggesting that estuarine wetland resilience is jointly shaped by sea-level conditions and climatic variability.

This study also has several methodological limitations that should be addressed in future work. First, the analysis was based primarily on remotely sensed indicators and gridded environmental covariates, which may not fully capture local ecological processes, field-scale vegetation dynamics, or management interventions. Second, although spatial cross-validation was used to reduce spatial overfitting, the correlative modeling framework cannot by itself establish causal relationships between environmental drivers and resilience loss. Third, differences in spatial and temporal resolution among MODIS, VIIRS, ERA5-Land, CMEMS, and terrain datasets may introduce scale-related uncertainty despite harmonized extraction procedures. Future research should integrate field observations, higher-resolution remote sensing, and process-based ecological or hydrodynamic models to better validate resilience indicators, clarify causal mechanisms, and assess how coastal wetland resilience may respond to future climate change, sea-level rise, and human disturbance. In addition, future efforts could link resilience-loss detection with enhanced monitoring, disturbance control, and adaptive intervention to facilitate the transition from resilience assessment toward operational early-warning and management systems.

Overall, our research identified hotspots of resilience loss and potential risk areas in key estuarine regions, and analyzed the main factors associated with resilience-loss hotspots classification. In the face of the dual pressures of climate change and human activities, our findings serve as a key step toward dynamic assessment frameworks for early degradation warning and differentiated ecological management of estuarine wetlands.

Statements

Data availability statement

The original contributions presented in the study are included in the article/Supplementary Material, further inquiries can be directed to the corresponding authors.

Author contributions

QD: Writing – original draft, Formal Analysis. YL: Writing – original draft, Validation, Data curation. WH: Data curation, Software, Writing – original draft, Formal Analysis. BA: Methodology, Writing – review and editing, Conceptualization. ZZ: Visualization, Validation, Writing – original draft. CX: Conceptualization, Writing – review and editing, Project administration. ST: Writing – review and editing, Methodology.

Funding

The author(s) declared that financial support was received for this work and/or its publication. This work was funded by the National Key R&D Program of China (2022YFF1301005).

Conflict of interest

The author(s) declared that this work 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) declared that generative AI was not 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.

Supplementary material

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fenvs.2026.1912598/full#supplementary-material

References

  • 1

    ApleyD. W.ZhuJ. (2020). Visualizing the effects of predictor variables in black box supervised learning models. J. R. Stat. Soc. Ser. B Stat. Methodol.82 (4), 1059–1086. 10.1111/rssb.12377

  • 2

    BalmfordA.GastonK. J.BlythS.JamesA.KaposV. (2003). Global variation in terrestrial conservation costs, conservation benefits, and unmet conservation needs. Proc. Natl. Acad. Sci.100 (3), 1046–1050. 10.1073/pnas.0236945100

  • 3

    BarendregtA.SwarthC. W. (2013). Tidal freshwater wetlands: variation and changes. Estuaries Coasts36 (3), 445–456. 10.1007/s12237-013-9626-z

  • 4

    BathianyS.BastiaansenR.BastosA.BlaschkeL.LeverJ.LorianiS.et al (2025). Ecosystem resilience monitoring and early warning using Earth observation data: challenges and outlook. Surv. Geophys.46 (2), 265–301. 10.1007/s10712-024-09833-z

  • 5

    BenjaminiY.HochbergY. (1995). Controlling the false discovery rate: a practical and powerful approach to multiple testing. J. R. Statistical Society Series B Methodol.57 (1), 289–300. 10.1111/j.2517-6161.1995.tb02031.x

  • 6

    BiggsR.CarpenterS. R.BrockW. A. (2009). Turning back from the brink: detecting an impending regime shift in time to avert it. Proc. Natl. Academy Sci.106 (3), 826–831. 10.1073/pnas.0811729106

  • 7

    BlainG. C. (2013). The Mann-Kendall test: the need to consider the interaction between serial correlation and trend. Acta Sci. Agron.35, 393–402. 10.4025/actasciagron.v35i4.16006

  • 8

    BoultonC. A.LentonT. M.BoersN. (2022). Pronounced loss of amazon rainforest resilience since the early 2000s. Nat. Clim. Change12 (3), 271–278. 10.1038/s41558-022-01287-8

  • 9

    ChenX.ZhangM.ZhangW. (2022). Landscape pattern changes and its drivers inferred from salt marsh plant variations in the coastal wetlands of the Liao river Estuary, China. Ecol. Indic.145, 109719. 10.1016/j.ecolind.2022.109719

  • 10

    ChenZ. L.YiY.ZhangH.LiP.WangY.YanZ.et al (2024). Differences in dissolved organic matter molecular composition along two plume trajectories from the Yangtze river estuary to the East China sea. ACS Environ. Au4 (1), 31–41. 10.1021/acsenvironau.3c00030

  • 11

    ChengY.ZhaY.TongC.HuM.DuD.FanY.et al (2020). Plant population dynamics in a degraded coastal wetland and implications for the carbon cycle. Wetlands40 (5), 1617–1625. 10.1007/s13157-020-01268-7

  • 12

    China Nature Reserve Specimen Resource Sharing Platform (2024). List and Vector Boundaries of Nature Reserves in China (Version v1).

  • 13

    CloernJ. E.FosterS.KlecknerA. (2014). Phytoplankton primary production in the world’s estuarine-coastal ecosystems. Biogeosciences11 (9), 2477–2501. 10.5194/bg-11-2477-2014

  • 14

    Robertson, H. A. (2025). “Global wetland outlook 2025: valuing, conserving, restoring and financing wetlands. Ramsar Convention on Wetlands. 10.69556/GWO-2025-eng

  • 15

    Group on Environmental Performance and Development. Group on the State of the Environment (1993). OECD core set of indicators for environmental performance reviews (No. 83). Organisation for Economic Co-operation and Development.

  • 16

    CruzM.RosmanP. C. C.AlmeidaC. M. R.IglesiasI. (2026). Modelling estuarine wetlands and their role as nature-based solutions: a review. Sci. Total Environ.1014, 181325. 10.1016/j.scitotenv.2025.181325

  • 17

    DakosV.KéfiS. (2022). Ecological resilience: what to measure and how. Environ. Res. Lett.17 (4), 043003. 10.1088/1748-9326/ac5767

  • 18

    DakosV.CarpenterS. R.BrockW. A.EllisonA. M.GuttalV.IvesA. R.et al (2012a). Methods for detecting early warnings of critical transitions in time series illustrated using simulated ecological data. PloS One7 (7), e41010. 10.1371/journal.pone.0041010

  • 19

    DakosV.van NesE. H.D’OdoricoP.SchefferM. (2012b). Robustness of variance and autocorrelation as indicators of critical slowing down. Ecology93 (2), 264–271. 10.1890/11-0889.1

  • 20

    DakosV.CarpenterS. R.van NesE. H.SchefferM. (2015). Resilience indicators: prospects and limitations for early warnings of regime shifts. Philosophical Trans. R. Soc. B Biol. Sci.370 (1659), 20130263. 10.1098/rstb.2013.0263

  • 21

    DayJ.AnthonyE.CostanzaR.EdmondsD.GunnJ.HopkinsonC.et al (2024). Coastal wetlands in the anthropocene. Annu. Rev. Environ. Resour.49, 105–135. 10.1146/annurev-environ-121922-041109

  • 22

    DidanK. (2021). MOD13Q1 MODIS/Terra vegetation indices 16-Day L3 global 250m SIN grid V061. NASA EOSDIS Land Process. Distrib. Act. Arch. Cent. (LP DAAC). 10.5067/MODIS/MOD13Q1.061

  • 23

    DormannC. F.ElithJ.BacherS.BuchmannC.CarlG.CarréG.et al (2013). Collinearity: a review of methods to deal with it and a simulation study evaluating their performance. Ecography36 (1), 27–46. 10.1111/j.1600-0587.2012.07348.x

  • 24

    Elsey-QuirkT.GrahamS. A.MendelssohnI. A.SneddenG.DayJ. W.TwilleyR. R.et al (2019). Mississippi river sediment diversions and coastal wetland sustainability: synthesis of responses to freshwater, sediment, and nutrient inputs. Estuar. Coast. Shelf Sci.221, 170–183. 10.1016/j.ecss.2019.03.002

  • 25

    ElvidgeC. D.ZhizhinM.GhoshT.HsuF.-C.TanejaJ. (2021). Annual time series of global VIIRS nighttime lights derived from monthly averages: 2012 to 2019. Remote Sens.13 (5), 922. 10.3390/rs13050922

  • 26

    FivashG. S.TemmermanS.KleinhansM. G.HeunerM.van der HeideT.BoumaT. J. (2023). Early indicators of tidal ecosystem shifts in estuaries. Nat. Communications14 (1), 1911. 10.1038/s41467-023-37444-6

  • 27

    ForzieriG.DakosV.McDowellN. G.RamdaneA.CescattiA. (2022). Emerging signals of declining forest resilience under climate change. Nature608 (7923), 534–539. 10.1038/s41586-022-04959-9

  • 28

    FriedmanJ. H. (2001). Greedy function approximation: a gradient boosting machine. Ann. Statistics29 (5), 1189–1232. 10.1214/aos/1013203451

  • 29

    GeldmannJ.ManicaA.BurgessN. D.CoadL.BalmfordA. (2019). A global-level assessment of the effectiveness of protected areas at resisting anthropogenic pressures. Proc. Natl. Acad. Sci.116 (46), 23209–23215. 10.1073/pnas.1908221116

  • 30

    GoldsteinA.KapelnerA.BleichJ.PitkinE. (2015). Peeking inside the black box: visualizing statistical learning with plots of individual conditional expectation. J. Comput. Graph. Statistics24 (1), 44–65. 10.1080/10618600.2014.907095

  • 31

    GuJ.LuoM.ZhangX.ChristakosG.AgustiS.DuarteC. M.et al (2018). Losses of salt marsh in China: trends, threats and management. Estuar. Coast. Shelf Sci.214, 98–109. 10.1016/j.ecss.2018.09.015

  • 32

    GuisanA.ZimmermannN. E. (2000). Predictive habitat distribution models in ecology. Ecol. Model.135 (2-3), 147–186. 10.1016/s0304-3800(00)00354-9

  • 33

    HanY.LinQ.HuangS.DuC.ShenJ.ZhangK. (2024). Human impacts dominate global loss of lake ecosystem resilience. Geophys. Res. Lett.51 (11), e2024GL109298. 10.1029/2024gl109298

  • 34

    HeJ.WangC.HanY.ZhangW.YangY.GuoL.et al (2025). Waterbird community response to wetland and climate changes in the Liaohe river Estuary wetlands, China. J. Environ. Manag.374, 124165. 10.1016/j.jenvman.2025.124165

  • 35

    HersbachH.BellB.BerrisfordP.BiavatiG.HorányiA.Muñoz SabaterJ.et al (2023). ERA5 monthly averaged data on pressure levels from 1940 to present, Copernicus Climate Change Service (C3S) Climate Data Store (CDS).

  • 36

    IbáñezC.CaiolaN.ElphickC. S.GrillasP. (2026). To what extent is coastal wetland biodiversity endangered by climate change? How can we boost the resilience of coastal ecosystems?Nat. Conserv.62, 63–94. 10.3897/natureconservation.62.154795

  • 37

    JaccardP. (1901). Étude comparative de la distribution florale dans une portion des Alpes et des Jura. Bull. la Société vaudoise Sci. Nat.37, 547–579. 10.5169/seals-266450

  • 38

    Jean-MichelL.EricG.RomainB.-B.GillesG.AngéliqueM.MarieD.et al (2021). The Copernicus global 1/12 oceanic and sea ice GLORYS12 reanalysis. Front. Earth Sci.9, 698876. 10.3389/feart.2021.698876

  • 39

    JiB.LiH.ZhangH.LiW.MengQ. (2025). Carbon fluxes in estuarine wetlands based on a material flow perspective. Cont. Shelf Res.295, 105571. 10.1016/j.csr.2025.105571

  • 40

    KennishM. J. (2002). Environmental threats and environmental future of estuaries. Environ. Conserv.29 (1), 78–107. 10.1017/s0376892902000061

  • 41

    KhojastehD.RaoS.McSweeneyS.IbacetaR.NichollsR. J.FrenchJ.et al (2025). Intermittent estuaries deserve global attention as vulnerable and vital ecosystems. Commun. Earth and Environ.6 (1), 443. 10.1038/s43247-025-02428-5

  • 42

    KirwanM. L.MegonigalJ. P. (2013). Tidal wetland stability in the face of human impacts and sea-level rise. Nature504 (7478), 53–60. 10.1038/nature12856

  • 43

    KongZ.LingH.DengM.LongA.GuoB.HuangS.et al (2026). Intensifying drought risk weakens vegetation resilience in global drylands. Commun. Earth and Environ.7 (1), 279. 10.1038/s43247-026-03303-7

  • 44

    LegendreP. (1993). Spatial autocorrelation: trouble or new paradigm?Ecology74 (6), 1659–1673. 10.2307/1939924

  • 45

    LentonT. M. (2020). Tipping positive change. Philosophical Trans. R. Soc. B375 (1794), 20190123. 10.1098/rstb.2019.0123

  • 46

    LentonT. M.AbramsJ. F.BartschA.BathianyS.BoultonC. A.BuxtonJ. E.et al (2024). Remotely sensing potential climate change tipping points across scales. Nat. Commun.15 (1), 343. 10.1038/s41467-023-44609-w

  • 47

    LevinS. A. (1992). The problem of pattern and scale in ecology: the Robert H. MacArthur award lecture. Ecology73 (6), 1943–1967. 10.2307/1941447

  • 48

    LiC.FangS.GengX.YuanY.ZhengX.ZhangD.et al (2023). Coastal ecosystem service in response to past and future land use and land cover change dynamics in the Yangtze river estuary. J. Clean. Prod.385, 135601. 10.1016/j.jclepro.2022.135601

  • 49

    LiuL.LinB.FangQ.JiangX. (2024). Effectiveness assessment of China’s coastal wetland ecological restoration: a meta-analysis. Sci. Total Environ.934, 173336. 10.1016/j.scitotenv.2024.173336

  • 50

    LundbergS. M.LeeS.-I. (2017). A unified approach to interpreting model predictions. Adv. Neural Inf. Process. Syst.30, 4765–4774.

  • 51

    MitchellS. B.JennerjahnT. C.VizziniS.ZhangW. (2015). Changes to processes in estuaries and coastal waters due to intense multiple pressures–An introduction and synthesis. Estuar. Coast. Shelf Sci.156, 1–6. 10.1016/j.ecss.2014.12.027

  • 52

    MoffettK. B.NardinW.SilvestriS.WangC.TemmermanS. (2015). Multiple stable states and catastrophic shifts in coastal wetlands: progress, challenges, and opportunities in validating theory using remote sensing and other methods. Remote Sens.7 (8), 10184–10226. 10.3390/rs70810184

  • 53

    MorrisJ. T.SundareshwarP.NietchC. T.KjerfveB.CahoonD. R. (2002). Responses of coastal wetlands to rising sea level. Ecology83 (10), 2869–2877. 10.1890/0012-9658(2002)083[2869:rocwtr]2.0.co;2

  • 54

    MunnsR.TesterM. (2008). Mechanisms of salinity tolerance. Annu. Rev. Plant Biol.59, 651–681. 10.1146/annurev.arplant.59.032607.092911

  • 55

    NeijnensF. K.SiteurK.van de KoppelJ.RietkerkM. (2021). Early warning signals for rate-induced critical transitions in salt marsh ecosystems. Ecosystems24 (8), 1825–1836. 10.1007/s10021-021-00610-2

  • 56

    NoharaY.MatsumotoK.SoejimaH.NakashimaN. (2019). “Explanation of machine learning models using improved shapley additive explanation,” in Proceedings of the 10th ACM international conference on bioinformatics, computational biology and health informatics.

  • 57

    PeinJ.StanevaJ. (2026). Effectiveness of estuarine adaptation strategies under future climate conditions. Sci. Rep.16 (1), 9655. 10.1038/s41598-026-43040-7

  • 58

    PekelJ.-F.CottamA.GorelickN.BelwardA. S. (2016). High-resolution mapping of global surface water and its long-term changes. Nature540 (7633), 418–422. 10.1038/nature20584

  • 59

    PettorelliN.LauranceW. F.O’BrienT. G.WegmannM.NagendraH.TurnerW. (2014). Satellite remote sensing for applied ecologists: opportunities and challenges. J. Appl. Ecol.51 (4), 839–848. 10.1111/1365-2664.12261

  • 60

    PichonB.DonnetS.GounandI.AnandM.KéfiS. (2026). A framework for integrating resilience indicators into risk assessments. Conserv. Biol., e70379. 10.1111/cobi.70379

  • 61

    ReuterH. I.NelsonA.JarvisA. (2007). An evaluation of void‐filling interpolation methods for SRTM data. Int. J. Geogr. Inf. Sci.21 (9), 983–1008. 10.1080/13658810601169899

  • 62

    Rodrigues-FilhoJ. L.MacêdoR. L.SarmentoH.PimentaV. R. A.AlonsoC.TeixeiraC. R.et al (2023). From ecological functions to ecosystem services: linking coastal lagoons biodiversity with human well-being. Hydrobiologia850 (12), 2611–2653. 10.1007/s10750-023-05171-0

  • 63

    RotbarthR.Van NesE. H.SchefferM.JepsenJ. U.VindstadO. P. L.XuC.et al (2023). Northern expansion is not compensating for southern declines in North American boreal forests. Nat. Commun.14 (1), 3373. 10.1038/s41467-023-39092-2

  • 64

    SaintilanN.WilsonN. C.RogersK.RajkaranA.KraussK. W. (2014). Mangrove expansion and salt marsh decline at mangrove poleward limits. Glob. Change Biol.20 (1), 147–157. 10.1111/gcb.12341

  • 65

    SchefferM.CarpenterS.FoleyJ. A.FolkeC.WalkerB. (2001). Catastrophic shifts in ecosystems. Nature413 (6856), 591–596. 10.1038/35098000

  • 66

    SchefferM.BascompteJ.BrockW. A.BrovkinV.CarpenterS. R.DakosV.et al (2009). Early-warning signals for critical transitions. Nature461 (7260), 53–59. 10.1038/nature08227

  • 67

    SchefferM.CarpenterS. R.LentonT. M.BascompteJ.BrockW.DakosV.et al (2012). Anticipating critical transitions. Science338 (6105), 344–348. 10.1126/science.1225244

  • 68

    SilvestriS.DefinaA.MaraniM. (2005). Tidal regime, salinity and salt marsh plant zonation. Estuar. Coast. Shelf Sci.62 (1-2), 119–130. 10.1016/j.ecss.2004.08.010

  • 69

    SmithM. D. (2011). An ecological perspective on extreme climatic events: a synthetic definition and framework to guide future research. J. Ecol.99 (3), 656–663. 10.1111/j.1365-2745.2011.01798.x

  • 70

    SmithT.BoersN. (2023). Reliability of vegetation resilience estimates depends on biomass density. Nat. Ecol. and Evol.7 (11), 1799–1808. 10.1038/s41559-023-02194-7

  • 71

    SmithT.TraxlD.BoersN. (2022). Empirical evidence for recent global shifts in vegetation resilience. Nat. Clim. Change12 (5), 477–484. 10.1038/s41558-022-01352-2

  • 72

    SudingK. N.GrossK. L.HousemanG. R. (2004). Alternative states and positive feedbacks in restoration ecology. Trends Ecol. and Evol.19 (1), 46–53. 10.1016/j.tree.2003.10.005

  • 73

    TörnqvistT. E.JankowskiK. L.LiY.-X.GonzálezJ. L. (2020). Tipping points of Mississippi Delta marshes due to accelerated sea-level rise. Sci. Adv.6 (21), eaaz5512. 10.1126/sciadv.aaz5512

  • 74

    TuckerC. J. (1979). Red and photographic infrared linear combinations for monitoring vegetation. Remote Sens. Environ.8 (2), 127–150. 10.1016/0034-4257(79)90013-0

  • 75

    van BelzenJ.van de KoppelJ.KirwanM. L.van der WalD.HermanP. M. J.DakosV.et al (2017). Vegetation recovery in tidal marshes reveals critical slowing down under increased inundation. Nat. Commun.8 (1), 15811. 10.1038/ncomms15811

  • 76

    van MarenD. S.SchrijvershofR. A.BeemsterJ.ZhuC.XieD.ZhouZ.et al (2025). Land reclamation impacts on tidal landscape evolution. Rev. Geophys.63 (3), e2024RG000860. 10.1029/2024rg000860

  • 77

    VerbesseltJ.UmlaufN.HirotaM.HolmgrenM.Van NesE. H.HeroldM.et al (2016). Remotely sensed resilience of tropical forests. Nat. Clim. Change6 (11), 1028–1031. 10.1038/nclimate3108

  • 78

    VermoteE. (2021). MODIS/terra surface reflectance 8-day L3 global 500m SIN grid V061. NASA EOSDIS Land Process. Distrib. Act. Arch. Cent. (DAAC) Data Set MOD09A01. 061. 10.5067/MODIS/MOD09A1.061

  • 79

    WanS.QinP.LiuJ.ZhouH. (2009). The positive and negative effects of exotic Spartina alterniflora in China. Ecol. Eng.35 (4), 444–452. 10.1016/j.ecoleng.2008.05.020

  • 80

    WangC.TemmermanS. (2013). Does biogeomorphic feedback lead to abrupt shifts between alternative landscape states? an empirical study on intertidal flats and marshes. J. Geophys. Res. Earth Surf.118 (1), 229–240. 10.1029/2012jf002474

  • 81

    WangW.LiuH.LiY.SuJ. (2014). Development and management of land reclamation in China. Ocean and Coast. Manag.102, 415–425. 10.1016/j.ocecoaman.2014.03.009

  • 82

    WangX.XiaoX.XuX.ZouZ.ChenB.QinY.et al (2021). Rebound in China’s coastal wetlands following conservation and restoration. Nat. Sustain.4 (12), 1076–1083. 10.1038/s41893-021-00793-5

  • 83

    WangC.DingX.LianH.MaX.LiC.LiZ.et al (2025). Wetland habitat quality assessment and improvement from a perspective of spatiotemporal heterogeneity. Wetlands45 (7), 94. 10.1007/s13157-025-01976-y

  • 84

    WangY.LinJ.YanY.ZengZ.LiuX.WangS.et al (2025). Spatial transmission of resilience loss in deltaic dynamic systems. Commun. Earth and Environ.6 (1), 530. 10.1038/s43247-025-02519-3

  • 85

    WangB.LiA.ZhuM.LiuW.LuZ.LiuW.et al (2026). Ecosystem multistability maintained by water resource management in drylands. Commun. Earth and Environ.7 (1), 17. 10.1038/s43247-025-03043-0

  • 86

    WangR.LinZ.ChenC.WuW. (2026). The dual role of human activities: how reclamation and biological invasion reshaped blue carbon dynamics in an urbanized Estuary. Ocean and Coast. Manag.271, 107982. 10.1016/j.ocecoaman.2025.107982

  • 87

    WangT.ZhouX.SuL.ZhengS.LiuS.LiuY.et al (2026). Slow recovery of stable carbon constrains long‐term soil carbon restoration in China’s restored wetlands. J. Appl. Ecol.63 (6), e70449. 10.1111/1365-2664.70449

  • 88

    WeiX.CaiS.ZhanW. (2021). Impact of anthropogenic activities on morphological and deposition flux changes in the pearl river Estuary, China. Sci. Rep.11 (1), 16643. 10.1038/s41598-021-96183-0

  • 89

    WiensJ. A. (1989). Spatial scaling in ecology. Funct. Ecol.3 (4), 385–397. 10.2307/2389612

  • 90

    WoodwardF. I. (1987). Climate and Plant Distribution. Cambridge University Press.

  • 91

    XuH. (2006). Modification of normalised difference water index (NDWI) to enhance open water features in remotely sensed imagery. Int. J. Remote Sens.27 (14), 3025–3033. 10.1080/01431160600589179

  • 92

    YuanF.ChenW.-l.LiL.YangL.-z. (2025). Changes in Mudflat Resources in the Pearl River Estuary and Ecological Restoration strate8gies.

  • 93

    ZhangM.SchwarzC.LinW.NaingH.CaiH.ZhuZ. (2023). A new perspective on the impacts of Spartina alterniflora invasion on Chinese wetlands in the context of climate change: a case study of the Jiuduansha Shoals, Yangtze Estuary. Sci. Total Environ.868, 161477. 10.1016/j.scitotenv.2023.161477

  • 94

    ZhangQ.FanB.ZhangZ.LiY. (2025). Increasing fragmentation of coastal wetlands unveils large protection and restoration gaps. Ocean and Coast. Manag.270, 107890. 10.1016/j.ocecoaman.2025.107890

  • 95

    ZhangZ.ChiY.LiuZ.ZhangX.ZhangG. (2026). National coastal wetland mapping over the last four decades: an annual classification with high accuracy. Sci. Data13 (1), 250. 10.1038/s41597-026-06570-3

  • 96

    ZhaoH.CuiB.ZhangH.FanX.ZhangZ.LeiX. (2010). A landscape approach for wetland change detection (1979-2009) in the Pearl river Estuary. Procedia Environ. Sci.2, 1265–1278. 10.1016/j.proenv.2010.10.137

  • 97

    ZhaoW.LinS.YanX.ZhongJ.SuL.WuS.et al (2026). Spatially differentiated restoration strategies optimize multiple ecosystem functions in tidal marshes after spartina eradication. Ocean and Coast. Manag.273, 108073. 10.1016/j.ocecoaman.2025.108073

Summary

Keywords

estuary ecosystems, explainable machine learning, pressure-state-response, resilience, time series

Citation

Dong Q, Liu Y, Hu W, Arani BMS, Zheng Z, Xu C and Teng SN (2026) Remotely sensed vegetation resilience in the four major estuarine ecosystems of eastern China. Front. Environ. Sci. 14:1912598. doi: 10.3389/fenvs.2026.1912598

Received

18 June 2026

Revised

01 September 2026

Accepted

07 September 2026

Published

28 September 2026

Volume

14 - 2026

Edited by

Chong Jiang, Guangdong Academy of Science (CAS), China

Reviewed by

Biao Liu, North University of China, China

Hui Fu, Hainan University, China

Updates

Copyright

*Correspondence: Chi Xu, ; Shuqing N. Teng,

Disclaimer

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.

Outline

Figures

Cite article

Copy to clipboard


Export citation file


Share article

Article metrics