Biophysical factors and management practices are key to shaping forest resilience

Domain: chemistry
Published: 2026-02-16
Editorial outcome: major_revision

Abstract
Forest management influences key ecosystem services essential for societal well-being. However, how these human management activities affect forest resilience under different environmental conditions remains poorly understood. This study investigates how different forest management practices interact with biophysical factors and explores their joint effect on forest resilience globally during 2001–2015. To address this gap, we develop a framework that integrates Critical Slowing Down theory, satellite-based vegetation indices, and machine learning. Results show that natural forests without management exhibit the highest resilience among all forest management types, followed by natural forests with management and planted forests. Overall, forest management practices weaken resilience prominently under intensive anthropogenic pressures. However, specific biophysical conditions can reverse this pattern. Planted forests show higher resilience than natural forests when the ratio of precipitation to potential evapotranspiration exceeds 1.5 in wet climates, whereas natural forests remain more resilient in drier climates below this threshold. This threshold reflects the point where water availability becomes sufficient to compensate for the detrimental effects originating from human pressures. Cold temperatures, dense vegetation, and high soil fertility may further enhance the resilience of managed forests to levels typical of natural forests.

Introduction
More than half of the world’s forest area is subject to human influences. Humans can protect forests by establishing conservation districts and mitigating the impacts of natural disturbances. However, anthropogenic activities such as logging and even-aged monocultures can degrade ecosystem functions and lead to unexpected detrimental effects. Regional afforestation has emerged as a proactive strategy to combat forest degradation. As a key outcome of global afforestation efforts, about 22% of the forest coverage gains between 2000 and 2015 are anthropogenically planted. These large-scale afforestation programs play a vital role in restoring essential ecosystem services, including carbon sequestration, biodiversity conservation, and soil fertility maintenance. Planted forests (PF) are expected to further expand at an annual rate of 1.8% until 2050 and to be increasingly exposed to climate-driven disturbances owing to global warming. Understanding whether these rapidly expanding and potentially endangered PF can sustain resilience comparable to natural forests (NF) has therefore become a critical issue.

Managed forests (e.g., PF) can differ substantially from NF in structural and compositional characteristics, such as species richness and root depth, which may affect their sensitivity to climate variability and anthropogenic pressures, ultimately influencing their resilience. Forests with low resilience are typically more vulnerable to external disturbances, including natural (e.g., droughts and wildfires) and human pressures (e.g., logging), and are more likely to experience regime shifts. This issue is particularly critical in light of the ongoing intensification of forest disturbances, which increasingly threatens the provision of ecosystem services. Existing studies have prominently focused on the impact of climate change on forest resilience, whereas the potential role of human activities, such as forest management, has been largely ignored. Addressing this research gap is essential for developing sustainable forest management policies and achieving ambitious climate targets, such as the European Green Deal, the Decade on Ecosystem Restoration (2021–2030), and the Kunming-Montreal Global Biodiversity Framework.

Recent studies have used process-based models to predict forest resilience with moderate accuracy in certain regions, but their global predictions still exhibit considerable uncertainty, particularly under complex environmental conditions and in areas exposed to forest management, such as timber harvesting, reforestation, afforestation, and thinning. Despite substantial advances in such modeling tools, key biophysical processes are still only partially represented. Observation-driven techniques may offer opportunities to explore the impact of management practices on forest resilience. Such approaches appear particularly promising in view of the expanding availability of high-resolution satellite observations, the growing sophistication of statistical learning algorithms and the recent advances in cloud-based data storage and processing platforms.

Beyond the datasets analyzed, a clear framework for defining and quantifying resilience is essential. Earlier studies introduced the concept of vegetation recovery in response to disturbances (i.e., engineering resilience). This concept is tailored to specific disturbances (e.g., forest recovery to extreme droughts is characterized by specific threshold), as examined in controlled ecosystem experiments. However, forest ecosystems are exposed to a spectrum of external perturbations, ranging from discrete, pulse-like disturbances to progressive climate stressors. To address these complexities, we adopt the concept of ecological resilience, defined as the capacity of a system to resist (or absorb) changes from disturbances and to maintain its present structure or functions. In other words, ecological resilience reflects the ability to tolerate perturbations without shifting to a qualitatively different system. According to the Critical Slowing Down (CSD) theory, as systems approach tipping points where system states change, their sensitivity to disturbances increases and the corresponding capacity to absorb perturbations declines, indicating a loss of resilience. These changes in resilience can be reflected by early-warning signals (EWS), such as increases in lag-1 autoregression (AR1) and variance (Var), as well as their revised formulations (\lambdaAR1 and \lambdaVar) that account for the underlying non-linear relationship with ecological resilience. Despite the growing popularity of such techniques, there is still a substantial lack of global-scale assessments specifically targeted at exploring the combined role of human and biophysical factors in shaping forest dynamics and the underlying resilience processes.

In this study, we investigate how different forest management practices interact with biophysical factors (e.g., climate and vegetation characteristics) and explore their joint effect on forest resilience globally during 2001–2015. To this end, we integrate global forest management datasets, multi-source satellite observations, and climate products within a CSD-related modeling framework. Forest resilience is quantified in terms of AR1, \lambdaAR1, and \lambdaVar, and its spatial and temporal variations analyzed at 1 km spatial resolution. To disentangle the impact of human activities and biophysical factors, we develop an approach based on a paired-grid algorithm that compares geographically proximate forest areas under different management practices but similar environmental conditions. Furthermore, we employ machine-learning techniques to elucidate the non-linear interactions between biophysical factors and resilience difference across forest management types.

Results
Forest resilience across management types
We derive the distribution of forest management types for 2000 (GFMD2000) as the reference map. We generate this reference by overlaying tree cover data from the Global Hansen Forest Change Map for 2000 (2000 GHFC) with the Global Forest Management Datasets for 2015 (2015 GFMD). Multiple satellite vegetation indices are collected for 2001–2015 to explore vegetation dynamics, including leaf area index (LAI), kernel normalized difference vegetation index (kNDVI), fraction of photosynthetically active radiation (FPAR), gross primary productivity from the Moderate Resolution Imaging Spectroradiometer (GPP_MODIS), and GPP from the Penman-Monteith-Leuning model (GPP_PML). Seasonality and long-term trends are removed from the original vegetation time series to ensure stationarity and then CSD indicators are calculated from the residual signals to characterize forest resilience. In this study, forest resilience is measured by three CSD indicators globally: AR1 (lag-1 autoregression), \lambdaAR1 (lambda lag-1 autoregression), and \lambdaVar (lambda variance). These indicators are computed over the entire observational period. Therefore, under similar ecological contexts, high (low) CSD values indicate low (high) average resilience. Additional details are provided in the Methods Section.

Results indicate that human-managed forest areas generally exhibit lower resilience than adjacent unmanaged natural forests (NF WOM). This pattern is particularly pronounced in the Amazon and Central Africa, where the \lambdaAR1 estimated from LAI (\lambdaAR1-LAI) in natural forests with management (NFWM) is substantially higher than that in adjacent NF WOM. Lower resilience is also observed in Indonesian oil palms, PF in the southeastern United States and agroforestry in Southeast Asia. Resilience estimated by other vegetation indices (kNDVI, FPAR, GPP_MODIS, GPP_PML) and alternative CSD indicators (AR1, \lambdaVar) exhibits similar spatial patterns.

To quantitatively evaluate the relationship between forest resilience and human management, we employ a paired-grid framework designed based on spatial moving windows. Within each moving window, pixels from two different management types are matched on specific thresholds for forest age, forest cover, and elevation, and their average resilience is then compared. This approach allows us to control for potential confounding covariates (e.g., climate and vegetation conditions) (see Methods). Furthermore, to better elucidate the role of human practices in shaping differences in resilience across forest management types, we classify the original forest management map into classes representing different levels of anthropogenic intervention, such as PF and NF. Thus, the difference in CSD indicators among forest classes (e.g., \Delta\lambdaAR1-LAI = \lambdaAR1-LAI - \lambdaAR1-LAI) is used to infer the underlying differences in resilience characteristics. For instance, \Delta\lambdaAR1-LAI > 0.01 indicates that the \lambdaAR1 derived from LAI in PF is 0.01 higher than that in NF on average, and thus PF exhibit lower resilience.

The largest resilience difference is estimated between NF WOM and Plantations (\Delta\lambdaAR1-LAI > 0.11), followed by NF WOM versus forests with management (FWM) (\Delta\lambdaAR1-LAI > 0.07), and NF WOM versus NF WM (\Delta\lambdaAR1-LAI > 0.06). The resilience difference between PF and NF (\Delta\lambdaAR1-LAI > 0.01) is lower than that observed between Plantations and NF (\Delta\lambdaAR1-LAI > 0.03). This suggests that Plantations, including all human-planted trees, such as agroforestry and oil palms, exhibit the lower resilience (than NF or PF) plausibly because they represent ecosystems particularly fragile to external disturbances. This is consistent with the general patterns that forest management types under more intensive human pressure are associated with lower resilience. Consistently, resilience in NF WOM is approximately 10% higher than in NF WM and about 20% higher than in Plantations.

Although high-impact management generally weakens forest resilience, targeted management practices with low anthropogenic pressures can have beneficial effects for enhancing forest resilience. For example, more than 40% of grids show higher resilience in PF than NF (i.e., \Delta\lambdaAR1-LAI < 0), and these grids are located prominently between 40°N and 60°N. Such divergent patterns in resilience difference between PF and NF also emerge at regional levels, for example, \Delta\lambdaAR1-LAI is below 0 (higher resilience in PF than in NF) in the northwestern United States, and above 0 (higher resilience in NF than in PF) in the southeastern United States.

The analyses described above have been replicated using 15 different combinations of CSD indicators and vegetation indices (5 vegetation indices (LAI, kNDVI, FPAR, GPP_MODIS, GPP_PML) \times 3 CSD indicators (AR1, \lambdaAR1, \lambdaVar)). The consistency across these combinations highlights the robustness of the observed resilience difference among forest management types.

Biophysical factors influencing forest resilience
Although NF generally have higher resilience than PF, this relationship varies with the biophysical context. This clearly emerges by exploring the interplay between resilience difference and dryness index, here represented by the ratio between precipitation and potential evapotranspiration (P/PET). We find that the resilience difference between PF and NF shows a consistent negative relationship with P/PET. This suggests that a wetter (drier) climate could minimize (amplify) the resilience gap between PF and NF or even reverse their relationship. When P/PET is less than 1.5, the resilience of PF (e.g., \lambdaAR1-LAI = -1.45) is weaker than that of NF (e.g., \lambdaAR1-LAI = -1.60). Conversely, when P/PET exceeds 1.5, the resilience of PF (e.g., \lambdaAR1-LAI = -2.31) surpasses that of NF (e.g., \lambdaAR1-LAI = -2.15). Results from various vegetation indices and CSD indicators suggest that the resilience transition occurs when P/PET is between 1.3 and 1.5. This threshold likely reflects the point where water availability becomes sufficient to compensate for the detrimental effects of human pressures, allowing PF to exhibit greater resilience than NF. However, this mitigating effect of wetter climates on resilience difference is not consistent across all forest-type pairs. It is important to note that the definition of PF varies across studies. For instance, some works have adopted a broad definition that includes all human-planted trees (i.e., Plantations), such as agroforestry, oil palms, PF with short rotation (rotation periods < 15 years), and PF with long rotation (rotation periods > 15 years). We test our methods using multiple aggregation levels for PF and observe similar patterns in the relationship between resilience difference and P/PET.

In addition, the resilience difference between PF and NF is also correlated with individual biophysical factors such as PET, P, and air temperature (Temp). P/PET exhibits negative correlations with resilience difference (r < -0.4), whereas PET and Temp are positively correlated (r > 0.2). Resilience difference shows weak correlation with P (-0.1 < r < 0.1). These findings indicate that in regions characterized by humid climates and low temperatures, the resilience gap between PF and NF is substantially reduced.

We further investigate the role of biophysical drivers in influencing the resilience difference between PF and NF (details for other managed types in Discussion) using an Extreme Gradient Boosting model (XGBT). Within this modeling framework, the targeted resilience difference is expressed as a function of a comprehensive pool of environmental predictors, including temperature (Temp), dryness index (P/PET), root depth (RD), above-ground biomass (AGB), below-ground biomass (BGB), cation exchange capacity (CEC), and soil texture (Clay and Sand). For each environmental variable, both the local background conditions (Location_Factors) and the difference computed between PF and NF (Delta_Factors) are considered as predictors. Different vegetation indices and CSD indicators have minimal impact on estimating resilience difference. Therefore, we focus here only on results obtained from differences in \lambdaAR1-LAI between PF and NF (\Delta\lambdaAR1-LAI). The XGBT model achieves an R² of 0.65 and a Kling-Gupta Efficiency (KGE) of 0.58 on the testing samples. Lower performance is found in predicting extreme high and low values of \Delta\lambdaAR1-LAI. We also build two additional models, each driven exclusively by Location_Factors or Delta_Factors. The model based on Location_Factors shows higher explanatory power (R²) compared to the model based on Delta_Factors. This suggests that local biophysical contexts are more important than the inter-class differences in predicting variations in resilience across forest management types.

The rankings of the Shapley Additive Explanations (SHAP) values suggest that climate variables, such as temperature and P/PET, are more important than vegetation factors, such as biomass and root depth. SHAP provides a diagnostic tool to elucidate the interplay between the target variable and a given predictor. The SHAP value reflects the effect of a given variable on the resilience difference (e.g., positive or negative), while the feature value indicates whether the variable is at a low or high level (e.g., low vs. high temperature) and shows how this level amplifies or dampens the resilience difference. For instance, the red points for Location_BGB (below ground biomass) correspond to negative SHAP values, implying that high below ground biomass strongly reduces the resilience difference between PF and NF. These results indicate that wetter and colder climatic conditions, higher above ground and below ground biomass, and higher soil cation exchange capacity (i.e., reflecting higher soil fertility) may minimize the resilience difference between PF and NF.

The CSD theory suggests that systems with lower resilience are more sensitive to external disturbances and therefore more prone to critical transitions. We therefore investigate the sensitivity of forests to water-related and energy-related anomalies, and link the emerging relationships to forest resilience patterns. Water-related factors include soil moisture (SM), Palmer Drought Severity Index (PDSI) and precipitation (P), whereas energy-related factors include potential evapotranspiration (PET) and air temperature (Temp).

Results show that forest sensitivity to water-related factors decreases with increasing P/PET, while sensitivity to energy-related factors increases. These patterns are consistent with the expected major controls for ecosystem dynamics: forests in arid regions are primarily constrained by water availability and more sensitive to water-related anomalies, whereas those in humid regions are energy-limited and influenced by energy-related anomalies.

Sensitivity differences between PF and NF are analyzed by applying the analogous paired-grid framework utilized for resilience difference. We find that PF exhibit higher sensitivity to water-related factors than NF (positive values). This suggests that PF are more vulnerable than NF when exposed to water-related anomalies (e.g., droughts), and such fragility can be associated with a general reduced capacity to absorb disturbances (low resilience). Such vulnerabilities are most pronounced in water-limited regions, where water availability is the dominant constraint on vegetation fluctuations. By contrast, although water-related sensitivity differences may be even larger in humid regions, water is not a limiting factor there, so a higher sensitivity is not necessarily associated with lower resilience.

In energy-limited conditions, the abundant water supply allows forests to absorb quickly from water-related anomalies, making vegetation–water sensitivity a less critical determinant in maintaining stability and resilience. Sensitivity to energy-related factors becomes more important for maintaining ecosystem stability and resilience. Since PF exhibit lower sensitivity to energy-related factors than NF (negative values), they are less affected by energy-related fluctuations, supporting their ability to maintain stability (high resilience).

The observed transition from water to energy limitation with increasing P/PET helps to explain why the resilience difference between PF and NF varies across environmental contexts. PF are more sensitive to water-related anomalies but less sensitive to energy-related anomalies than NF, making them less resilient in water-limited regions but more resilient in energy-limited regions. Overall, these findings clearly show the key role of the regional biophysical context, reflected by P/PET, in modulating the resilience patterns of PF and NF.

Resilience gains and losses in planted and natural forests
Whereas the climate sensitivity analysis clarifies the potential mechanisms driving PF-NF resilience difference, it remains unclear whether these differences are stable or evolving over time. To explore this issue, we estimate the resilience trends from 2003 to 2013 using a four-year moving window. Results reveal a general decline in the resilience of global forests, consistent with previous studies. We investigate the temporal changes in resilience difference between PF and NF. Specifically, in drier regions (P/PET < 1.5), PF initially exhibit lower resilience than NF (\Delta\lambdaAR1-LAI > 0), but this positive difference decreases at a rate of -0.003 yr⁻¹. In wetter regions (P/PET > 1.5), PF have higher resilience than NF (\Delta\lambdaAR1-LAI < 0), but this negative difference shows an increasing trend of +0.002 yr⁻¹. Although the trends differ in direction across dry and wet regions, both indicate that the resilience gap between PF and NF became smaller from 2003 to 2013. Results by different vegetation indices and resilience indicators exhibit similar patterns. In water-limited regions (P/PET < 1.5), both PF and NF show a decline in \lambdaAR1-LAI from 2003 to 2008, followed by an increase from 2008 to 2013. The recent upward trend is stronger in NF (with a \lambdaAR1-LAI trend of +0.013 yr⁻¹) than in PF (+0.005 yr⁻¹), indicating that NF have become increasingly vulnerable to resilience loss in drylands. These results highlight the urgent need for targeted protection of natural forests in such regions.

In contrast, in wetter regions (P/PET > 1.5), both PF and NF experienced a decrease in \lambdaAR1-LAI during 2003–2013, with larger negative trends in PF (-0.012 yr⁻¹) compared to NF (-0.003 yr⁻¹). This suggests that resilience is increasing more prominently in PF than in NF under humid conditions. This pattern further suggests that afforestation in wetter regions can help narrow the resilience gap with NF.

Methods
Overall methodological framework
We explore the joint effects of forest management practices and biophysical factors on forest resilience at the global scale over the period 2000–2015 within a data-driven framework. The approach is structured into sequential modeling steps. (1) We combine forest management and tree cover maps to identify stable forest areas over the observational period. (2) We collect multiple satellite-based vegetation indices and extract residuals by removing seasonality and long-term trends. The resulting signals are then used to calculate forest resilience based on CSD indicators. (3) We implement a paired-grid moving window method to compare resilience between different forest management types under similar environmental conditions. (4) We employ machine learning tools to quantify the marginal effects of biophysical factors on resilience difference. (5) Monthly water- and energy-related climate variables are extracted, and vegetation–climate sensitivity coefficients are derived using bivariate linear regression to assess their relationship with forest resilience and identify possible threshold mechanisms. (6) A suite of sensitivity analyses is conducted to test the robustness of our results under varying definitions, modeling settings, and datasets. The workflow and abbreviations used in this study are shown in Supplementary Fig. 34 and Supplementary Table 1, respectively.

Forest management map
We use the Global Forest Management Data (GFMD) to classify major forest management types. GFMD combines crowdsourced points from the Geo-Wiki platform and satellite images to classify global forests into six management types at 100 m spatial resolution for the reference year 2015. The six forest management types include agroforestry, oil palms, short-rotation planted forests, long-rotation planted forests, natural forests with human management, and natural forests without management (Supplementary Table 2). Agroforestry is a forest ecosystem that integrates trees with crops. Oil palms are tropical monoculture plantations for oil production. Short-rotation planted forests use fast-growing tree species and are typically harvested within 15 years for timber, pulp, or bioenergy production. Long-rotation planted forests have rotations exceeding 15 years and are primarily intended for ecological restoration, carbon sequestration, and biodiversity enhancement. Natural forests with management experience human interventions such as logging, thinning, and cutting, while natural forests without management are pristine ecosystems (e.g., Amazon forest) that preserve their original biodiversity and ecological processes. Further details on specific forest management types are provided in Supplementary Table 2. The GFMD map is resampled to the 1 km spatial resolution to match the spatial resolution of vegetation data utilized to retrieve forest resilience.

We use the global tree cover for the year 2000 derived from the Global Hansen Forest Change product (GHFCv1.10), originally provided at 30 m spatial resolution and here resampled to 1 km, as a reference vegetation map for the start of our observational period. A pixel is considered forest if tree cover exceeds 50%. Following previous studies, we overlay the GHFC map (forest cover in 2000) with the GFMD map (forest management types in 2015) to identify stable forest areas during the 2000–2015 period. We only retain pixels that are classified as forest in 2000 and managed forest (of any kind) in 2015, assuming that the same management practice is applied throughout the entire period. This process ensures consistency in the forest types analyzed over 2000–2015 and rules out areas subject to land-cover change. We reclassify forest types based on a scheme specifically designed to compare resilience across different levels of human intervention (Table 1). Reclassified forest types include PF (long-rotation planted forests), NF (natural forests), NF WOM (natural forests without management), NF WM (natural forests with management), F WM (forests with management) and Plantations (human-established forests for productive or specific purposes).

Climate factors and auxiliary land surface products
To estimate the vegetation sensitivity to climate factors and explore the relationship between climate and forest resilience (detailed in the following sections), we use monthly mean precipitation (P), maximum and minimum air temperature at 2m, soil moisture (SM), PDSI, and potential evapotranspiration (PET, based on the Penman-Monteith equation). These climate variables are derived from the TerraClimate database at 4 km spatial resolution and for the period 2001–2015. Furthermore, we derive the mean temperature at 2m (Temp) as the average of maximum and minimum air temperatures and the dryness index as the ratio of P to PET (P/PET). TerraClimate integrates multiple meteorological datasets, including the Climatic Research Unit (CRU), the Japanese Reanalysis (JRA-55), and WorldClim, through meteorological interpolation, and its accuracy and reliability have been widely documented in previous studies. TerraClimate is selected as the reference climate dataset for this study as it provides a higher spatial resolution compared to other commonly used datasets, such as ERA5-Land, ensuring a spatial detail comparable with satellite-derived vegetation observations.

Furthermore, we explore additional environmental variables (static datasets) as the potential drivers of resilience. The root depth (RD) is derived from the Hydrologic regulation of plant rooting depth dataset provided at 1 km resolution. Above- and below-ground biomass (AGB and BGB, respectively) are derived from the NASA ORNL Global Biomass dataset for the year 2010 at 300 m resolution. Soil properties, including clay fraction (Clay), sand fraction (Sand), and cation exchange capacity (CEC), are obtained from the SoilGrids 2.0 at 250 m resolution.

All climate and auxiliary land surface data are resampled to 1 km through bilinear interpolation.

Vegetation dynamics
Different vegetation indices are used in this study as state variables to characterize spatial and temporal changes in forest resilience (methodological details about the resilience framework are presented in the next section). We select five vegetation indices including: leaf area index (LAI), fraction of photosynthetically active radiation (FPAR), gross primary productivity from Penman-Monteith-Leuning model (GPP_PML), gross primary productivity from Moderate Resolution Imaging Spectroradiometer (GPP_MODIS), and kernel-based normalized difference vegetation index (kNDVI). These vegetation indices are exploited to monitor the dynamic changes in vegetation from 2001 to 2015. LAI and FPAR are obtained from MYD15A2H.061 (Aqua) and MOD15A2H.061 (Terra) MODIS observations at 500 m spatial resolution and 8-day temporal resolution. GPP_PML is derived from the PML_V2.0.1.7 model. GPP_MODIS is derived from MOD17A2H.061 (Terra) and MYD17A2H.061 (Aqua) MODIS observations. Both GPP products (PML and MODIS) are released at 500 m spatial resolution and 8-day temporal resolution. The metric kNDVI, recently developed to address the saturation of traditional spectral-based vegetation indices in dense vegetation regions, is derived from NDVI based on the formulation:

kNDVI = tanh(NDVI^2)

In this work, NDVI data are derived from MOD13A2.061 (Terra) and MYD13A2.061 (Aqua) MODIS observations at 1 km spatial resolution and 16-day temporal resolution. Additional details on these products can be found in Supplementary Table 3.

We exploit vegetation retrievals from both MODIS satellites (Terra and Aqua) to mitigate the uncertainties caused by clouds, snow, and shadowing. Furthermore, we only retain high-quality pixels, according to the quality flag of the specific products, and synthesize monthly data using the maximum value composite (MVC) method. Specifically, for each calendar month, the maximum value of the vegetation index (e.g., LAI) is selected from all available observations within that month (e.g., from both Terra and Aqua). Each composite value is then assigned to the first day of the corresponding month, resulting in a temporally regular monthly time series, regardless of the actual observation date of the maximum value. To ensure consistency with the forest management map, all datasets are resampled to the common 1 km spatial resolution.

Resilience indicators
To estimate forest resilience, we use multiple indicators of EWS related to the CSD theory. The scientific context and the theoretical approach documenting the interplay between these early-warning indicators and resilience are detailed in Supplementary Note 9. The estimation for EWS requires stationary time series of the state variable, in this study expressed in terms of LAI, FPAR, GPP_PML, GPP_MODIS, and kNDVI. For this purpose, the seasonal signal is removed from the original time series by subtracting the monthly mean values, and the remaining trends are subsequently removed through linear regression. In line with previous studies, before applying such detrending and deseasonalization, we do not perform any gap filling of missing values to avoid potential bias in resilience indicators. All operations are conducted separately for each vegetation index.

Within the CSD framework, lag-1 autoregression (AR1) and variance are typically used as ecological resilience indicators. However, variance tends to be highly susceptible to environmental noise, and therefore potentially less reliable in detecting resilience changes. More recently, \lambdaAR1 and \lambdaVar have been proposed as complementary resilience indicators, to account more explicitly for the observational spatiotemporal resolution and observational noise. \lambdaAR1 and \lambdaVar are statistical indicators of resilience that are more closely aligned with the system’s actual ecological conditions. The mathematical formulations of \lambdaAR1 and \lambdaVar are described in the Supplementary Note 10. Considering the specific characteristics of the above indicators, we choose to evaluate forest resilience in terms of AR1, \lambdaAR1, and \lambdaVar.

To account for potential uncertainty in individual satellite vegetation products and CSD indicators, we generate an ensemble of resilience assessments based on the combination of five vegetation indices (LAI, FPAR, GPP_PML, GPP_MODIS, kNDVI) and three CSD metrics (AR1, \lambdaAR1, \lambdaVar). Such values have been used to robustly monitor forest resilience at the global scale over the period 2001–2015.

Forest resilience is quantified as a long-term ecological characteristic and derived from the CSD indicators computed over the whole time series (2001–2015) at each pixel (1-km) (e.g., 2001–2015 \lambdaAR1-LAI). This estimate describes the long-term average resilience of forest ecosystems. We complement this analysis by exploring the temporal evolution of forest resilience, deriving CSD indicators over 4-year moving windows, as proposed in recent studies. This results in a time series of resilience values over the period 2003–2013 at each pixel. Resilience dynamics are then expressed as the linear trend of the time-varying indicators (e.g., \lambdaAR1-LAI trend).

Paired grids using moving windows
To explore resilience difference across forests subject to varying levels of human interference, we implement a modeling framework to identify grids of geographically proximate forests under similar environmental conditions but different management practices. To this end, we apply a moving window method, which enables us to disentangle possible confounding environmental factors and determine the effects attributable to different management practices. To implement this method, we follow three key steps, as outlined below (taking PF vs. NF as an example).

(1) A spatial window of 25×25 pixels (approximately 0.25°×0.25°) is initially identified as an area characterized by overall similar background conditions. Here, pixels refer to the 1-km resolution units within the window (or the grid, used as a synonym in the text). PF and NF pixels within the window are then selected. Within each window, the CSD-based resilience metrics are first computed at the pixel level for all PF and NF pixels, respectively. The average of these pixel-level resilience values is then calculated for each forest type and used to represent the window-level resilience of the given forest type. To enable full spatial coverage and maintain continuity across the study area, a moving window approach is applied. The windows are shifted by 13 pixels in both x (longitude) and y (latitude) directions (e.g., half-overlapped), resulting in a partial overlap between adjacent windows, consistent with previous studies.

(2) Forest age strongly influences vegetation structure, while elevation affect microclimate and associated ecological processes. To mitigate the potential confounding effects of varying forest age and elevation on vegetation resilience, we ensure that the differences in these variables within each grid are restricted to within 10 years and 500 m, respectively. The global forest age data at 1 km resolution is obtained from Besnard et al., while elevation data at 30 m resolution is acquired from the Shuttle Radar Topography Mission (SRTM) v4 dataset.

(3) To ensure sufficient sample size, both compared management types (PF and NF) are required to each occupy more than 20% of the 25×25 pixels window (i.e., ≥ 125 pixels per type).

After obtaining the resilience values of different forest management types within each window, we compare them by calculating their differences (e.g., \Delta\lambdaAR1-LAI = \lambdaAR1-LAI - \lambdaAR1-LAI). We design six experiments to examine the resilience difference between pairs of different forest management types: NF vs. PF, NF vs. Plantations, NF WOM vs. F WM, NF WOM vs. NF WM, and NF WOM vs. Plantations (definitions of forest management types reported in Table 1). Notably, the number of grids may affect the robustness of our analysis. For each type of comparison, a minimum of 500 grids over the globe are analyzed to ensure statistical robustness and comprehensive coverage. Some combinations of forest management types have not been explored due to insufficient sample size (e.g., PF vs. NF WOM).

To facilitate the interpretation of results and more clearly capture possible emerging patterns associated with background climate conditions, resilience difference is explored across a P/PET gradient and correlated with additional climate drivers, such as temperature and precipitation.

The approach described above is replicated for each satellite vegetation product and CSD-based resilience indicator. The analysis exploits metrics of spatial and temporal variations in resilience, expressed in terms of long-term average CSD indicators (e.g., 2001–2015 \lambdaAR1-LAI) and trend (e.g., \lambdaAR1-LAI trend). Specific experiments to assess the sensitivity of our results with respect to thresholds, window width settings and additional environmental controls (e.g., disturbance regimes and ecological context) are described in a dedicated section later on (Sensitivity tests).

Climatesensitivity of forest ecosystems
Ecosystems that show high climate sensitivity, thus strong responses to water- or energy-related factors, generally tend to be more vulnerable to associated disturbances and possibly less able to maintain stable functions or structures (i.e., lower resilience). To investigate this issue, we examine the climate sensitivity of forest ecosystems and link these results to the emerging patterns in resilience difference between PF and NF across P/PET gradients.

We extract monthly water- and energy-related variables from the TerraClimate dataset from 2001 to 2015, including soil moisture (SM), mean temperature at 2m (Temp), precipitation (P), PDSI, and PET. Following the same data processing applied to vegetation indices, we remove seasonality and long-term trends from the climate variables. Because vegetation and climate variables differ in units and magnitude, we further standardize them to enable direct comparison of sensitivities. Each variable is standardized by dividing its monthly mean anomaly by the corresponding monthly standard deviation. Consistent with previous studies, we apply a series of bivariate linear regressions, each one using a standardized vegetation index (e.g., LAI) as dependent variable (Y) and one standardized climate factor (e.g., SM) as independent variable (X). The regression equation takes the form:

Y = b + aX

where a represents the sensitivity of the forest ecosystem to the selected climate factor, and b is the intercept of the regression. Following the same methodological approach used to derive resilience difference, sensitivity differences between PF and NF are estimated using the paired-grid approach over sliding windows. This method allows us to disentangle possible confounding factors and attribute the differences in sensitivity to variations in forest management practices. Sensitivity differences are then explored across a P/PET gradient and compared to the analogous observed patterns in resilience difference to elucidate possible underlying mechanisms.

Machine learning
We explore the biophysical drivers of resilience difference between PF and NF within a machine learning framework by employing Extreme Gradient Boosting (XGBT) and SHAP values. XGBT is a powerful ensemble learning algorithm widely utilized in ecosystem studies that combines gradient-boosting and tree-boosting steps. The XGBT framework developed in this study uses input data spatially aggregated at the window scale (25×25 pixels, approximately 0.25°×0.25°). The difference in long-term average resilience (e.g., \Delta\lambdaAR1-LAI) between PF and NF represents the target variable. Because the resilience differences derived from various vegetation indices and CSD indicators are strongly correlated (Supplementary Fig. 4), we focus our analysis on \lambdaAR1-LAI. As predictors, we consider eight biophysical factors, represented based on two computing schemes: (1) their average values within each window (Location_Factors), and (2) their between-type differences, calculated as the mean value in PF pixels minus that in NF pixels within the same window (Delta_Factors). The eight environmental factors include temperature (Temp), dryness index (P/PET), root depth (RD), above-ground biomass (AGB), below-ground biomass (BGB), clay fraction (Clay), sand fraction (Sand), and cation exchange capacity (CEC). These factors, including key climatic, vegetation and soil characteristics, are selected for their documented influence on forest resilience.

Based on the paired-grid method, we select 8220 grids reporting resilience difference values between PF and NF over the globe. The dataset is randomly split into 70% for model training and 30% for testing. To prevent XGBT overfitting, we employ 5-fold cross-validation to identify the optimal hyperparameters. Model performance is evaluated in terms of R-square (R2) and Kling-Gupta efficiency (KGE) index. To compare the relative importance of background biophysical factors and difference factors in predicting \Delta\lambdaAR1-LAI, we develop two additional XGBT models: the first one is driven only by Location_Factors, whereas the second one is driven only by Delta_Factors (Supplementary Fig. 9). We develop a set of additional XGBT models to explore resilience difference among other combinations of forest management types, including NF vs. Plantations, NF WOM vs. F WM, NF WOM vs. NF WM, and NF WOM vs. Plantations (Supplementary Fig. 17).

SHAP values are utilized to elucidate the data-driven relationships emerging from the XGBT model and quantify the marginal contribution of each biophysical factor. Positive (negative) SHAP values indicate a positive (negative) effect of the predictor on the target. The direction and magnitude of the effect depend not only on the SHAP value but also on the value of the feature. For example, when a feature has a high value and its corresponding SHAP value is positive, this means that higher levels of this feature are associated with a larger resilience difference. Conversely, a high feature value with a negative SHAP value indicates that higher levels of this feature tend to reduce resilience difference.

Sensitivity tests
To ensure the robustness and generalizability of our findings, we conduct a series of sensitivity tests, as briefly described in the following paragraphs.

Thresholds in the paired-grid method
We test the impact of different moving window widths (0.25°, 0.35°, and 0.45°), thresholds of forest age differences (5, 7, and 10 years), and forest coverage (10%, 15%, and 20%) on the resilience difference (e.g., \Delta\lambdaAR1-LAI). Results are shown in Supplementary Figs. 18–20.

Definition of planted forests
The narrow definition of planted forests used in our assessment includes only long-rotation PF. We recalculate the resilience difference by adopting a broad definition of planted forests (i.e., human-established forests for productive or specific purposes), including agroforestry, oil palms, long-rotation and short-rotation planted forests. Results are shown in Supplementary Fig. 7.

Effects of disturbance regimes and ecological contexts on the paired-grid method
We test the paired-grid method, and thus the assessment of the resilience difference between forest management types, with respect to varying disturbance regimes and ecological contexts. For this purpose, we exploit the Global Forest Cover Change (GFCC) product to retrieve the disturbance intensity at the pixel level. GFCC provides forest cover at 30 m spatial resolution over the period 2000–2015 at 5-year interval. The overall change in forest cover during 2000–2015 (e.g., FCC = FC2015 - FC2000) is used to derive the reference disturbance regime map. Furthermore, we use AGB retrievals derived from NASA/ORNL at 300 m spatial resolution, as a representative characteristic of the ecological background. Both disturbance regime and ecological context maps are resampled to the common 1-km spatial resolution. Similar to the constraints applied to elevation and forest age in the paired-grid design, we further constrain FCC differences between the two paired forest management types to within 3–10% and AGB differences within 5–10 Mg/ha. For example, we select the grids where the difference in FCC between PF and NF is within 10%, 5%, or 3%, or where the difference in AGB between PF and NF is within 10 Mg/ha or 5 Mg/ha. Results are shown in Supplementary Fig. 21.

Missing values and data pre-processing
We test the sensitivity of resilience estimation to different strategies for temporal aggregation, sensor selection, and missing data. Based on the pre-processing in previous studies, we test the impact of three different pre-processing methods to characterize vegetation dynamics: (1) retrieval of the monthly maximum vegetation data via MVC from both MODIS Terra and Aqua satellites; (2) original temporal resolution of vegetation data from Terra satellite; and (3) original temporal resolution of vegetation data from Terra satellite restricted to the growing season, defined as months with monthly temperature > 5°C. Applying temporal lags to growing-season-only data creates one missing value per year as the first growing-season month shifts beyond the data range (i.e., missing data). The analysis focuses on the resilience difference among different forest management types as well as the spatial patterns of resilience, and results are compared against the original setup based on monthly MVC data. Results are shown in Supplementary Figs. 22–25.

Standardization and normalization methods
We test three procedures for \lambdaAR1 retrievals based on: (1) kNDVI anomaly, obtained by subtracting the long-term monthly mean from each monthly value; (2) kNDVI Z-scores, obtained by subtracting the long-term monthly mean and dividing by the corresponding monthly standard deviation; and (3) normalized kNDVI anomaly, obtained by subtracting the long-term monthly mean and then rescaling each monthly anomaly to the [0, 1] range based on the minimum and maximum anomalies of the same month. Results are shown in Supplementary Fig. 26.

Spatial resolution of datasets
Analyses for resilience estimation and differences among forest management types originally conducted at 1-km spatial resolution are replicated at 500 m for all different vegetation products. Specifically, we repeat the resilience estimation and moving window analysis at 500 m for the full set of vegetation indices (LAI, FPAR, GPP_PML, GPP_MODIS, and kNDVI), GFMD, and GHFC datasets. For kNDVI (NDVI), which is originally available at 1 km resolution, we additionally employ MOD13A1 and MYD13A1 products at 500 m resolution aggregated to monthly values for the sensitivity test. Results are shown in Supplementary Fig. 29.

Timespan and temporal frequency
The effects of time span and temporal frequency on resilience estimation are evaluated by comparing AR1-kNDVI values retrieved from 16-day and monthly datasets, as well as from different study periods (2001–2010, 2001–2015, 2001–2020, and 2001–2024). Results are shown in Supplementary Fig. 30.

Uncertainty in the forest management map
The GFMD is potentially subject to uncertainty originating from the underlying classification algorithm. To test the impact of such a source of uncertainty, we conduct a series of experiments by excluding pixels characterized by low classification accuracy (tested accuracy thresholds: 50%, 60%, and 70%). For example, we remove vegetation management types with classification accuracies below 70% and evaluate the corresponding resilience difference. Results are shown in Supplementary Figs. 27 and 28. It is important to note that the GFMD map is derived from a classification algorithm that uses various factors, including LAI, which is also employed in this study as a vegetation index to assess forest resilience. This circularity may undermine the independence between types of forest management and resilience. To better evaluate this potentially critical issue, we collect ~210,000 crowdsourced validation samples (N_NFWOM = 39,534, N_NFWM = 82,358, N_PFwithoutrotation = 8357, N_PFwithrotation = 20,247, N_Agroforestry = 47,751, and N_oilpalms = 11,691) from crowdsourced datasets in 2015. These ground-truth reference samples are obtained via expert-verified and crowdsourced campaigns on the Geo-Wiki platform. We exclude crowdsourced datasets characterized by non-forest pixels in the year 2000 based on GHFC. The remaining crowdsourced datasets are then used to analyze the variation in resilience across P/PET ranges among different forest management types. Results are compared with analogous results obtained from the GFMD and are shown in Supplementary Figs. 31–33.

Moving window method
We explore an alternative method based on the nearest-neighbor matching algorithm to select two paired forest management types and explore resilience difference, using globally crowdsourced datasets. For instance, each PF point is matched with its n nearest NF points (n=5, 10, 15) based on environmental similarity. This approach provides an alternative perspective to the moving window method by focusing on pairwise comparisons. Results are shown in Supplementary Fig. 33.

Across all tests, the resilience difference between PF and NF and its relation with P/PET appear consistent with results obtained from the benchmark modeling setup described in previous sections. Minor deviations emerge from the nearest-neighbor pairing with crowdsourced samples, plausibly due to sample imbalance/label heterogeneity (see Discussion for details). Overall, the sensitivity experiments conducted confirm the robustness of the proposed methodology and the validity of results. Additional details on sensitivity tests are described in Supplementary Notes 3-6.

Peer review rounds
[{"round": 1, "reviewer_comments": {"reviewer_1": "## Overall impression\nThe manuscript attempts an analysis of the sensitivity of forest resilience to human disturbance and environmental drivers in\ndifferent forest management categories, arguing that natural forests have higher resilience than managed forests in 60% of\nthe cases. While the subject is of great interest, the manuscript is unfortunately not convincing. The techniques used are\ninappropriately explained and the limitations of the underlying data are not well exposed. Furthermore, the overall quality\nand style of the manuscript is not in-line with what is expected in such a journal. References are incomplete, lacking journal\nnames, sometimes pages. There are quite some errors in language that should have been picked up by any spelling\ncorrector. Some phrases seem unfinished (no verb, or no object, see L93). Then technically too, various crucial\nmethodological details are not explained in the main text, making understanding and following the text quite difficult, and\nundermining the credibility of the claims that the manuscript tries to put forward. Given these various points, the\noverall recommendation would be rejection. However, if the authors want to improve the manuscript,\nhopefully this review can help them. Below is a list of more general issues followed a more detailed list of specific points.\n## General issues\n- The abstract gives no clue on the nature of the work. How is resilience quantified? Model simulations? satellite metrics?\nground-data?\n- Consider that the methodology section comes at the end, so any concept, method, dataset, or abbreviation needs to be\nintroduced before being dropped on the reader without context. A few examples are listed in the detailed comments.\n- In the introduction, four warning signals are presented: AR1, variance, λAR1, and λVar. Yet, the results are based only on\nthree, excluding the variance metric. This needs to be motivated or risk feeling like cherry-picking.\n- Check that the references to the figures are correct throughout the document. It felt that some were not, which made\nfollowing along sometimes very confusing. Also, Figures and supplementary Figures should be always presented in order,\nand all should be referenced in the text.\n- Statements based on studies from the literature should probably be moved all to the discussion, as in the results it is not\nalways clear if the statements made are based on own results or the literature. Alternatively, it should be made very clear\nwhenever a reference is cited as basis for a statement, or as support of own conclusions.\n- The definition of forest types seems pretty central to the research question, so it should also be central in the text, but it is\nnow only relegated to the supplementary materials. Either tables S1 and/or S2 should be moved to the main manuscript, or\nthe most important forest type definitions should be clearly stated in the text.\n- The arguments about the independence between forest management and resilience of the GFMD data should be more\ncentral in the discussion, as it is potentially a huge underlying limitation. Therefore, the independent tests on crowdsourcing\ndata should be highlighted more.\n## Specific points\nL17: P/PEF should be defined\nL49: correct “bring” to “bringing”.\nL60: “However, despite models appearing to predict stably, forest dynamics in observation remain complex and diverse”\nunclear sentence.\nL62: correct “plating” to “planting”.\nL64: remove double comma after “Furthermore”; clarify what is meant with “early stable predictions”.\nL67: “Recent study report that…” correct singular/plural form (either “A recent study reports that…” or “Recent studies report\nthat…”.\nL68: correct “dectect” to “detect”\nL69-70: “…it will be better to…” this formulation sounds very strong. Both observational and modelling frameworks have their\nuses. I would present them as alternatives to one another, rather than stating that one is superior.\nL84-86: it might be worth spending a few more words/sentences to introduce and define the concepts of AR1, variance,\nλAR1, and λVar, as they seem rather central to the presented results.\nL90-93: “Motivated by the issue of what impacts of human activities especially forest management including PF and NF on\nvegetation resilience, this study used the latest global forest management dataset, multi-source satellite observations, and\nclimate datasets to investigate. Here, this study mainly to:” unclear sentences.\nL94: it might be worth introducing the “paired method” before mentioning it in the research questions.\nL106: “In Spatial”, meaning?\nL107: “decline” as per meaning become smaller, fewer, weaker, etc.\n(https://www.oxfordlearnersdictionaries.com/definition/english/decline_1) is an intransitive verb, meaning that it cannot have\nan object. In your formulation: “…human forest management activities have declined the resilience” you want to use a\ntransitive verb such as “decreased”. Please check the use of decline vs decrease throughout the manuscript.\nL108: Explain what the GHFC and GFMD are.\nL111: Rationale for these metrics is needed. We have no explanation of what they do, how they are constructed and what\nassumptions hide behind. This needs to be clear in the main text, not guessed or ignored. Also, in the introduction you talked\nabout four warning signals: AR1, variance, λAR1, and λVar. Why are the results based only on three, excluding the variance\nmetric?\nL112: is the starting time of the analyzed data 2000 or 2001? This seems inconsistent.\nL112-113: “The results revealed…” what results? It is not clear to me what analyses were performed. Are the statements in\nthis paragraph purely based on visual inspection, or is there more to it that is not mentioned? This should be clarified.\nL117-118: “Forest fragments experienced higher levels of human disturbances and pressures, such as logging, clearing, or\nagriculture 38” It should be clarified that these are results from the literature, if so.\nL123-124: “These… influence…” check consistency of singular/plural forms; “management…” practices?\nL125-127: Again, not entirely clear if these statements are based solely on studies from the literature or if they simply agree\nwith the literature.\nFig 1: confusing and not vizually clear. strangely structured. Axis titles are missing. Lat/lon indicators are missing. So far, LAI\nhas never been mentioned in the text, so the subplots b and d come as a bit of a surprise. Background grid lines could aid\nthe comparison of subplots b, d, e., and the x-axes could be aligned for more consistency and ease of interpretation. In the\nfigure description: “as an example of a specific global grid point”, you could say what point you are presenting (e.g.\ncentral coordinates). Finally, the last part of the figure caption probably belongs to the main text: “The other examples in Central\nAfrica… …found in Figure S2, showing similar spatial patterns”.\nL149: it is here we discover there are paired grids. But how are they paired? some info needs to percolate in the main text\nL149: “paired grids” again, it should be briefly mentioned before what this method is about.\nL168-169: “over 40% of paired grids show higher resilience in PF than in NF” this is quite high! If this is the case, I wonder if\nthe statements about management decreasing resilience is founded or whether you should be more conservative.\nConsidering the relevant uncertainties, is the 60%-40% ratio really different than a 50%-50% ratio based on random noise?\nL174-175: “P/PET (the Ratio of Precipitation to Potential Evapotranspiration)” it is good that you define P/PET, but you\nshould maybe also mention what it is since it is the first time you mention it in the text, i.e. an aridity index.\nL179: “these 15 resilience calculations…” not entirely clear what “these” refers to.\nL183: “Fig. S7”, what happened to figure S6?\nFigure 2: the color palettes chosen to represent the resilience difference and the datasets are very similar (browns and\ngreens especially). This makes the interpretation of the figure difficult, so the colors should be differentiated more.\nAlternatively, subplots b, d, and f should come with their own color legend. In the figure description (here and in other\nfigures) you mention a significance test, but it is not specified what was tested. The significant difference from 0?\nL221-223: “since… which…” correct the structure of the sentence. Currently the main clause is missing.\nL245: Mention what “background factors” are.\nL249: Correct “focued” with “focused”.\nL256: “explanatory” power?\nL259: “other factors”, which ones? Maybe make a few examples.\nFigure 4: I am a bit confused about SHAP values vs feature values. What do they represent and why do we need both? Are\nthe feature values similar to effect coefficient? It should be made clearer even for readers that are not very familiar with the\nmethodology.\nL281-285: It is a bit strange that subplots e (PDSI) and i (temperature) are presented before a, b, c, etc.\nL285: “…increases as the enhancement of…”, you mean “…increases with the increase of…”?\nL287-288: “In this regard, results further suggests that PF are more sensitive to water resources than NF (Fig. 5e-5j)” this is a\nbit hard to follow, since the results are presented via sensitivity difference. You could add a sentence to explain why this is\nthe case.\nL291: what happened to Figure S10?\nL292: sensitivity to “soil water” or “soil moisture” might be more appropriate than “water” in general, as SM is what is shown\nin Fig. 5c).\nL296-298: “Because the energy sensitivity difference between PF and NF gradually increases along with the wetter climate,\nand the variations of vegetation greenness is energy-limited in humid areas” missing a main clause.\nL298-299: “As the enhancement of P/PET” does not seem grammatically correct. Change to “With increasing P/PET” or “As\nP/PET increases”.\nFigure 5: Results descriptions and references to supplementary materials should probably go in the main text (results) and\nnot in the figure description. Clarify what the significance test is testing. Also, “Vegetation greenness is limited by water in\ndry conditions and by energy in wet conditions, which may be related to shifts in resilience differences between PF and NF”:\nto me, this sounds related to actual resource limitations rather than resilience differences.\nL327: λPF – λNF is a newly introduced metric/abbreviation. Explain what it is. It is also not found in the referenced figure 5.\nL332: “…PF and NF are reducing in both…” should be “decreasing” instead of “reducing”.\nL359: “However, a notable ambiguity arises from the definition of managed forests also came to this study” unclear sentence.\nL388: “CEC” is not defined so far.\nL399-401: repeated sentence.\nL416-418: unclear sentence.\nL420-422: “More interestingly, vegetation greenness is water-limited in dry regions (PF is more sensitive to water) and\nenergy-limited in wet regions (NF is more sensitive to energy)”. To me, this sounds rather obvious. It is well known and\nintuitive that vegetation is limited by water where water is scarce, and energy becomes the limiting resource where water is\nabundant. Maybe you want to reformulate this to frame it more in the context of resilience and sensitivity.\nL438: I believe one needs to reach this point to learn for the first time that these indices are made from LAI\nL451: “decline” change to “decrease”.\nL458: “the primary focus is on the differential impacts of forest management between PF on resilience” unclear sentence.\nL480: “…and reduce uncertainties from land use change…” and following sentence: to me, this sounds more like land use\nchange was ignored since you are only assuming that pixels do not change, but in reality, they might be changing. Am I\nmissing something?\nL491: “ASCE” define/explain what it means.\nL503: change “photosynthetic” to “photosynthesis”.\nL506: “PML” and “MODIS” you may want to explain the abbreviations or mention the source (which for MODIS you do a few\nlines below).\nL511: this is a wrong use of kNDVI. see https://doi.org/10.1016/j.isprsjprs.2022.12.019\nL516: “which are mainly from [MODIS]” and where are the other products from?\nL518-519: “these two observed vegetation greenness or biomass are similar in a day” give a reference stating this.\nL519: “But two times may be affected by the cloud” unclear sentence.\nL524-526: You mention two contrasting spatial resolutions for the vegetation indices (500 m and 1 km). Which one is it?\nL527: why are you limited to 1km by computational constraints? it seem very likely it could all be done at 500m\nL530: you still have not explained how to calculate resilience. why not start with that before?\nL531: If not already mentioned in the introduction, the “paired grids” concept should really be explained in detail here, as I\nam still wondering what it means.\nL538: what is Simon et.al?\nL541-543: mention supplementary figures explicitly.\nL546: why monthly data? to do a robust time series analysis it is better to take all data available, not picking a single MVC\nvalue per month\nL550: you need to explain what these are\nL552: It is still unclear to me what the difference between AR1 and λAR1 is.\nFormulas: try to keep the font of parameters consistent between the formulas and the main text.\nL571: “same as the Xn” unclear.\nL591-593: correct typing mistakes: “however”, “maxium”, “northern”.\nL591: Explain what “MVC” stands for.\nL592: Explain what “month datasets” are, or reformulate.\nL593: “influencing” reformulate, unclear.\nL595: “some CSD studies did not deal with missing values”, does that mean that missing values are\nexcluded or kept as they are? Either way, it is not possible to “not deal with missing values”.\nL603-605: “Not only because of the strong relationship between resilience difference between PF and NF and background\nvariables, …” missing main clause. Also, this part specifically sounds as if you chose the method based on the results\nbecause “it looks good”, instead of basing the decision on hypotheses or research gaps.\nL607: “widely utilized in ecosystem studies”, then maybe give a few more references instead of only one.\nL612: what does “its” refer to?\nL615” what does “their” refer to?\nL621: “datasets availability section” should be corrected to “Data and code availability section”.\nL671: “low probability pixels” what does this refer to, low probability of what?\nL681-683: “GFMD is the result of a classification algorithm that takes in input various factors including vegetation data, the\nsame data to calculate resilience. It may undermine the independence between forest management and resilience”. This\nsounds like a potentially huge issue. If the input data are not independent, how can you say anything about the results? The\nfollowing parts were you control with independent data should really be mentioned more centrally.\nL702-703: “This is aligned with our main conclusion. The above analysis results can be seen in the discussion and\nsupplementary section” good, as mentioned above, I think this confirmation should be more central to the results. Also, I\ncould not find it in the discussion, where it should really be, rather than hidden at the end of the methods or in the\nsupplementary materials.\nL707: the data here means the data produced by the paper, not just the data used as input (although the latter is also\nneeded)", "reviewer_2": "I co-reviewed this manuscript with one of the reviewers who provided the listed reports. This is part of the Nature Communications initiative to facilitate training in peer review and to provide appropriate recognition for Early Career Researchers who co-review manuscripts.", "reviewer_3": "This study explores the resilience of forests under varying management intensities, ranging from unmanaged natural forests\nto managed, planted forests and plantations. The research examines how resilience, as defined by the critical slowing down (CSD)\ntheory, differs among forests with different levels of management. The study concludes that natural, unmanaged\nforests exhibit greater resilience, while planted forests are more vulnerable and thus show lower resilience. Additionally,\nthe study investigates the factors influencing these differences, reporting that wetter regions are less vulnerable and,\nconsequently, more suitable for planted forests.\nI find the study both interesting and novel, contributing important findings beyond previous similar studies. However, several\npoints need to be addressed to strengthen the results. The authors conducted multiple sensitivity analyses to demonstrate\nthat the findings are consistent, considering factors such as the sliding window width, forest coverage threshold, definitions\nof planted (PF) and natural forests (NF), reliability of forest management maps, and preprocessing methods of satellite\nvegetation data. However, one crucial aspect that should be considered is forest age: natural and managed forests are\ntypically older than planted forests and plantations, and resilience is known to vary among different age groups. Therefore, I\nrecommend including an additional sensitivity analysis to compare planted forests with younger unmanaged stands, to verify\nthat the results are not merely due to age differences between PF and NF. Furthermore, in the final chapter discussing\npotential limitations, the authors should expand on the limitations of the CSD theory.\nAside from these points, I have only a few minor comments detailed below:\nAbstract: The abstract lacks specific mention of the aspect of resilience studied. I suggest mentioning the CSD theory in the\nabstract.\nIntroduction: The introduction is well-prepared, effectively outlining the concept of resilience, identifying the research gap,\nand clearly stating the study's objectives.\nMethods: The authors utilized validated datasets, such as GFMD, the TerraClimate database, and MODIS. Although the\nmethods involve some subjective choices, the sensitivity analysis provides reassurance about the robustness of the results.\nThe classification of resilience into stable, increasing, or decreasing is based on a well-established theory and is therefore\njustified.\nDetailed Comments:\nLine 106: The phrase \"In Spatial\" is unclear—please clarify.\nLine 108: Introduce the abbreviations GHFC and GFMD here (or in the Introduction) as they are used for the first time.\nLine 176: Could you provide an explanation for this pattern already here? Why are the results different in these regions?\nLine 183: Please check for grammatical accuracy—consider changing to \"showing similar patterns\" without \"the.\"\nLine 242: Instead of \"cycle,\" did you mean \"circle\"?\nLine 320: Please provide references to the \"previous studies\" that also documented a global decline in resilience.\nLine 425: I suggest adding a sensitivity analysis comparing younger NF and younger PN here. Additionally, consider\nmentioning the limitations of the CSD theory.\nLine 538: The reference \"Simon et al.\" is missing from the reference list.\nLines 614-615: Why didn't you use age as a predictor variable? Age could significantly impact resilience and potentially\naffect your findings. Please elaborate on this. Since PF are generally younger than NF, I recommend an additional\nsensitivity analysis comparing PF and NF of similar ages.\nLines 617-621: Although you mention that \"soil and vegetation variables\" are detailed in the dataset availability section, I\nbelieve these should be briefly described in the Methods section, as they are critical inputs for your model.\nLine 650: I couldn't find the reference to \"Taylor et al.\"—please ensure it's included in the reference list.", "reviewer_4": "Review report\nThe study is focused on the forest management practices (human factors) and the (natural factors) driving forest ecosystem\nresilience. There are sustainable forest management practices that support forest ecosystem resilience while other\nunsustainable practices undermine forest ecosystem resilience. This is a very interesting topic, especially as the resilience\nof forest ecosystems is critical to the supply if ecosystem services benefiting people. The effect of natural factors that also\ninfluence forest ecosystem resilience, and that is a key objective of the manuscript is not capture in the manuscript title?\nThe resilience of forest ecosystems is not only explained by forest management practices (type and intensity) but also more so,\nby biophysical factors, including climate change, and as such should be captured in the title.\nBesides, I have some general major concerns.\nComparing planted forests and natural forests seems odd for me, as the difference in their resilience is quite obvious from\ntheir management practices. The results stated in the abstract are quite obvious: For example, “forest resilience suffers\ngreater losses under intensive management, as seen in human plantations like agroforestry, which are much less resilient\nthan unmanaged natural forests”. And “Such resilience difference between PF and NF is related to background conditions\nsuch as local climate, vegetation, and soil”. Topography also plays a role? These are obvious well documented results and\ndo not provide any novelty, even where a quantified difference is provided. Is the resilience of well managed planted forests\nat the same level as or even more than unmanaged natural forests?\nMy suggestion is to compare similar forest ecosystem types (e.g., managed and unmanaged natural forests; planted\ncommercial forests and planted forests mainly aimed to restore ecosystem functions and service, etc.). Comparing similar\nforest ecosystem types with different or similar forest management practices and biophysical conditions would be more\nintriguing in understanding differences in their resilience and inform restoration and conservation decision making.\nComparing within planted forests with primary objective for restoration functions and services and planted forests under the\ncommercial forestry approach can enhance our understanding of the resilience of these ecosystems in terms of supply of\necosystem services. The resilience of planted forests monocultures and multiple tree species forests can be compared to\nfurther understand why forest monocultures are discouraged.\nComparing managed and unmanaged natural forests can provided important insight on the effect of forest silvicultural\nmanagement practices on forest ecosystem resilience. Unmanaged old growth forest may sequester less carbon due to\nolder age trees compared to young secondary regrowth forests. Do unmanaged natural forests need some sought of\nmanagement intervention to increase resilience or should they be kept intact?\nThe resilience of forest ecosystems is not only explained by forest management practices (type and intensity) but also more\nso, by biophysical factors. But it is clear this was not investigated in this study, making the scope so one dimensional. This is\nclearly acknowledged by the authors. For example:\nIn line 446 to 448, the authors clearly stated that “…the narrowing disparity in resilience between PF and NF may be\nattributed to multiple factors, including vegetation physiology, climate change, and the intensity of human activities. This\ncomplex issue merits further investigation”.\nIn lines 459 to 460, “… but there is less exploration of whether the climate, soil, and vegetation background can mitigate the\neffects of other forest managements”."}, "author_response": ""}, {"round": 2, "reviewer_comments": {"reviewer_1": "## Overall impression\nThe authors have made a remarkable effort to make the manuscript more legible. They have effectively answered several\nconcerns of the reviewers. However, now that it is clearer (notably on the metrics used), my overall feeling is that there is\ninherent flaws in this study. Unfortunately, these flaws confirm to me the first impression that it should not be accepted. I will\nelaborate below my main concern, and the expose some of the more specific issues I encountered while reading the\ndocument.\n## General issues\nThe main problem I have with this study relates to the use of the CSD theory specifically on these global vegetation time\nseries of remote-sensing derived indicators in view of characterizing alleged differences in resilience between natural and\nplanted forests (i.e. the main goal of the study). As the authors state (L85) \"According to CSD theory, as systems approach\ntipping points where system states\nchange, their sensitivity to disturbances increases, and recovery from disturbances slows down, indicating a loss of\nresilience\". To begin with, it would be necessary to state specifically how these concepts are specifically mapped to the\nstudy case at hand of forests. First, what are the tipping points we are considering? Is it \"collapse\" of the forest? how is that\ndefined specifically for PF and NF? Defining the tipping point seems necessary to me so as to say what signal we should be\nusing to monitor the disturbances for establishing the characterization of resilience. Then, the next question is: what\nspecifically are we considering as \"disturbances\"? Windthrows? Fires? bark beetle infestations? clear-cuts? Later in the\ndiscussion it states (L500) \"abrupt events and human disturbances lead to non-stationary conditions\" and that this limits the\naccuracy. So I do not exactly understand what disturbances we are trying to see in the signal. This adds questions then to\nwhether the signal used are appropriate in terms of (1) what they measure, (2) over what spatial supports, and (3) with what\ntemporal frequency. For point (1), all metrics used fundamentally relate to temporal variations in greenness (whether they be\nLAI, kNDVI, fAPAR, GPP, they are reflect mainly greenness and the information comes mostly from NIR and Red\nreflectances). Are changes in greenness over 5 by 5km (point (2)) detectable with a frequency of one or twice a month (point\n(3)) really capturing the target \"disturbances\" that are not man-made and not abrupt? Is this the approapriate tool for this\nquestion?\nTo bring a bit more of a specific detail of why this is problematic for this question, consider that managed forest will be more\nsubjected to clear cuts by definition. These will entail stronger drops in greenness than for natural forests, so when this is\naveraged of coarse spatial resolutino of 5km, the mixed signal will be more strongly affected (on average) for the managed\nforest than natural forest. First, is any normalization used here? will the AR1 metric not be artificially affected by this?\nThen, the recovery could be faster in a planted forest if fast-growing species are planted (and assisted, with protection from\nherbivory for instance) than if natural regeneration occurs. This would be reflected in the metric by suggesting managed\nforest would have higher resilience (as faster recovery), and yet this is counter to what the results suggest. If such clear-cuts\nevents should not be considered (as suggested by some phrases in the text cited above), then how is this prevented from\naffecting the results and the conclusions?\nI can go on with more concerns, such as the questioning of whether there are enought samples in time to have sufficient\nfrequency to sample the stable state as suggested in Figure 2c and 2d of Scheffer et al (2015), mentioned by the authors as\nthe reference. With monthly data over 20 years, I have some serieous doubts that this can apply so seamlessly. I would think\nyou would then need daily data from satellite, which is of course available, but then you need to deal with changes in viewed\nangles and missing data, which I suppose incapacitate the use of these AR techniques. Furthermore, the MVC technique to\ncomposite means the actual time stamps of each month is not regularly separated from the former or from the next, so this\nfurther means the AR techique is applied in sup-optimal conditions (i.e. non regular timesteps). All these issues are not\naddressed.\nI know others have used these metrics on such similar data before (e.g. Forzieri et al.), but that does not convince me either,\nas I believe the same flaw applies there. I would require a proper rebuttal plus explanation to address these specific issues\nto accept the premises of this manuscript. Below I further provide more (non exhaustive) specific points to help the authors\nimprove the text.\n## Specific points\nL73: Management should not be plural. Perhaps change to \"types of management\"\nL82: CSD should be spelt out here\nL91: These new metrics should be explained, even if briefly, not just cited. The reader should get here the diffrence between\nAR1 and lambda_AR1. What does the lambda stand for? What is Lambda variance?\nL94: what is \"variance\" by itself? you need to be more specific... the temporal variance of the signal... and over how long of a\nperiod?\nL97: In this paragraph you speak of statistical metrics, but you should have also mentioned the signal on which the metrics\nneed to be applied, otherwise it is pointless. They cannot be applied to a signal that does not speak about resilience. So first\nmention how RS signals can effectiveely (or not) be sensitive to properties that speak to resilience, and only then introduce\nthe statistical metrics.\nL104: Not clear why ML comes into play? Why would it be needed? needs to be introduced. Also the \"impacts of biophysical\nfactors\". What are these \"biophysical factors\"?\nL110: \"biophysical factors\" should have been mentioned and described before\nL117: why tree cover from 2000 and management for 2015? the former is available for 2000\nL122: kNDVI is with lowercase k\nL127: again, the acronyms are meaningless and useless if we do not know what they measure\nL499: Here is a line that puzzles me. If we need stationary behaviour, and \"abrupt events and human disturbances lead to\nnon-stationary conditions\"\nL502: This scale dependence is crucial to take into account in my opinion as in this case, natural forest will likely have\ndifferent (smaller) spatial extents of perturbations than managed forest (where, by definition, there would be more clear cuts).\nL508: GPP is not a biomass index. And the way it is calculated (from MODIS) is very likely affected by the same saturation\neffect of VIs.\nL614: MVC only makes sense for NDVI, not for (downstream) other products (e.g. LAI, fAPAR, GPP, etc.), so be precise on\nwhether you use it only for NDVI or not, and what are you using for the other signals. Also, the reference of MVC should be\nthe original, Holben et at.\nL668: doesn't this window width also depend on the frequency of observations?\nL672: Why MVC?\nL674: you should be more specific, as here what you are doing (with MVC), is to downsample the frequency of actual\nobservations. if the monthly values are effectively a single value in the time series (as with MVC), the time between\nconsecutive observations is not a month, but van vary. So what does that do for CSD?\nL690: How good is the Besnard et al data to differentiate forest age data over a 0.25 by 0.25dd area? it is a requirement that\nthis is very high if this is used in a moving window approach, otherwise this does not provide robust info\n### The following are 3 other specific points regarding the response to reviewers.\n# 1\nFigure 1. Action taken. (2) Background grid lines have also been added to subplots b, d, and e to aid comparison (add the\nLat/Lon in the figures).\nI do not see any background grid lines. Maybe they were lost due to the low-resolution of the figure?\n# 2\nL168-169: “over 40% of paired grids show higher resilience in PF than in NF” this is quite high! If this is the case, I wonder if\nthe statements about management decreasing resilience is founded or whether you should be more conservative.\nConsidering the relevant uncertainties, is the 60%-40% ratio really different than a 50%-50% ratio based on random noise?\nR: Thank you for this important question. Our intention is not only to highlight that certain management activities can\ndiminish resilience, but also to emphasize the complexity and variability of these effects. While forest management activities\nassociated with high human pressure—such as intensive plantations—can significantly reduce forest resilience, there are\nalso management practices characterized by lower anthropogenic pressure that may enhance resilience under specific\nconditions. For instance, as illustrated in Figure 2, the most intensive form of human planting (i.e., plantations) shows the\ngreatest negative impact on resilience, whereas less intensive afforestation practices (i.e., planted forests) have a relatively\nweaker effect on reducing resilience. Action taken. Nonetheless, we recognize that the reader may be misled here, so we\nhave rephrased the paragraph topic, now reading as (Lines 184~186). “Although forest management activities with strong\nhuman pressures have weakened forest resilience, managements with lower anthropogenic pressures can enhance\nresilience under specific conditions (Fig. 2).”\nNot sure if this is answered. My original comment referred to the fact that the percentage of pixels where one or the other\nconclusion could be drawn is potentially indistinguishable from random noise. The authors state that “40% of paired grids\nshow higher resilience in PF than in NF” (L186-188). This means that the main outcome is that 60% of paired grids show\nlower resilience in PF than in NF. Given the many confounding factors affecting natural systems, how confident are you that\nthese results are different than a null hypothesis where 50% of paired grids show higher resilience in PF than in NF, and\n50% show lower resilience in PF than NF? While I agree with the explanation provided by the authors about the\nmanagement effects on forests, this misses to answer the raised concern about the percentages and random noise, which is\na statistical issue.\n# 3\nSupplementary Figures. Given the updated figure specifications (e.g. color scales), try to also update the supplementary\nfigures to be consistent with the figures in the main text. E.g.: the maps in Figure 2 and Figure S4 have now a different color\nbar.", "reviewer_2": "I co-reviewed this manuscript with one of the reviewers who provided the listed reports. This is part of the Nature Communications initiative to facilitate training in peer review and to provide appropriate recognition for Early Career Researchers who co-review manuscripts.", "reviewer_3": "The authors have adequately addressed my primary concern regarding the effect of age on the main conclusions of their\nstudy. Additionally, they have satisfactorily resolved other issues raised in my initial review. After carefully rereading the\nrevised manuscript, I find no further comments or concerns.", "reviewer_4": "Comments for authors\nThe focus of the study on forest management types, biophysical factors and management practices, and forest ecosystem\nresilience” is timely and interesting, as the results will help to inform decision-making on forest management best practices\nfor various forest management types to build resilient forest ecosystems in a changing climate.\nI took the time to review the author response to all the reviewers’ comments, supplementary materials and main document,\nand I think the authors have greatly improved on the manuscript in terms of technical and scientific methodology and clarity\nin the presentation of the results. In particular, the classification of the forest management types by the authors is much more\nimproved and clearer, enabling better comparison between forest management types. The conceptual figure 2(I) that\nillustrates the decreasing vegetation resilience across the different forest management types due to human impact is much\nappreciated. The objectives of the study are clearer and are well answered.\nThe role of forest management practices (e.g., logging rotation cycle, restoration efforts - for managed natural and planted\nforests) and some biophysical factors in shaping the resilience of forest types/ecosystems has been addressed in greater\ndetail in the text and figures, which has significantly strengthen the ecological context and applicability of the findings.\nFurthermore, the authors have provided clearer explanations of the models and their validation procedures, and essentially\nspecified how the models were assessed, including performance metrics, validation techniques, and potential sources of\nbias and uncertainty in the spatial data. Therefore, I agree with the methods used and results reached by the authors.\nHowever, I still have some major comments/suggestions\nMajor comments\nI think that the results are mainly confirmatory of exiting factual knowledge and merely qualitative in the way that they are\npresented in the text. I think that the results could be better quantified, precisely when comparing the resilience between\nforest management types. In particular, instead of using “higher” to qualify unmanaged natural forest’s superior resilience\nover managed forests and planted forest, for example, the authors could quantify the resilience difference in terms of how\nmuch divergence or variation that the study found (resilience difference) that may differ from previous global studies.\nFor example, how many times or by how many folds is one forest management type more resilient than the other?\nThis is an important quantification information that provides more context for comparison than merely using the words “higher” or\n“greater” and “lower” resilience to express the resilience difference between forest management types or forest ecosystems.\nFor example, if an unmanaged natural forest is subject to intensive management practices, how many folds will be the\ndecrease in the resilience factor? This is because, it is factual knowledge that unmanaged NF are older and intact, and thus,\nare more resilient than managed NF and PF in all climatic and biogeographical environments, and what this means in the\ncontext of climate change. Thus, providing the magnitude of the resilience difference in terms of numbers (number of times or\nfolds resilience increased or decreased between forest management types due to forest management practices) should be\none of the research gaps to be filled in the existing literature. As a demonstration, in forest carbon management studies that\ncompare natural and regenerating forests, for example, results are presented in a way that is more captivating such as “old\ngrowth (unmanaged) natural forest sequester 10 to 20 times more carbon than young secondary or regenerating forest\n(managed forest)”; is such quantitative comparison possible with the current study between forest management types or\nforest ecosystems?\nMinor comments\nThe title is still unclear in terms of the structural presentation of the independent and dependent variables and the\ninterrelationship. The biophysical factors and forest management practices determine the resilience of forest management\ntypes. The focus of the study is on the resilience differences of forest management types as influenced by biophysical factors\nand forest management practices? I suggest that the authors should provide a more thoughtful and informative title.\nIn the abstract, natural forest should be used instead of native forest, as this is the terminology used throughout the\nmanuscript; the efficient “supply” or “provision” (choose one term) of ecosystem services should be used instead of\n“functioning”; external disturbances (i.e. ecosystem resilience) is unclear, since external disturbance does not mean the\nsame thing as ecosystem resilience; the former influences the latter.\nThe authors should categorically state in the abstract that unmanaged natural forest demonstrated higher (quantified in terms\nof magnitude of resilience difference) resilience than managed natural and planted forests (or all the other forest types). The\nauthors clearly stated this result in the text as follows “… the results by the nearest neighbor method showed that\nunmanaged NF demonstrate higher resilience than managed”, but not in the abstract. This is the key result and should be\nstated in the abstract first, before stating how these may vary under different biophysical conditions such as wet climate\n(water availability/deficiency modulating the resilience difference effect of forest management types). I think that unmanaged\nnatural forest is more resilient under all climate conditions, though I agree that adapted best management practices can\nimprove resilience under favourable biophysical conditions.\nIn the abstract, “These differences are influenced by biophysical contexts such as local climate, vegetation, and soil”.\nHowever, I think the authors should add one sentence on the limitation of the study, especially in terms of the limited\nvariables included in the model in the abstract.\nThe full stop (.) mark comes after the quotation (”) mark.\nThere are too many abbreviations and acronyms in the text, please reduce them to a minimum to improve on the readability\nof the manuscript. For example, all the forest management types could be written in full in the text.\nThere are still some issues with grammatical accuracy that need greater attention and revision before publication. The\nauthors can seek the help of a native English speaker colleague in the same field of research.\nFigure S1. A clearer font (without text effects) should be used to label the forest management types, a single font should be\nused for all labels in all the figures.\nFigure S12. The shadows in the labels or lettering should be remove.\nFigure S22. The inclination (e.g., slant vertical) of the letterings of forest management types below the figure should be\nimproved for better readability.\nFigure S27. The labels of forest management types below the figure should be restricted to fall outside the figure.\nDoes differences in species composition between forest ecosystems influenced the ecosystem resilience in this study? Was\nthis factor considered as species richness, evenness and diversity under the vegetation variable, in addition to AGB\nvariable? At least, using more vegetation variables in the model would help in detecting the actual influence of biophysical\nfactors on forest ecosystem resilience.\n“This discrepancy is potentially due to the imbalance of samples in the crowdsourced points, where samples of forest\nmanagements like PF are less as compared to NF” – Please provide precise numbers in the manuscript text instead of using\n“less” and “more” qualifiers."}, "author_response": "The authors sincerely thank the anonymous reviewers for their valuable and\nconstructive feedback. We have carefully addressed all concerns and revised the manuscript\naccordingly. Below, we provide detailed responses to each comment from the Editor and reviewers.\nIn the revised manuscript, reviewer comments are shown in black, our revisions are highlighted in\nred, and our responses are indicated in blue. Please note that all line numbers refer to the revised\nversion of the manuscript. Additionally, the references cited within the red-colored text are available\nin the tracked manuscript, while those in the blue responses are included in this document.\nResponse to Reviewer #1\nReviewer#1\n## Overall impression\nThe manuscript attempts an analysis of the sensitivity of forest resilience to human disturbance and\nenvironmental drivers in different forest management categories, arguing that natural forests have\nhigher resilience than managed forests in 60% of the cases. While the subject is of great interest, the\nmanuscript is unfortunately not convincing. The techniques used are inappropriately explained and\nthe limitations of the underlying data are not well exposed. Furthermore, the overall quality and\nstyle of the manuscript is not in-line with what is expected in such a journal. References are\nincomplete, lacking journal names, sometimes pages. There are quite some errors in language that\nshould have been picked up by any spelling corrector. Some phrases seem unfinished (no verb, or\nno object, see L93). Then technically too, various crucial methodological details are not explained\nin the main text, making understanding and following the text quite difficult, and undermining the\ncredibility of the claims that the manuscript tries to put forward. Given these various points, the\noverall recommendation would be rejection. However, if the authors want to improve the manuscript,\nhopefully this review can help them. Below is a list of more general issues followed a more detailed\nlist of specific points.\nR: We appreciate your review and constructive feedback on our manuscript. Your comments have\nprovided useful insights that have guided us in revising the manuscript. We have addressed all your\nspecific concerns as detailed in the following lines. We have incorporated additional analyses to\ninvestigate the impact of varying spatial resolutions (i.e., 500m or 1000m) on our findings. The\nIntroduction and Method Section have been substantially revised to better differentiate between\nmultiple indices within the CSD framework, with particular emphasis on justifying our focus on the\nλAR1, λVar, and AR1. In addition, we added a discussion of sensitivity analysis based on\ncrowdsourced points related to the independence of GFMD and LAI data. The abstract has been\nrestructured to more prominently highlight resilience quantification. We have thoroughly reviewed\nthe reference format and added missing journal names, volume numbers, and page numbers. The\ndetailed reversion can be found in the response and tracked manuscript.\n## General issues\n1. - The abstract gives no clue on the nature of the work. How is resilience quantified? Model\nsimulations? satellite metrics? ground-data?\nR: We appreciate your valuable comments regarding the clue of the abstract and its clarity.\nResilience is quantified based on the Critical Slowing Down (CSD) theory, which posits that as a\nsystem approaches a tipping point, its recovery rate from perturbations slows down. Forest resilience\nis captured through three resilience metrics (AR1, λAR1, and λVar) derived from satellite-based\nvegetation indices (LAI, FPAR, KNDVI, GPP_MODIS, GPP_PML).\nAction Taken. In response, we have made the following revisions to address your concerns: (1)\nDescriptions of CSD Theory, (2) Clarification of Resilience Quantification, and (3) Clarification of\nmethod and datasets from satellite-based metrics in the abstract. Here is the modified Abstract (Lines\n4~22):\n“Policymakers are concerned that newly planted trees may not provide ecosystem services as\neffectively as native forests. The efficient functioning of ecosystem services relies on the capacity\nof forest ecosystems to recover from external disturbances (i.e., ecosystem resilience). This study\nutilizes the critical slowing down (CSD) theory to quantify resilience differences among forest\nmanagement types, particularly between planted forests (PF) and natural forests (NF). CSD\nresilience indicators such as AR1, λAR1, and λVar are estimated from satellite-based vegetation\nindices. This study further investigates the biophysical drivers for the resilience difference between\nPF and NF by machine learning. Results show that most forest management activities weaken forest\nresilience. Forest resilience experiences greater losses under intensive managements However, some\nforest managements can enhance resilience under specific conditions, such as wet climates. PF\ndisplay higher resilience than NF in P/PET (the ratio of precipitation to potential\nevapotranspiration) > 1.5 (wetter climates), while NF is more resilient in P/PET < 1.5 (drier\nclimates). We further detect such resilience differences in all climate zones between PF and NF\ndecreasing from 2001 to 2015. These differences are influenced by biophysical contexts such as\nlocal climate, vegetation, and soil. Regions characterized by cold, wet conditions, dense vegetation,\nand higher soil fertility minimize differences in resilience, suggesting that rational site selection for\nafforestation can enhance forest resilience compared to native forests.”\n2. - Consider that the methodology section comes at the end, so any concept, method, dataset, or\nabbreviation needs to be introduced before being dropped on the reader without context. A few\nexamples are listed in the detailed comments.\nR: We agree that concepts, methods, datasets, and abbreviations should be introduced before being\nmentioned in the text. Action taken. According to the reviewer’s comment, we have introduced key\nmethodological aspects throughout the manuscript and better contextualized their rationale before\npresenting and interpreting results. Here are the modifications:\n(1) Introduction and explanation of paired grid method\nWe have introduced and explained the concept and methodology of paired grids in both the\nIntroduction and Result Section to ensure readers understand the approach:\n\"To differentiate the impacts of human activities from natural variability, we employ a paired\nmethod that compares geographically proximate forest areas under different management\npractices.\" (In the Introductions, Lines 102~104)\n“Specifically, grids are paired using the sliding window method, which moves a fixed-size window\nacross the study area. Within each window, paired grids are selected based on similar forest age\n(differences of ≤10 years), elevation (differences of ≤500 meters), and a minimum of 20% forest\ncoverage for each management type (more details in the Method Section). We point out that the\neffects of potential covariates have been ruled out within the sliding window approach.” (In the\nResults, Lines 159~165)\n(2) Clarification of abbreviations and datasets\nWe have provided the full interpretations of all abbreviations such as variables involved in the\nXGBT model:\n“Abbreviations of biophysical factors include: Temp-temperature, P/PET-the ratio of precipitation\nto potential evapotranspiration, AGB-above-ground biomass, BGB-below-ground biomass, CEC-\ncation exchange capacity, RD-root depth, Sand-sand fraction, and Clay-clay fraction.” (In the\ncaption of Figure 4, Lines 304~307)\n“The eight environmental factors include annual temperature (Temp derived from TeerraClimate),\ndryness index (P/PET derived from TeerraClimate), root depth (RD92), above-ground biomass\n(AGB93), below-ground biomass (BGB93), clay fraction (Clay94), sand fraction (Sand94), and cation\nexchange capacity (CEC94). These factors are selected for their known influence on forest resilience,\nincluding both climatic and land-surface variables.” (In Method, Lines 735~740)\nWe ensured that all dataset abbreviations were introduced with their full names upon first\nmentioning. Additionally, a full-text check was conducted to confirm the proper introduction of all\nabbreviations.\n“Time series of multiple vegetation indices were retrieved for the selected forest pixels over the\nperiod 2001-2015, including leaf area index (LAI), kernel normalized difference vegetation index\n(KNDVI), fraction of absorbed photosynthetically active radiation (FAPAR), gross primary\nproductivity from MODIS (GPP_MODIS), and GPP from the Penman-Monteith-Leuning model\n(GPP_PML).” (Lines 120~124)\n“We derived the reference distribution of forest managements for the year 2000 (Fig. 1a) by\noverlapping tree cover data available for the year 2000 derived from the Global Hansen Forest\nChange Map (GHFC) and the spatial mapping of forest management types retrieved from the Global\nForest Management Datasets (GFMD) for the year 2015.” (Lines 116~119)\nWe mentioned the abbreviations of different management forests in the caption of Figure 2. In\naddition, the details of different management forests are shown in the Method Section and the\nsupplementary.\n“(k) Further comparisons of resilience are performed across various forest management types using\nthe same methodologies. These comparisons include NF vs. PF, NF vs. Plantations, NF WOM vs.\nF WM, NF WOM vs. NF WM, and NF WOM vs. Plantations. The notation “X vs. Y” indicates the\nresilience difference, calculated as the resilience of Y minus that of X. Statistical significance tests\nare performed to assess whether these resilience differences are statistically distinguishable from\nzero. Bars marked with an asterisk (*) denote statistical significance at p < 0.05, while bars labeled\nwith 'NS' denote non-significant differences. (l) A conceptual figure illustrates different forest\nmanagement types (e.g., NF WOM). Abbreviations used in the figure include: NF-Natural Forests,\nPF-Planted Forests, NF WOM-Natural Forests without Managements, NF WM-Natural Forests with\nManagements, and F WM-Forests with Managements.” (In the caption of Figure 2, Lines 206~216)\n(3) Extending resilience metrics concept\nWe extended the introduction and quantification of resilience metrics, clarifying the differences\namong them, particularly emphasizing why the metric of variance is less used. Changes in the\nrevised text are in Lines 88~97 and Lines 621~641.\n3. - In the introduction, four warning signals are presented: AR1, variance, λAR1, and λVar. Yet, the\nresults are based only on three, excluding the variance metric. This needs to be motivated or risk\nfeeling like cherry-picking.\nR: We agree with the reviewer’s comments. It needs to be clear explanation of excluding the\nvariance metric. In the revised manuscript, we have provided a clear justification. Although variance\nserves as a warning signal for potential losses in resilience or stability, its application is limited for\nseveral reasons: (1) Sensitivity to Environmental Noise: Variance is highly influenced by external\nenvironmental factors (i.e., climate change or human activities), which can introduce noise and\nobscure the true system dynamics, particularly in observational datasets1. (2) Lower Sensitivity to\nCritical Transitions: Variance is less effective in detecting critical transitions compared to AR1-\nbased metrics, as it does not capture recovery dynamics as robustly1. A recent study showed that\nAR1 increased, while the variance decreased for global vegetation resilience2 (the signs of variance\nshowed the resilience increasing at a global scale) (Fig. R1). Many other studies have reported that\nglobal vegetation resilience has declined in three decades2-4. (3) Physics interpretation: Recent\nstudies have increasingly focused on the AR1, λAR1, and λVar due to their closer alignment with\nthe dynamic concept of resilience 4-6. Whereas the concept of variance is more closely related to the\nfluctuation of the system. (4) Previous studies: Recently many studies have begun to focus on the\napplication of λAR1 and λVar4,7,8, due to their improvement from variance and AR1. According to\nthe formula, the newly proposed λVar calculation process includes not only variance but also AR1,\nwhich may be a more comprehensive index. λAR1 is derived by applying a logarithmic\ntransformation to AR1, giving it a more direct physical interpretation related to recovery rates.\nFigure R1. (a) Median AR1 for vegetation optical depth (VOD) data (blue; median taken over all\nAR1 series globally). (c) Same as (a) but for the variance. The blue line is the VOD datasets from\nthe satellite, while the red line indicates results generated by satellite-based data with random noise\n(Signal-to-Noise Ratio (SNR) = 0.1). This figure is derived from a previous study2.\nSmith et.al applied the λAR1 and λVar, which indicate the rate of recovery from disturbances and\nprovide deeper insights into system resilience4. AR1 (autocorrelation) indicators mainly capture a\nsystem’s “memory effects” or delays in returning to equilibrium, while variance focuses on the\noverall magnitude of deviations from the mean1. In contrast, λAR1 and λVar quantify the post-\ndisturbance recovery rate, directly characterizing how quickly the system returns from a non-\nequilibrium state to a stable equilibrium4. In other words, they provide not just a statistical proxy\nfor resilience, but a measure more closely aligned with its true ecological or dynamical meaning\n(the recovery rate after disturbance)3. Furthermore, variance is easily influenced by external\nenvironmental noise, making it challenging to distinguish genuine critical transition signals from\nrandom fluctuations9. By directly indicating recovery rates, λAR1 and λVar help isolate system-\nintrinsic dynamics and true recovery properties, rather than being confounded by noise-induced\ntemporal correlations and variance changes7. While an increase in AR1 is often regarded as an early\nwarning signal of an impending critical transition, studies have shown that the variance indicator\nmay not display similarly pronounced changes in certain contexts2. The direct measurement of\nrecovery rates via λAR1 and λVar allows for the detection of subtle declines in resilience more\nsensitively, providing clearer and more robust evidence of resilience deterioration even when\nvariance does not exhibit a noticeable increase or decrease7. Additionally, our analysis based on\nλVar has demonstrated similar patterns to results by λAR1, indicating that the choice of indicators\ndoes not significantly affect our results. Furthermore, from a formula perspective, λVar incorporates\nboth AR1 and variance information in the quantification process, making it a more comprehensive\nindex. In summary, our study has excluded variance metric, instead focusing on AR1, λAR1, and\nλVar—with our analysis by the three metrics showing similar patterns (see Supplementary Fig. S3-\n\nS6 and Supplementary Fig. S17–S21).\nAction taken. In response to your comments, we modified some texts in the Introduction to clarify\nthe reasons for using three metrics but exclude the variance (briefly introduce). Here are the\nmodified texts: (Lines 88~97)\n“The AR1 measures how strongly the current state of a system is influenced by its previous state18,\nand the variance measures the magnitude of fluctuations around the system's equilibrium43. In\naddition, new metrics such as λAR1 and λVar have been developed, which are more closely related\nto the rate of vegetation recovery from disturbances, further aligning with the concept of resilience\nin dynamic systems19, 44, 45. Compared to the other three metrics, variance is susceptibly affected by\nenvironmental noise, and thus is challenging to detect critical transitions. As a result, recent studies\nhave increasingly focused on AR1, λAR1, and λVar due to their stronger alignment with the concept\nof forest resilience45, 46, 47.”\nMoreover, we have specifically introduced and explained it in the Method Section. (Lines 621~637)\n“Based on CSD theory16, various indicators, including lag-1 autoregression (AR1) and variance,\nλAR1, and λVar, have been widely employed to quantify ecosystem resilience19, 44. However,\ncompared to other indicators, variance is highly susceptible to environmental noise, leading to less\nreliability in detecting critical tipping points43. For instance, while global vegetation resilience is\nreported to be in decline, as reflected by significantly increasing AR1 values and abrupt\ndisturbances18, variance has instead shown a decreasing trajectory68. As a result, recent studies have\nincreasingly favored AR1 as an alternative to variance, as it can detect resilience and its changes\nacross various ecosystems, including forests, wetlands, and grasslands, before regime shifts occur39,\n44, 45. Additionally, in light of the influence of observational space-temporal resolution and\nobservational noise, complementary indicators such as λAR1 and λVar have been introduced45.\nThese metrics quantify recovery rates, measuring how quickly the system returns from a non-\nequilibrium state to a stable equilibrium85. In other words, they provide not only a statistical proxy\nfor resilience but also a measure more closely aligned with its true ecological state. Consequently,\nwe evaluated forest resilience using AR1, λAR1, and λVar.”\n4.- Check that the references to the figures are correct throughout the document. It felt that some\nwere not, which made following along sometimes very confusing. Also, Figures and supplementary\nFigures should be always presented in order, and all should be referenced in the text.\nR: Action taken. According to the reviewer’s comment, we have made the following improvements\nin references: (1) Figure and Supplementary Figure References: We have ensured that all figures\nand supplementary figures are correctly referenced in the text and are presented in the proper order.\nAny discrepancies or missing references have been corrected, ensuring clarity and consistency\nthroughout the document (see References Section). (2) References in the Reference List: We have\nchecked the entire reference list for completeness and accuracy. Missing details, such as journal\nnames, have been added.\n5.- Statements based on studies from the literature should probably be moved all to the discussion,\nas in the results it is not always clear if the statements made are based on own results or the literature.\nAlternatively, it should be made very clear whenever a reference is cited as basis for a statement, or\nas support of own conclusions.\nR: We agree that it is essential to distinguish between statements derived from our results and those\nbased on the literature. Action taken. To address this issue, we have taken the following steps:\n(1) Most literature-based statements previously included in the Result Section have been moved to\nthe Discussion Section, where they are more appropriately contextualized. Here are the modified\ntexts, now reading as: (Lines 382~394)\n“Human activity datasets, generated by road maps, nightlight data, and forest cover change datasets,\nhave been used to delineate intact and managed forests in the previous study. However, a notable\nambiguity arises from the definition of managed forests, as these studies and datasets did not further\ndifferentiate between various types of forest management practices49. Despite this limitation, both\nthe previous study18 and our findings illustrate the damage of human management practices on\nvegetation resilience (Fig. 2). Managed forest fragments experience human disturbance and pressure,\nincluding forest planting, rotation, fragmentation, and irrigation. These factors directly influence\nvarious aspects such as tree species composition, forest age, and physiological processes50,\ndisrupting forest recovery processes. Furthermore, these management practices also affect how\nvegetation responds to climate change, ultimately diminishing the resilience of forest ecosystems in\nrecovering from disturbances18, 51.”\n(2) For statements in the Result Section, we have ensured it is explicitly stated whether they are\nbased on our findings or supported by the literature. Proper citations and clear connectors such as\n\"consistent with\" have been added for clarity. Here are the modified texts:\n“The CSD theory indicates that systems with lower resilience are more sensitive to external\ndisturbance42.” (Lines 309~310)\n“This analysis reveals a general decline in the resilience of global forests, which is consistent with\nprevious studies18, 48.” (Lines 348~349)\n6.- The definition of forest types seems pretty central to the research question, so it should also be\ncentral in the text, but it is now only relegated to the supplementary materials. Either tables S1\nand/or S2 should be moved to the main manuscript, or the most important forest type definitions\nshould be clearly stated in the text.\nR: We agree with your comment that these definitions are critical to the study and should be more\nprominently presented in the main text. Action taken. To address this, we have moved the Table S1\ninto the manuscript (Table 1), and briefly introduced the specific forest management types. Here are\nthe briefly introduced the specific forest management types: (Lines 541~551)\n“Agroforestry is a forest ecosystem that integrates trees with crops. Oil palms are tropical\nmonoculture plantations for oil production, identifiable by their distinctive crown shapes. Shorter-\nrotation planted forests are established with fast-growing tree species for rapid biomass production,\ntypically involving rotations of up to 15 years aimed at producing timber, pulp, or bioenergy.\nLonger-rotation planted forests are planted for rotations exceeding 15 years, focusing on ecological\nrestoration, carbon sequestration, and biodiversity enhancement. Natural forests with human\nmanagements are subject to human interventions such as logging, trimming, and cutting, while\nnatural forests without human disturbance are pristine ecosystems (e.g. Amazon Forest) that\npreserve their original biodiversity and ecological processes. The details for specific forest\nmanagement types can be found in Table S11.”\nHere are the introduced of Table 1: (Lines 562~571)\n“We further re-classified forest types based on a classification of forest management types to\ncompare resilience differences (Table 1). For instance, we compared resilience between natural\nforests (NF) and Plantations. Plantations are defined as all human-planted forests including planted\nforests (i.e., with longer rotating years and mainly for ecological restoration), short-rotation planted\nforests (up to 15 years), agroforestry, and oil palms. Reclassified forest types include NF WOM\n(natural forests without managements), NF WM (natural forests with managements), F WM (forests\nwith managements), Plantations, PF and NF (Table 1). This classification allowed a detailed\ncomparison of resilience across varying levels of human interventions.”\nTable 1. The definition and description of abbreviations for forest managements\nAbbreviation Definition Description\nIt concludes the natural\nforests without management\nNF Natural forests\nand natural forests with\nmanagement.\nIt concludes the planted\nPF Planted forests forests (longer rotating time\nover 15).\nNatural forests without It concludes the natural\nNF WOM management/ Unmanaged forests without management\nnatural forests (original or pristine forests)\nNatural forests with It concludes the natural\nNF WM management/ Managed forests with management such\nnatural forests as logging and clearing up.\nF WM Forests with It concludes the planted\nmanagement/ Managed forests, planted forests with\nforests rotating up to 15 years,\nagroforestry, and oil palms\nIt concludes the planted\nEconomic Plantations/\nforests, planted forests with\nPlantations Board definition of planted\nrotating up to 15 years,\nforests\nagroforestry, and oil palms\n7.- The arguments about the independence between forest management and resilience of the GFMD\ndata should be more central in the discussion, as it is potentially a huge underlying limitation.\nTherefore, the independent tests on crowdsourcing data should be highlighted more.\nR: We agree that this is an important aspect to involve in the Discussion Section. Action taken. We\nhave added a new paragraph describing analysis based on crowdsourced data, and discussed the\nconsistency and inconsistency with results by paired grids. Here are modified texts, now reading as:\n(Lines 469~491)\n“The GFMD map generated by factors like LAI may affect the independence of our findings. This\nstudy extracts the resilience of crowdsourced points overlapped with the GHFC map. Results show\nthat intensive human managements reduce resilience, aligning with our initial GFMD-based\nfindings. For example, PF with rotations up to 15 years (i.e., wood plantations/shorter-rotation PF)\nexhibit lower resilience compared to those with rotation periods over 15 years (i.e., planted forests\nfor restoration/longer-rotation PF) (Supplementary Fig. S26). These shorter-rotation PF are\neconomically composed of a single species and are often partially harvested at regular intervals66.\nAs a result, species richness and biodiversity in shorter-rotation PF are generally lower compared\nto longer-rotation PF, reducing the capacity to recover from disturbances52. Our study also supports\nthe indications that the increase in tree species richness can further enhance forest resilience\n(Supplementary Fig. S26). It encourages multi-species plantations as a policy recommendation for\nthe enhancement of forest resilience for afforestation. Additionally, the results by the nearest\nneighbor method showed that unmanaged NF demonstrate higher resilience than managed NF\n(Supplementary Fig. S27), aligning with our initial GFMD-based findings. However, some findings\ndiffer from GFMD-based outputs. For example, the resilience difference between PF and NF\nremains below 0 until P/PET < 2.5 based on the crowdsourced datasets (Supplementary Fig. S28),\nwhile the GFMD-based findings in Figure 3 show this threshold occurs until P/PET < 1.5. This\ndiscrepancy is potentially due to the imbalance of samples in the crowdsourced points, where\nsamples of forest managements like PF are less as compared to NF. Overall, caution is needed when\ninterpreting specific thresholds and effects.”\n8.## Specific points\nR: We sincerely appreciate your detailed comments and suggestions. We have addressed the relevant\npoints one by one in our responses table and texts. To ensure clarity, we have categorized the\nresponses into two parts: Part 1: the comments and suggestions about the figures, and Part 2: the\ncomments and suggestions about the texts.\nPart 1: the comments and suggestions about the figures\nThe comments and suggestion about the figures:\nFig 1: confusing and not vizually clear. strangely structured. Axis titles are missing. Lat/lon\nindicators are missing. So far, LAI has never been mentioned in the text, so the subplots b and d\ncome as a bit of a surprise. Background grid lines could aid the comparison of subplots b, d, e., and\nthe x-axes could be aligned for more consistency and ease of interpretation. In the figure description:\n“as an example of a specific global grid point”, you could say what point you are presenting (e.g.\ncentral coordinates). Finally, the last part of the figure caption probably belongs to the main text:\n“The other examples in Central Africa… …found in Figure S2, showing similar spatial patterns”.\nR: We have corrected Figure 1 in the revised version of the manuscript according to the reviewer’s\nsuggestion. Action taken. Below are our responses and the corresponding changes made: (1) We\nhave restructured the figure layout for better organization such as moving the color label to the\nbottom of the figure. (2) Background grid lines have also been added to subplots b, d, and e to aid\ncomparison (add the Lat/Lon in the figures). (3) The x-axes have been aligned for more consistency\nfor b, d, and e. We added the specific global grid point in the figure for explanation (-88.877, 31.905).\n(4) We added the introduction of LAI and resilience metrics (Lines 116~128 and Lines 129~135).\nWe have revised the main text to introduce LAI earlier and explain its relevance to the analysis. We\nalso introduce the resilience increases with the decrease of AR1/λAR1/λVar. (5) We agree that the\nlast part of the figure caption belongs to the main text (Lines 139~142). This content has been\nrelocated to the main text, and the figure caption has been revised according. Here are the modified\ntexts:\n“We derived the reference distribution of forest managements for the year 2000 (Fig. 1a) by\noverlapping tree cover data available for the year 2000 derived from the Global Hansen Forest\nChange Map (GHFC) and the spatial mapping of forest management types retrieved from the Global\nForest Management Datasets (GFMD) for the year 2015. Time series of multiple vegetation indices\nwere retrieved for the selected forest pixels over the period 2001-2015, including leaf area index\n(LAI), kernel normalized difference vegetation index (KNDVI), fraction of absorbed\nphotosynthetically active radiation (FAPAR), gross primary productivity from MODIS\n(GPP_MODIS), and GPP from the Penman-Monteith-Leuning model (GPP_PML). Seasonal period\npatterns and long-term trends were removed from the original vegetation time series to make them\nstationary (Fig. 1b, d) and CSD indicators were retrieved on the residual signal to characterize the\nforest resilience. In this study, we express resilience in terms of AR1, λAR1, and λVar (additional\ndetails in the Method Section).” (Lines 116~128)\n“Results show that areas characterized by human activity, as reflected in the forest management\nmap, generally have lower resilience than natural forest without management (NF WOM) (Fig. 1).\nThis emerges particularly evident for instance in the Amazon and Central African tropical forests\nwhere the AR1 estimated from LAI (AR1-LAI) in natural forests with management (NF WM) is\nsubstantially higher than that one quantified in NF WOM (Fig. 1j, k and Supplementary Fig. S1).\nThe higher AR1/λAR1/λVar indicates lower forest resilience.” (Lines 129~135)\n“Resilience estimated by using other vegetation indices (KNDVI, FAPAR, GPP_MODIS,\nGPP_PML) and alternative CSD indicators (λAR1, λVar) exhibits similar spatial patterns\n(Supplementary Fig. S2) and confirm the robustness of our assessment.” (Lines 139~142)\nFig. 1 Spatial Distribution of Global Forest Managements and Forest Resilience. (a) The spatial\ndistribution of global forest management data (GFMD) in 2000 is generated by overlaying GFMD\nin 2015 onto GHFC in 2000 (details in the Method Section). (c) The spatial distribution of forest\nresilience displayed here is calculated using LAI and AR1 from 2001 to 2015 (AR1-LAI). (b-d) The\noriginal LAI time series is processed with detrending and deseasonalization methods to calculate\nresilience, using the specific global grid point (-88.877, 31.905) as an example. (e) Vegetation\nresilience in this study is estimated by AR1, λAR1, and λVar, showing similar changes at the\nspecified grid point. The higher values of AR1, λAR1, and λVar indicate lower forest resilience. (f–\nk) Forest management activities are observed to reduce vegetation resilience in geographically\nsimilar regions, such as (f, g) the Southeastern United States, (h, i) South Asia, and (j, k) the Amazon\nRainforest.\nFigure 2: the color palettes chosen to represent the resilience difference and the datasets are very\nsimilar (browns and greens especially). This makes the interpretation of the figure difficult, so the\nso the\ncolors should be differentiated more. Alternatively, subplots b, d, and f should come with their own\ncolor legend. In the figure description (here and in other figures) you mention a significance test,\nbut it is not specified what was tested. The significant difference from 0?\nR: We have corrected Figure 2 in the revised version of the manuscript according to the reviewer’s\nsuggestion. Action taken. Below are our responses and the corresponding changes made: (1) We\nagree that the original color palettes for resilience differences and datasets were visually similar,\nmaking interpretation challenging. To address this, we have updated the color palettes (now the\ncolor palette is blue and red) in Figure 2. (2) In addition, we have added separate legends for b, d,\nand f subplots for reference. (3) Thank you for pointing out the need for clarification regarding the\nsignificance test. In the revised figure description, we have specified that the significance test\nevaluates whether the differences are significantly different from 0. Here are the modified texts and\nFigure 2:\nFig. 2 Vegetation Resilience Differences Under Forest managements. (a-f) Spatial patterns of\nresilience difference between PF and NF are estimated using LAI and three measurements: (a, b)\nAR1, (c, d) λVar, and (e, f) λAR1. (b, d, f) These resilience differences are also assessed using\nmultiple vegetation indices—GPP_MODIS (blue), GPP_PML (green), LAI (orange), KNDVI\n(blue-black), and FPAR (earthy color)—which display similar patterns across the latitudinal\ngradient. (g-j) Paired grids in key regions including the Southeastern U.S., Northwestern U.S.,\nEurope, and East Asia are highlighted to showcase regional variations. (k) Further comparisons of\nresilience are performed across various forest management types using the same methodologies.\nThese comparisons include NF vs. PF, NF vs. Plantations, NF WOM vs. F WM, NF WOM vs. NF\nWM, and NF WOM vs. Plantations. The notation “X vs. Y” indicates the resilience difference,\ncalculated as the resilience of Y minus that of X. Statistical significance tests are performed to assess\nwhether these resilience differences are statistically distinguishable from zero. Bars marked with an\nasterisk (*) denote statistical significance at p < 0.05, while bars labeled with 'NS' denote non-\nsignificant differences. (l) A conceptual figure illustrates different forest management types (e.g.,\nNF WOM). Abbreviations used in the figure include: NF-Natural Forests, PF-Planted Forests, NF\nWOM-Natural Forests without Managements, NF WM-Natural Forests with Managements, and F\nWM-Forests with Managements.\nFigure 4: I am a bit confused about SHAP values vs feature values. What do they represent and why\ndo we need both? Are the feature values similar to effect coefficient? It should be made clearer even\nfor readers that are not very familiar with the methodology.\nR: Thank you for your question regarding SHAP values and feature values. (1) SHAP values\nrepresent the contribution of each feature to the model’s prediction for a specific observation (i.e.,\nthe ranking of importance). (2) They also quantify how much a particular feature value increases or\ndecreases the predicted resilience difference. Positive SHAP values indicate a positive contribution,\nwhile negative SHAP values indicate a negative contribution. Feature values are the actual input\nvalues (e.g., temperature, root depth) for a given feature in the dataset. They represent the observed\nor measured characteristics of the biophysical or forest management factors in the study. While\nSHAP values explain the contribution of features to the model prediction, feature values provide\ncontext for interpreting these contributions. For instance, higher temperature with SHAP values > 0\nindicates its positive impact on resilience differences, but the feature value (e.g., a specific\ntemperature range like higher temperature values) helps us understand the conditions under which\nthis impact occurs. Thus, high temperatures have positive SHAP values and low temperatures have\nnegative SHAP values, indicating that as the temperature increases, ΔλAR1 between PF and NF\nwill also increase.\nAction taken. We agree that providing a clearer explanation is essential for readers who may not\nbe familiar with the methodology. Below, we clarify their roles in the Result and Method Section:\n“The ranking of Shapley Additive Explanations (SHAP) values suggests that climate variables, such\nas temperature and P/PET, are more important than other factors like biomass and root depth (i.e.,\nrelated to the vegetation context) (Fig. 4b). The SHAP value reflects the effects of variables on the\npredicted resilience difference (e.g., positive or negative). While the feature value interprets whether\na feature is at a level (e.g., high temperature) that amplifies or diminishes resilience differences.\nFor instance, negative SHAP values associated with high belowground biomass indicate that high\nbelowground biomass shows negative effects on resilience differences (Fig. 4b).” (Lines 280~287\nin Result Section)\n“Specifically, SHAP values show how environmental variables (e.g., temperature and root depth)\ninfluence the predicted resilience difference. Positive values indicate the positive effects, and vice\nversa. Feature values, like observed temperatures, interpret the effects under specific conditions.\nThe direction and magnitude of the variable effects depend not only on the SHAP value itself but\nalso on the value of the feature. A high feature value with a positive SHAP value indicates that as\nthe feature rises to a higher value, it contributes positively.” (Lines 720~727 in Method Section)\nWe also modified the related descriptions of SHAP values in the figure’s caption.\nFig. 4 Analyzing the Impact of Biophysical Factors on Vegetation Resilience Differences\nBetween PF and NF. (a, c) Observed versus predicted ΔλAR1-LAI are presented for both\n(PF, NF)\ntraining and testing samples. The coefficient of determination (R²) and Kling-Gupta Efficiency\n(KGE) for each sample are indicated in the labels, with the frequency distribution depicted by the\ncolor gradient. (b) The SHAP values for environmental factors affecting the resilience difference\nare displayed. In this panel, red and blue points represent high and low feature values, respectively.\nSHAP values quantify the contribution of each feature to the model's prediction (e.g. positive or\nnegative). For example, red points on the SHAP value plot that are greater than zero indicate the\nenhancing effects produced under high values. The “Location_Factor” indicates the local contexts\nof paired grid cells, while the “Delta_Factor” represents the difference in biophysical factors\nbetween PF and NF. Abbreviations of biophysical factors include: Temp-temperature, P/PET-the\nratio of precipitation to potential evapotranspiration, AGB-above-ground biomass, BGB-below-\nground biomass, CEC-cation exchange capacity, RD-root depth, Sand-sand fraction, and Clay-clay\nfraction.\nFigure 5: Results descriptions and references to supplementary materials should probably go in the\nmain text (results) and not in the figure description. Clarify what the significance test is testing. Also,\n“Vegetation greenness is limited by water in dry conditions and by energy in wet conditions, which\nmay be related to shifts in resilience differences between PF and NF”: to me, this sounds related to\nactual resource limitations rather than resilience differences.\nR: Here, we respond to the reviewer’s comment by addressing separately the two sections.\nSection1: the comments on the caption of Figure 5\nWe have corrected Figure 5 in the revised version of the manuscript according to the reviewer’s\nsuggestion. Action taken. We agree that result descriptions and references to Supplementary\nMaterials should be placed in the main text rather than in the figure caption (Lines 330~332).\nFurthermore, we have clarified the significance test in the caption, with the test of the significance\nbetween resilience differences and 0. We have updated the figure description to focus solely on the\nfigure's content. Here are the modifications:\n“The conceptual diagram illustrates the reversed resilience patterns of PF and NF, driven by\nvegetation sensitivity to climate factors (Supplementary Fig. S12).” (Lines 330~332)\n“Fig. 5 Climate Sensitivity of Vegetation for PF and NF. (a, b) Spatial patterns of vegetation\nsensitivity to (a) SM and (b) PET for PF and NF are displayed. (c, d) Vegetation sensitivity is\nexamined along a gradient of P/PET, emphasizing the variations in response across different climatic\nconditions. (e–j) Differences in climate sensitivity between PF and NF are quantified using five\nvegetation indices. Statistical tests are conducted to evaluate if these differences are statistically\ndistinct from zero. Bars marked with an asterisk (*) denote significant differences (p < 0.05), while\nthose labeled with \"NS\" indicate non-significant results. The climate variables analyzed include\nboth water-related and energy-related factors: precipitation (Pre), soil moisture (SM), Palmer\nDrought Severity Index (PDSI), potential evapotranspiration (PET), and temperature (Temp).\nStandardization across these variables and vegetation indices ensures comparability of forest\nsensitivity to these factors.” (Lines 334~344)\nSection 2: the relationship between vegetation sensitivity and resilience.\nThanks for your insightful comment. We agree that the statement about vegetation greenness being\nlimited by water in dry conditions and by energy in wet conditions does indeed reflect resource\nlimitations, rather than directly addressing resilience differences between PF and NF. However, our\nintention was to highlight how these resource constraints might influence the resilience of the\necosystems. This difference in sensitivity to resource availability could help explain the observed\nresilience differences between PF and NF. (1) Vegetation sensitivity refers to the degree to which\nvegetation greenness responds to changes in external conditions, such as precipitation. For example,\nhigh sensitivity to precipitation implies that a reduction in rainfall or the occurrence of drought\nwould lead to a significant decline in vegetation greenness. High sensitivity may indicate that\nvegetation greenness is more susceptible to external climate changes or sudden climatic events,\nresulting in greater fluctuations and lower stability (weaker resilience). For instance, PF are more\nsensitive to water factors like precipitation than natural forests. When external disturbances occur\n(such as drought events), the biomass reduction in PF is greater than in NF. Therefore, plantations\nrequire a longer time to recover to their initial state, resulting in lower resilience compared to NF.\n(2) Additionally, existing studies have demonstrated that vegetation sensitivity can also be used to\npredict tipping points10. Sensitivity is considered a key functional trait of vegetation, and changes\nin sensitivity are often accompanied by shifts in resilience indicators, such as AR1. Critical Slowing\nDown (CSD) theory11 indicates that as a system approaches a tipping point (TP), its sensitivity to\nexternal disturbances increases. This means that near the critical point, even small perturbations can\ncause significant changes or large fluctuations in the system (i.e., Because the sensitivity of the\nsystem is very high, it’s state can easily shift), a phenomenon known as \"critical slowing down.\"\nIn this context, the recovery time of the system is prolonged, indicating a decrease in the recovery rate.\nThe reduction in recovery rate reflects the system's slowed pace in returning to its original stable\nstate, making it more vulnerable to external changes.\nAction taken. In response, we have modified texts in discussion to explain the relationship between\nsensitivity and resilience.\n“Existing studies have indicated the relation between resilience and climate sensitivity16, 60. For\ninstance, forests with higher sensitivity to water variability experience greater declines in biomass\nwhen subjected to external disturbances like droughts, prolonging the recovery process and making\nit more challenging for forests to recover their original state.” (Lines 427~431)\nL527: why are you limited to 1km by computational constraints? it seem very likely it could all be\ndone at 500m\nR: Action taken. Thanks for your important comments. In response, we tested the main conclusions\nat different resolutions such as 500m and 1000m. The results show that the resolution has a minor\nimpact on our main conclusions. For example, the resilience estimates for parts of the African forest\nby 500m and 1km resolutions exhibit similar spatial patterns. The conclusions at 500m resolution\nalso indicate that when P/PET < 1.5, the resilience of NF (Natural Forest) is greater than PF (Planted\nForest); whereas when P/PET > 1.5, the resilience of PF is greater than NF. These findings are also\naligned with the results by 1km. Overall, the spatial resolution of analysis datasets may not\nsignificantly influence our main findings (Fig .R2).\nFig R2. The results by 500m and 1000m spatial resolution. (a-d) The spatial distribution of λAR1\nfor the African forests by 500m and 1000m. (e, f) The resilience difference between PF and NF\n(λAR1 and λVar) by 500m and 1000m. (g, h) The resilience difference between PF and NF across\nthe P/PET gradient.\nPart 2: the comments and suggestions about the texts\nL17: P/PEF should be defined\nR: Action taken. We have added the definition of P/PET in the 15-17 lines of revised manuscript:\n“PF display higher resilience than NF in P/PET (the ratio of precipitation to potential\nevapotranspiration) > 1.5 (wetter climates), while NF is more resilient in P/PET < 1.5 (drier\nclimates).”\nL49: correct “bring” to “bringing”.\nR: Action taken. We rewrite this sentence. The updated sentence reads (Lines 48~50): “This is\nparticularly critical in view of the ongoing intensification of natural and anthropogenic disturbance\nregimes that could affect the long-term resilience of global forests.”\nL60: “However, despite models appearing to predict stably, forest dynamics in observation remain\ncomplex and diverse” unclear sentence.\nR: Action taken. We have corrected the text in the revised version of the manuscript according to\nthe reviewer’s suggestion. The revised sentence now states (Lines 58~63): “Although these models\nare widely used for large-scale impact assessments, they only partially capture the complexity of\nforest dynamics, particularly the ecosystem's response to evolving human factors. Real-world\nforests may indeed undergo sudden and abrupt changes from a range of unexpected disturbances\n(e.g. forest clear-up by humans) that exceed model forecasts27, 28.”\nL62: correct “plating” to “planting”.\nR: Action taken. We have corrected the text in the revised version of the manuscript according to\nthe reviewer’s suggestion. The updated sentence reads as (Lines 64~66): “activities like planting\nnon-native species in arid zones not only resulted in unexpectedly lower tree survival rates but also\nled to a serious depletion of water resources29.”\nL64: remove double comma after “Furthermore”; clarify what is meant with “early stable\npredictions”.\nR: Action taken. We have addressed the issue by removing the double comma after “Furthermore”\nand clarifying the meaning of \"early stable predictions\" (i.e., the predictions of stable growth in\nAmazon forests from ecological models) for better readability. The updated sentence now reads as\n(Lines 66~67): “Furthermore, despite predictions of stable growth from ecological models, signs of\ndeath or degradation have been observed in the Amazon rainforest30.”\nL67: “Recent study report that…” correct singular/plural form (either “A recent study reports that…”\nor “Recent studies report that…”.\nL68: correct “dectect” to “detect”\nR: Action taken. We have corrected the text in the revised version of the manuscript according to\nthe reviewer’s suggestion. The updated sentence reads as (Lines 67~69): “A recent study reports\nthat land surface models (i.e., Dynamic Global Vegetation Models) rarely detect the resilience trend\nin carbon sinks31.”\nL69-70: “…it will be better to…” this formulation sounds very strong. Both observational and\nmodelling frameworks have their uses. I would present them as alternatives to one another, rather\nthan stating that one is superior.\nR: Thank you for your suggestion. We agree that both observational and modeling frameworks have\ntheir respective strengths and should be presented as complementary approaches rather than one\nbeing superior to the other. Action taken. The updated sentence now reads as (Lines 74~76):\n“Therefore, given the parameter and structure uncertainty of models, observational datasets offer an\nalternative to assess the impact of forest managements on the resilience of forests35.”\nL84-86: it might be worth spending a few more words/sentences to introduce and define the\nconcepts of AR1, variance, λAR1, and λVar, as they seem rather central to the presented results.\nR: Action taken. We agree that these resilience metrics need to be introduced. We have corrected\nthe text in the revised version of the manuscript according to the reviewer’s suggestion. The updated\nsentence now reads as (Lines 88~97): “The AR1 measures how strongly the current state of a system\nis influenced by its previous state18, and the variance measures the magnitude of fluctuations around\nthe system's equilibrium43. In addition, new metrics such as λAR1 and λVar have been developed,\nwhich are more closely related to the rate of vegetation recovery from disturbances, further aligning\nwith the concept of resilience in dynamic systems19, 44, 45. Compared to the other three metrics,\nvariance is susceptibly affected by environmental noise, and thus is challenging to detect critical\ntransitions. As a result, recent studies have increasingly focused on AR1, λAR1, and λVar due to\ntheir stronger alignment with the concept of forest resilience45, 46, 47.”\nL90-93: “Motivated by the issue of what impacts of human activities especially forest management\nincluding PF and NF on vegetation resilience, this study used the latest global forest management\ndataset, multi-source satellite observations, and climate datasets to investigate. Here, this study\nmainly to:” unclear sentences.\nL94: it might be worth introducing the “paired method” before mentioning it in the research\nquestions.\nR: Action taken. We have corrected the text in the revised version of the manuscript according to\nthe reviewer’s suggestion. We re-construct the sentence for clarification and introduce the paired\ngrids method before the research questions. The revised version of the sentence is as follows (Lines\n98~106): “In this study, we investigate how different forest management practices (e.g., PF, NF)\ninteract with biophysical factors to affect global forest resilience. To achieve this, we integrate\ncomprehensive global forest management datasets, multi-source satellite observations, and climate\nproducts. Forest resilience is expressed in terms of AR1, λAR1, and λVar from Critical Slowing\nDown (CSD) theory. To differentiate the impacts of human activities from natural variability, we\nemploy a paired method that compares geographically proximate forest areas under different\nmanagement practices. Machine-learning techniques are used to quantify the impacts of biophysical\nfactors on resilience differences among forest management types.”\nL106: “In Spatial”, meaning?\nL107: “decline” as per meaning become smaller, fewer, weaker, etc.\n(https://www.oxfordlearnersdictionaries.com/definition/english/decline_1) is an intransitive verb,\nmeaning that it cannot have an object. In your formulation: “…human forest management activities\nhave declined the resilience” you want to use a transitive verb such as “decreased”. Please check\nthe use of decline vs decrease throughout the manuscript.\nR: Action taken. The phrase \"In Spatial\" was unclear and unnecessary (in the original texts, we\nwant to express “the results are based on spatial analysis”). We have removed it. We agree that\n\"decline\" is an intransitive verb and cannot take an object. Furthermore, we carefully reviewed the\nmanuscript to change “decline” to “decrease”/”reduce” or change the “decline” to “…lower than…”.\nThe revised sentence now reads (Lines 129~135):\n“Results show that areas characterized by human activity, as reflected in the forest management\nmap, generally have lower resilience than natural forest without management (NF WOM) (Fig. 1).\nThis emerges particularly evident for instance in the Amazon and Central African tropical forests\nwhere the AR1 estimated from LAI (AR1-LAI) in natural forests with management (NF WM) is\nsubstantially higher than that one quantified in NF WOM (Fig. 1j, k and Supplementary Fig. S1).\nThe higher AR1/λAR1/λVar indicates lower forest resilience.”\nL108: Explain what the GHFC and GFMD are.\nR: Action taken. We have added following explanations in the 116~119 lines: “We derived the\nreference distribution of forest managements for the year 2000 (Fig. 1a) by overlapping tree cover\ndata available for the year 2000 derived from the Global Hansen Forest Change Map (GHFC) and\nthe spatial mapping of forest management types retrieved from the Global Forest Management\nDatasets (GFMD) for the year 2015.”\nL111: Rationale for these metrics is needed. We have no explanation of what they do, how they are\nconstructed and what assumptions hide behind. This needs to be clear in the main text, not guessed\nor ignored. Also, in the introduction you talked about four warning signals: AR1, variance, λAR1,\nand λVar. Why are the results based only on three, excluding the variance metric?\nR: The analysis excludes the variance metric for the following reasons: Variance is highly sensitive\nto external noise, such as climate variability and human activities, which can introduce bias.\nDuring critical transitions, variance can sometimes decrease instead of increasing, making it less\nreliable for detecting system state changes. Variance is more closely associated with system\nfluctuations, while AR1, λAR1, and λVar are more aligned with the dynamic concept of resilience\nand recovery rates in complex systems. Recently many studies have begun to focus on the\napplication of λAR1 and λVar4,7,8, due to their improvement from variance and AR1. Action taken.\nWe have revised the manuscript to include detailed explanations and comparisons of AR1, variance,\nλAR1, and λVar in the Introduction Section (Lines 88~97).\nL112: is the starting time of the analyzed data 2000 or 2001? This seems inconsistent.\nR: The analysis starting from 2001 is based on the following reasons: (1) Sliding Window Analysis:\nMODIS Terra data in 2000 do not cover a complete calendar year, as observations began mid-year.\nThis incomplete temporal coverage could affect the consistency of sliding window analysis and\nintroduce biases in estimating seasonal variations, as different months would have unequal sample\nsizes. (2) Forest Cover Data Constraints: The Hansen Forest dataset provides forest cover data only\nfor the year 2000. Estimation of forest cover has been applied more on Hansen’s data and has a\nhigher precision12. Considering these factors, we assume that the differences in long-term resilience\nestimation between 2001–2015 and 2000–2015 are negligible, and thus we chose to start the analysis\nfrom 2001 instead of changing the forest cover dataset.\nL112-113: “The results revealed…” what results? It is not clear to me what analyses were performed.\nAre the statements in this paragraph purely based on visual inspection, or is there more to it that is\nnot mentioned? This should be clarified.\nR: Action taken. We rewrite these two paragraphs according to the reviewer’s suggestion. The first\nparagraph now primarily describes the main methodology, clearly outlining the analytical\nframework and steps. The second paragraph focuses on discussing how forest management practices\n(e.g., plantation forestry, interventions in natural forests, and artificial forests) impact forest\nresilience in specific regions. Additionally, we have supported the conclusions based on visual\ninspection with relevant figures to enhance clarity and support the arguments. Here are the\nmodifications (Lines 116~142):\n“We derived the reference distribution of forest managements for the year 2000 (Fig. 1a) by\noverlapping tree cover data available for the year 2000 derived from the Global Hansen Forest\nChange Map (GHFC) and the spatial mapping of forest management types retrieved from the Global\nForest Management Datasets (GFMD) for the year 2015. Time series of multiple vegetation indices\nwere retrieved for the selected forest pixels over the period 2001-2015, including leaf area index\n(LAI), kernel normalized difference vegetation index (KNDVI), fraction of absorbed\nphotosynthetically active radiation (FAPAR), gross primary productivity from MODIS\n(GPP_MODIS), and GPP from the Penman-Monteith-Leuning model (GPP_PML). Seasonal period\npatterns and long-term trends were removed from the original vegetation time series to make them\nstationary (Fig. 1b, d) and CSD indicators were retrieved on the residual signal to characterize the\nforest resilience. In this study, we express resilience in terms of AR1, λAR1, and λVar (additional\ndetails in the Method Section).\nResults show that areas characterized by human activity, as reflected in the forest management\nmap, generally have lower resilience than natural forest without management (NF WOM) (Fig. 1).\nThis emerges particularly evident for instance in the Amazon and Central African tropical forests\nwhere the AR1 estimated from LAI (AR1-LAI) in natural forests with management (NF WM) is\nsubstantially higher than that one quantified in NF WOM (Fig. 1j, k and Supplementary Fig. S1).\nThe higher AR1/λAR1/λVar indicates lower forest resilience. Furthermore, lower resilience is also\nobserved in Indonesian oil palm plantations and PF of southeastern United States compared to\nsurrounding NF WOM (Fig. 1f, g and Supplementary Fig. S1). In Southeast Asia, NF WM and\nagroforestry demonstrate lower resilience relative to NF WOM (Fig. 1h, i). Resilience estimated by\nusing other vegetation indices (KNDVI, FAPAR, GPP_MODIS, GPP_PML) and alternative CSD\nindicators (λAR1, λVar) exhibits similar spatial patterns (Supplementary Fig. S2) and confirm the\nrobustness of our assessment.”\nL117-118: “Forest fragments experienced higher levels of human disturbances and pressures, such\nas logging, clearing, or agriculture 38” It should be clarified that these are results from the literature,\nif so.\nL123-124: “These… influence…” check consistency of singular/plural forms; “management…”\npractices?\nR: Action taken. We have revised the sentence to explicitly indicate that these findings are derived\nfrom previous studies. Then, we have moved the related literature discussion from the result to the\ndiscussion part. The revised sentences now read as (Lines 386~394):\n“Despite this limitation, both the previous study18 and our findings illustrate the damage of human\nmanagement practices on vegetation resilience (Fig. 2). Managed forest fragments experience\nhuman disturbance and pressure, including forest planting, rotation, fragmentation, and irrigation.\nThese factors directly influence various aspects such as tree species composition, forest age, and\nphysiological processes50, disrupting forest recovery processes. Furthermore, these management\npractices also affect how vegetation responds to climate change, ultimately diminishing the\nresilience of forest ecosystems in recovering from disturbances18, 51.”\nL125-127: Again, not entirely clear if these statements are based solely on studies from the literature\nor if they simply agree with the literature.\nR: Action taken. We agree that it is essential to distinguish between conclusions drawn from our\nanalysis and statements supported by the literature. To address this, we have revised the text as\nfollows (Lines 131~134):\n“This emerges particularly evident for instance in the Amazon and Central African tropical forests\nwhere the AR1 estimated from LAI (AR1-LAI) in natural forests with management (NF WM) is\nsubstantially higher than that one quantified in NF WOM (Fig. 1j, k and Supplementary Fig. S1).”\nFor results based on our analysis, we have now explicitly referenced the supporting figures and\ntables. For statements supported by the literature, we have included proper citations and used clear\nconnectors such as \" consistent with \" to explain their context (Lines 348~349). We have moved the\nrelated literature discussion from the result to the discussion part (Lines 382~394).\nL149: it is here we discover there are paired grids. But how are they paired? some info needs to\npercolate in the main text\nL149: “paired grids” again, it should be briefly mentioned before what this method is about.\nR: Action taken. According to your comments, we have revised the text to provide a clearer\ndefinition of the paired grids method in the Introduction Section (Lines 102~104). Additionally, we\nhave specifically introduced the details of the method in Result Section. The revised texts in the\nResult Section now reads (Lines 159~165):\n“Specifically, grids are paired using the sliding window method, which moves a fixed-size window\nacross the study area. Within each window, paired grids are selected based on similar forest age\n(differences of ≤10 years), elevation (differences of ≤500 meters), and a minimum of 20% forest\ncoverage for each management type (more details in the Method Section). We point out that the\neffects of potential covariates have been ruled out within the sliding window approach.”\nL168-169: “over 40% of paired grids show higher resilience in PF than in NF” this is quite high! If\nthis is the case, I wonder if the statements about management decreasing resilience is founded or\nwhether you should be more conservative. Considering the relevant uncertainties, is the 60%-40%\nratio really different than a 50%-50% ratio based on random noise?\nR: Thank you for this important question. Our intention is not only to highlight that certain\nmanagement activities can diminish resilience, but also to emphasize the complexity and variability\nof these effects. While forest management activities associated with high human pressure—such as\nintensive plantations—can significantly reduce forest resilience, there are also management\npractices characterized by lower anthropogenic pressure that may enhance resilience under specific\nconditions. For instance, as illustrated in Figure 2, the most intensive form of human planting (i.e.,\nplantations) shows the greatest negative impact on resilience, whereas less intensive afforestation\npractices (i.e., planted forests) have a relatively weaker effect on reducing resilience.\nAction taken. Nonetheless, we recognize that the reader may be misled here, so we have rephrased\nthe paragraph topic, now reading as (Lines 184~186). “Although forest management activities with\nstrong human pressures have weakened forest resilience, managements with lower anthropogenic\npressures can enhance resilience under specific conditions (Fig. 2).”\nL174-175: “P/PET (the Ratio of Precipitation to Potential Evapotranspiration)” it is good that you\ndefine P/PET, but you should maybe also mention what it is since it is the first time you mention it\nin the text, i.e. an aridity index.\nR: Action taken. We have corrected the text in the revised version of the manuscript according to\nthe reviewer’s suggestion (Line 16). In addition, we check similar issues throughout the manuscript.\nL179: “these 15 resilience calculations…” not entirely clear what “these” refers to.\nR: Action taken. We have corrected the text in the revised version of the manuscript according to\nthe reviewer’s suggestion. The modified version of the sentence now reads as (Lines 193~198): “To\nbuild confidence in these results, analyses have been produced for an ensemble of different\ncombinations of CSD indicators and satellite vegetation indices (15 combinations = 5 vegetation\nindices × 3 CSD indicators). The overall high consistency emerging from the ensemble results\nconfirms the robustness of the differences in resilience amongst the different forest management\ntypes investigated (Fig. 2a-f and Supplementary Fig. S4).”\nL183: “Fig. S7”, what happened to figure S6?\nR: Action taken. We have ensured that all figures and supplementary figures are correctly\nreferenced in the text and are presented in the proper order.\nL221-223: “since… which…” correct the structure of the sentence. Currently the main clause is\nmissing.\nR: Action taken. We have corrected the text in the revised version of the manuscript according to\nthe reviewer’s suggestion. The modified version of the sentence now reads (Lines 231~233):\n“Notably, the number of paired grids in arid regions with P/PET less than 1 is insufficient (only 20%\nof grids are in the P/PET < 1 region), which may bring some uncertainties to the estimation of\nresilience.”\nL245: Mention what “background factors” are.\nR: Action taken. We have corrected the text in the revised version of the manuscript according to\nthe reviewer’s suggestion. We change the “background factors” to “biophysical factors” based on\nthe comment from you and Reviewer #4. We further explain it. The modified version of the sentence\nnow reads as (Lines 265~267):\n“The biophysical contexts including climate, vegetation, and soil variables may affect the resilience\ndifference between PF and NF (details and models for other paired types in Discussion).”\nIn addition, we further described all related biophysics factors in the Figure 4 caption (Lines\n304~307):\n“Abbreviations of biophysical factors include: Temp-temperature, P/PET-the ratio of precipitation\nto potential evapotranspiration, AGB-above-ground biomass, BGB-below-ground biomass, CEC-\ncation exchange capacity, RD-root depth, Sand-sand fraction, and Clay-clay fraction.”\nL249: Correct “focued” with “focused”.\nR: Action taken. It has been corrected. The modified version of the sentence now reads (Lines\n269~272): “Different vegetation indices and resilience estimation metrics show minimal impact on\nestimating the resilience difference (Supplementary Fig. S4); thus, we focus on the results for the\nλAR1 difference via LAI between PF and NF (∆λAR1-LAI ).”\n(PF, NF)\nL256: “explanatory” power?\nR: Action taken. The explanatory power of model is demonstrated by its ability to predict outcomes\naccurately like R2. We have corrected the text in the revised version of the manuscript according to\nthe reviewer’s suggestion. The modified version of the sentence now reads (Line 27):\n“with Model 1 having higher R2 than Model 2 (Supplementary Fig. S9).”\nL259: “other factors”, which ones? Maybe make a few examples.\nR: Action taken. We have corrected the text in the revised version of the manuscript according to\nthe reviewer’s suggestion. The modified version of the sentence now reads (Lines 280~282):\n“The ranking of Shapley Additive Explanations (SHAP) values suggests that climate variables, such\nas temperature and P/PET, are more important than other factors like biomass and root depth (i.e.,\nrelated to the vegetation context) (Fig. 4b).”\nL281-285: It is a bit strange that subplots e (PDSI) and i (temperature) are presented before a, b, c,\netc.\nL285: “…increases as the enhancement of…”, you mean “…increases with the increase of…”?\nR: Action taken. We have adjusted the text for clarity and logical flow. Subplot references now\nfollow a consistent order, and the “phrasing regarding the increase” has been corrected (Lines\n310~315):\n“We investigate the sensitivity of forests to water-related and energy-related factors such as soil\nmoisture (SM) and PET (Fig. 5 and Supplementary Fig. S10). Results show that forest sensitivity\nto SM decreases with the increase of P/PET, while the sensitivity of forests to PET increases with\nthe increase of P/PET (Fig. 5a-d). It indicates that forests in humid regions are limited by energy,\nwhile those in arid areas are limited by water resources.”\nL287-288: “In this regard, results further suggests that PF are more sensitive to water resources than\nNF (Fig. 5e-5j)” this is a bit hard to follow, since the results are presented via sensitivity difference.\nYou could add a sentence to explain why this is the case.\nR: Action taken. We have revised the text to explicitly explain why PF are more sensitive to water\nresources than NF, based on the sensitivity difference. The updated sentence now reads as (Lines\n315~318):\n“Furthermore, results reveal that the difference in vegetation sensitivity to water-related factors\nbetween PF and NF is over 0, suggesting that PF are more sensitive to water resources than NF in\nregions with P/PET < 1.5 (Fig. 5e-j).”\nL291: what happened to Figure S10?\nR: Action taken. We checked the full text and rearranged the figure’s references in order. The\nrevised sentences now read as (Lines 310~312):\n“We investigate the sensitivity of forests to water-related and energy-related factors such as soil\nmoisture (SM) and PET (Fig. 5 and Supplementary Fig. S10).”\nL292: sensitivity to “soil water” or “soil moisture” might be more appropriate than “water” in\ngeneral, as SM is what is shown in Fig. 5c).\nR: Action taken. We have corrected the text in the revised version of the manuscript according to\nthe reviewer’s suggestion. The modified version of the sentence now reads as (Lines 321~324):\n“Although the sensitivity difference of forests to water-related factors between PF and NF is over 0\nin regions with P/PET > 1.5 (Supplementary Fig. S11), forests sensitivity to water-related factors\nlike SM in wetter climates is below 0 (Fig. 5c).”\nL296-298: “Because the energy sensitivity difference between PF and NF gradually increases along\nwith the wetter climate, and the variations of vegetation greenness is energy-limited in humid areas”\nmissing a main clause.\nR: Action taken. Based on our overall considerations, we have removed this sentence. We have\nprovided a detailed description of the association between water sensitivity and resilience and then\nextended to energy factors. The modified version of the sentence now reads (Lines 321~327):\n“Although the sensitivity difference of forests to water-related factors between PF and NF is over 0\nin regions with P/PET > 1.5 (Supplementary Fig. S11), forests sensitivity to water-related factors\nlike SM in wetter climates is below 0 (Fig. 5c). It suggests that water resources are not the primary\nfactor influencing LAI in wetter climates. In contrast, NF shows a higher sensitivity to energy-\nrelated factors in climates where P/PET exceeds 1.5 (Fig. 5e-j), aligning with that PF may exhibit\nhigher resilience than NF in these wetter climates (Fig. 3).”\nL298-299: “As the enhancement of P/PET” does not seem grammatically correct. Change to “With\nincreasing P/PET” or “As P/PET increases”.\nAction taken. We have corrected the text in the revised version of the manuscript according to the\nreviewer’s suggestion. The modified version of the sentence now reads (Lines 327~332):\n“Forests shift from being water-limited to energy-limited as P/PET increases, and this shift may\nexplain the transition in resilience between PF and NF across varying climatic zones (Fig. 5 and\nSupplementary Fig. S10-12, details in Discussion). The conceptual diagram illustrates the reversed\nresilience patterns of PF and NF, driven by vegetation sensitivity to climate factors (Supplementary\nFig. S12).”\nL327: λPF – λNF is a newly introduced metric/abbreviation. Explain what it is. It is also not found\nin the referenced figure 5.\nL332: “…PF and NF are reducing in both…” should be “decreasing” instead of “reducing”.\nR: Action taken. We have changed the λPF – λNF to ΔλAR1-LAI . ΔλAR1-LAI\n(PF, NF) (PF, NF)\nindicates the resilience difference between PF and NF, that is estimated by λAR1 and LAI index.\nWe also update these using in the full texts (like Lines 167~176). We have modified the “reducing”\nto “decreasing”. We have addressed the following issues, now reading as (Lines 353~360):\n“Specifically, the ΔλAR1-LAI decreases in drier regions (P/PET < 1.5) with a trend of -0.003\n(PF, NF)\nyr-1, while it increases in wetter areas (P/PET > 1.5) with a trend of +0.002 yr-1 (Fig. 6c-e).\nConsidering the direction of resilience differences between PF and NF (ΔλAR1-LAI >0 in\n(PF, NF)\nP/PET < 1.5, and ΔλAR1-LAI <0 in P/PET > 1.5), these findings indicate that the absolute\n(PF, NF)\nvalue of ΔλAR1-LAI decreases in both drier and wetter zones from 2003 to 2013 (Fig. 5). This\nimplies that the difference between PF and NF reduces over time.”\nL359: “However, a notable ambiguity arises from the definition of managed forests also came to\nthis study” unclear sentence.\nR: Action taken. We have corrected the text in the revised version of the manuscript according to\nthe reviewer’s suggestion (we removed the literature discussion in the Discussion Section). The\nmodified version of the sentence now reads (Lines 382~386):\n“Human activity datasets, generated by road maps, nightlight data, and forest cover change datasets,\nhave been used to delineate intact and managed forests in the previous study. However, a notable\nambiguity arises from the definition of managed forests, as these studies and datasets did not further\ndifferentiate between various types of forest management practices49.”\nL388: “CEC” is not defined so far.\nR: Action taken. According to the reviewer’s suggestion, we have described the effects of CEC on\nresilience difference between PF and NF in the Result Section, now reading as (Lines 288~291):\n“Results of SHAP values indicate that a wetter and colder climate, higher levels of aboveground\nand belowground biomass, and a higher soil cation exchange capacity (reflecting higher soil fertility)\nmay minimize the resilience differences between PF and NF (Fig. 4b).”\nIn addition, we have added the definition of the abbreviation in the 304~307 lines of the Figure 4\ncaption.\n“Abbreviations of biophysical factors include: Temp-temperature, P/PET-the ratio of precipitation\nto potential evapotranspiration, AGB-above-ground biomass, BGB-below-ground biomass, CEC-\ncation exchange capacity, RD-root depth, Sand-sand fraction, and Clay-clay fraction.”\nL399-401: repeated sentence.\nR: Action taken. We have deleted the repeated sentence.\nL416-418: unclear sentence.\nR: Action taken. We have rewritten the text in the revised version of the manuscript according to\nthe reviewer’s suggestion. The modified version of the sentence now reads (Lines 436~441): “Trees\nin PF have minimal height variation, which enables them to receive similar levels of solar radiation\nand create a stable microclimate to counter temperature changes63. It allows PF to maintain stable\nphysiological processes under fluctuating energy conditions. NF exhibit significant height\ndifference, with shorter trees often unable to access adequate radiation due to shading by taller\ntrees64.”\nL420-422: “More interestingly, vegetation greenness is water-limited in dry regions (PF is more\nsensitive to water) and energy-limited in wet regions (NF is more sensitive to energy)”. To me, this\nsounds rather obvious. It is well known and intuitive that vegetation is limited by water where water\nis scarce, and energy becomes the limiting resource where water is abundant. Maybe you want to\nreformulate this to frame it more in the context of resilience and sensitivity.\nR: Action taken. We have corrected the text in the revised version of the manuscript according to\nthe reviewer’s suggestion. We have added an example to explain the relationship between sensitivity\nand resilience. The modified version of the sentence now reads as (Lines 427~431):\n“Existing studies have indicated the relation between resilience and climate sensitivity16, 60. For\ninstance, forests with higher sensitivity to water variability experience greater declines in biomass\nwhen subjected to external disturbances like droughts, prolonging the recovery process and making\nit more challenging for forests to recover their original state.”\nL438: I believe one needs to reach this point to learn for the first time that these indices are made\nfrom LAI\nR: Action taken. We have corrected the text in the revised version of the manuscript according to\nthe reviewer’s suggestion. To address this, we have revised the Result Section to better introduce\nthese metrics and vegetation products (Lines 116~128).\nL451: “decline” change to “decrease”.\nR: Action taken. We have corrected the text in the revised version of the manuscript according to\nthe reviewer’s suggestion. Since decline is an intransitive verb, we have checked the whole\nmanuscript and replaced “decline” with “decrease” or “reduce”.\nL458: “the primary focus is on the differential impacts of forest management between PF on\nresilience” unclear sentence.\nR: Action taken. We have rewritten the text in the revised version of the manuscript according to\nthe reviewer’s suggestion, now reading as (Lines 514~520): “This study prominently examines the\ndifferences in resilience between PF and NF, further focusing on their response to biophysical\nfactors. Certain factors (e.g., landforms and slope), that influence the vegetation recovery process,\nmay not have been adequately accounted for in the model70. Furthermore, we also built XGBT\nmodels to investigate the influence of biophysical contexts on the resilience differences between NF\nand Plantations, NF WM and NF WOM, NF WOM and F WM, and between NF WOM and\nPlantations (Supplementary Fig. S29).”\nL480: “…and reduce uncertainties from land use change…” and following sentence: to me, this\nsounds more like land use change was ignored since you are only assuming that pixels do not change,\nbut in reality, they might be changing. Am I missing something?\nR: Thank you for your important question. Vegetation changes can influence the comparison of\nresilience under different forest management practices. Therefore, it is necessary to minimize this\ninfluence. We achieved this by overlaying the GFMD and GFHC maps. If a pixel was classified as\na forest in both 2000 and 2015, we assume/consider that no land-use change occurred in the pixel.\nThen, we only analyze pixels with no landcover change. Action taken. We have corrected the text,\nnow reading as (Lines 558~560): “This process ensured the spatial consistency of the pixels during\n2000-2015, and pixels with land cover change may not be accounted for in the analysis.”\nL491: “ASCE” define/explain what it means.\nR: Action taken. We have corrected the text in the revised version of the manuscript according to\nthe reviewer’s suggestion, now reading as (Lines 576~580): “we used monthly mean precipitation,\nmaximum and minimum air temperature at 2 m, soil moisture, PDSI (Palmer Drought Severity\nIndex), and potential evapotranspiration (PET that based on the American Society of Civil Engineers\nPenman-Monteith equation75), acquired from the TerraClimate database at 4-km (1/24th degree) for\n2001-2015 period76. ”\nL503: change “photosynthetic” to “photosynthesis”.\nR: Action taken. We have corrected the text in the revised version of the manuscript according to\nthe reviewer’s suggestion, now reading as (Lines 590~591): “Different vegetation indices were used\nto capture the forest biomass, greenness, and photosynthesis.”\nL506: “PML” and “MODIS” you may want to explain the abbreviations or mention the source\n(which for MODIS you do a few lines below).\nR: Action taken. We have corrected the text in the revised version of the manuscript according to\nthe reviewer’s suggestion, now reading as (Lines 591~596): “Here we selected five vegetation\nindices including LAI (leaf area index), FPAR (fraction of photosynthetically active radiation),\nGPP_PML (gross primary productivity from Penman-Monteith-Leuning model), GPP_MODIS\n(gross primary productivity from Moderate Resolution Imaging Spectroradiometer), and KNDVI\n(kernel-based normalized difference vegetation index) to monitor the dynamic changes in vegetation\nfrom 2001-2015. ”\nL511: this is a wrong use of kNDVI. see https://doi.org/10.1016/j.isprsjprs.2022.12.019\nR: There are various methods for calculating KNDVI. The formula for calculating the KNDVI is:\n(KNDVI = tanh ((NIR - Red)/2σ)^2) 13. Where NIR is the near-infrared reflectance, Red is the red\nband reflectance, and σ is a scaling parameter that normalizes the input. Previous studies have\ndemonstrated that σ = 0.5(NIR + red) is a reasonable choice for related calculations14. Then, the\nequation can be simplified as KNDVI = tanh(NDVI²). This formulation ensures that the scaling\nparameter σ adapts dynamically to the reflectance values of NIR and Red, thereby preserving the\nphysical interpretation of the vegetation index while minimizing sensitivity to noise in the\nreflectance data. The inclusion of σ accounts for variability in reflectance magnitudes, making the\nindex robust to differences in illumination and sensor calibration. In addition, many previous studies\nused the simplified equation of KNDVI4,15. Action taken. We have corrected the text in the revised\nversion, now reading as (Lines 597~598): “It is calculated using the simplified formula KNDVI =\ntanh(NDVI²),”\nL516: “which are mainly from [MODIS]” and where are the other products from?\nR: Action taken. We rewrite the description of the vegetation products according to the reviewer’s\nsuggestion, now reading as (Lines 599~608):\n“These vegetation indices were derived from various vegetation products. Specifically, GPP_PML\nis derived from the PML_V2.0.1.7 model82, and GPP_MODIS is derived from MODIS Terra\n(MOD17A2H.061) and Aqua (MYD17A2H.061) observations. These datasets offer a spatial\nresolution of 500 meters and an 8-day temporal resolution. LAI and FPAR are obtained from\nMYD15A2H.061 and MOD15A2H.061 MODIS Aqua and Terra observations, also at a 500 m\nspatial resolution and 8-day temporal resolution. NDVI is derived from MOD13A2.061 and\nMYD13A2.061 MODIS Terra and Aqua observations, with a 1 km spatial resolution and 16-day\ntemporal resolution. The details of these products can be found in Supplementary Materials (Table\nS2).”\nL518-519: “these two observed vegetation greenness or biomass are similar in a day” give a\nreference stating this.\nL519: “But two times may be affected by the cloud” unclear sentence.\nR: Action taken. We have corrected the text in the revised version of the manuscript according to\nthe reviewer’s suggestion, now reading as (Lines 608~611):\n“Furthermore, the Terra and Aqua satellites observe the earth's surface at 10:30 am and 1:30 pm,\nand their observations of vegetation greenness or biomass are similar on a given day83. These\ndatasets may be affected by clouds, snow, and shadowing.”\nL524-526: You mention two contrasting spatial resolutions for the vegetation indices (500 m and 1\nkm). Which one is it?\nR: Action taken. We have corrected the text in the revised version of the manuscript according to\nthe reviewer’s suggestion. We have revised the text to ensure clarity and explicitly state the temporal\nand spatial resolutions for all datasets (Lines 599~608).\nL530: you still have not explained how to calculate resilience. why not start with that before?\nR: Action taken. We have corrected the text in the revised version of the manuscript according to\nthe reviewer’s suggestion. We have placed the Estimating Vegetation Resilience subsection after the\ndataset’s subsection, followed by matching methods (Lines 617~706).\nL531: If not already mentioned in the introduction, the “paired grids” concept should really be\nexplained in detail here, as I am still wondering what it means.\nAction taken. We have corrected the text in the revised version of the manuscript according to the\nreviewer’s suggestion. We agree that this is a critical methodological element and have ensured it is\nexplained in detail across the manuscript, specifically in the Introduction, Result, and Method\nSections.\nIn the Introduction:\nWe introduced the paired grids concept as part of the methodology in the Introduction Section (Lines\n102~104):\n“To differentiate the impacts of human activities from natural variability, we employ a paired\nmethod that compares geographically proximate forest areas under different management practices.”\nIn the Result Section:\nWe elaborated on how paired grids were constructed, specifying the use of the sliding window\napproach to ensure similarity in environmental conditions (Lines 159~165):\n“Specifically, grids are paired using the sliding window method, which moves a fixed-size window\nacross the study area. Within each window, paired grids are selected based on similar forest age\n(differences of ≤10 years), elevation (differences of ≤500 meters), and a minimum of 20% forest\ncoverage for each management type (more details in the Method Section). We point out that the\neffects of potential covariates have been ruled out within the sliding window approach.”\nIn the Method Section:\nWe have added the definition and significance of paired grids. In the Method Section, we have\nintroduced the details of the paired grid method (Lines 679~706).\nL538: what is Simon et.al?\nR: Action taken. We have corrected the text in the revised version of the manuscript according to\nthe reviewer’s suggestion. The modified texts now read as (Lines 690~691): “The global forest age\ndata are obtained from Besnard et al.88, while elevation data were derived from the Shuttle Radar\nTopography Mission (SRTM) v4 dataset89.”\nL541-543: mention supplementary figures explicitly.\nR: Action taken. We have corrected the text in the revised version of the manuscript according to\nthe reviewer’s suggestion. In response, we have reviewed the manuscript thoroughly to ensure all\nreferences to supplementary figures are explicitly stated and correctly cited.\nL546: why monthly data? to do a robust time series analysis it is better to take all data available, not\npicking a single MVC value per month\nR: We agree that robust time series analysis benefits from utilizing as much temporal data as possible,\nand we appreciate your suggestion. To address this concern, we conducted a sensitivity analysis to\nexplore the impact of different temporal preprocessing methods on resilience estimation (Lines\n765~781 for the methods). We tested the impact of the three preprocessing methods on the results,\nincluding (1) obtaining monthly maximum vegetation data via MVC from MODIS Terra and Aqua\nsatellites, (2) acquiring original temporal resolution vegetation data from Terra satellite, (3) and\nobtaining original temporal resolution vegetation data from Terra satellite only for the growing\nseason. Our results indicate that the choice of temporal preprocessing method does not significantly\naffect the main conclusions of our study. The related sensitivity test can be found in Supplementary\nFig. S20-S23, and the related discussion can be found in Lines 457~462 and Supplementary Note 2.\nAction taken. We have corrected the text in the revised version to mention the sensitivity test (Lines\n674~677).\n“Some CSD studies even suggest that missing values should be excluded in resilience estimating,\nas this could induce new uncertainty by interpolation methods45. In this regard, we conducted\nadditional tests to evaluate the impact of data pre-processing methods on our key findings\n(Supplementary Fig. S20-S23).”\nL550: you need to explain what these are\nR: Action taken. We have corrected the text in the revised version of the manuscript according to\nthe reviewer’s suggestion, now reading as (Lines 636~639): “Consequently, we evaluated forest\nresilience using AR1, λAR1, and λVar. Notably, these λ metrics (λAR1 and λVar) capture resilience\nwithin a linearized framework of the system’s dynamics near a stable equilibrium, which can be\nformally modeled as an Ornstein–Uhlenbeck process.”\nL552: It is still unclear to me what the difference between AR1 and λAR1 is.\nR: In fact, AR1 and λAR1 are closely related. Mathematically, λAR1 is simply the logarithm of AR1.\nIn other words, λAR1 excludes cases where AR1 is less than 0. In ecological systems, recovery rate\nindicators typically assume that the system exhibits positive autocorrelation, meaning that past data\ninfluence current or future conditions (for example, if a drought leads to a decline in greenness in\nan earlier period, the drought’s effect on greenness can still be observed in later data). However, in\nthe dynamic real world, conditions can be complex and variable. Therefore, we retain both AR1 and\nλAR1 for analysis (because we also need to test this assumption). Here are the original texts to\ndescribe the relation between AR1 and λAR1 (Line 665): “λAR1 and λVar may result in null values\n∅2\ndue to logarithmic constraints (i.e., AR1 < 0 or >1).”\n𝑉𝑎𝑟[𝑋]\nFormulas: try to keep the font of parameters consistent between the formulas and the main text.\nR: Action taken. We have addressed this issue by reformatting all formulas to ensure consistency.\nThe font for parameters in the formulas has been standardized to Cambria Math, with a font size\nequivalent to \"12 pt\". We have reviewed the entire manuscript to ensure that the parameters in the\nformulas match those in the main text, both in terms of notation and style (see equations).\nL571: “same as the Xn” unclear.\nR: Action taken. We checked the full text for italicized notations about formulas and variables.\nWe have corrected the text in the revised version of the manuscript according to the reviewer’s\nsuggestion, now reading as (Lines 650~651): “Where X denotes the system state at discrete time\nn+1\nstep n+1 (where X is the system state at discrete time step n)”\nn\nL591-593: correct typing mistakes: “however”, “maxium”, “northern”.\nR: Action taken. We have corrected the text in the revised version of the manuscript according to\nthe reviewer’s suggestion.\nL591: Explain what “MVC” stands for.\nL592: Explain what “month datasets” are, or reformulate.\nL593: “influencing” reformulate, unclear.\nR: Action taken. We have corrected the text in the revised version of the manuscript according to\nthe reviewer’s suggestion. We have addressed the following points: MVC stands for Maximum\nValue Composite (we first introduce the MVC in Line 614), a method commonly used to handle\nmissing data in satellite-derived vegetation indices. This has been clarified in the revised text. The\nphrase \"month datasets\" has been reformulated for clarity as \"monthly vegetation datasets.\"\nThe\nterm \"influencing\" has been replaced with \"affect,\" and the sentence structure has been adjusted to\nenhance readability and clarity, now reading as (Lines 671~672):\n“Furthermore, to mitigate the impact of missing values on resilience estimation, we applied the\nMVC method to generate monthly vegetation data.”\nL595: “some CSD studies did not deal with missing values”, does that mean that missing values are\nexcluded or kept as they are? Either way, it is not possible to “not deal with missing values”.\nR: Action taken. We agree that the phrase \"not deal with missing values\" is unclear and could lead\nto misunderstanding. We have corrected the text in the revised version of the manuscript according\nto the reviewer’s suggestion. The revised sentence now reads (Lines 674~677):\n“Some CSD studies even suggest that missing values should be excluded in resilience estimating,\nas this could induce new uncertainty by interpolation methods45. In this regard, we conducted\nadditional tests to evaluate the impact of data pre-processing methods on our key findings\n(Supplementary Fig. S20-S23).”\nL603-605: “Not only because of the strong relationship between resilience difference between PF\nand NF and background variables, …” missing main clause. Also, this part specifically sounds as if\nyou chose the method based on the results because “it looks good”, instead of basing the decision\non hypotheses or research gaps.\nR: Action taken. We have corrected the text in the revised version of the manuscript according to\nthe reviewer’s suggestion. Here are the modified texts, now reading as (Lines 708~713):\n“We first explored the biophysical drivers for resilience differences between PF and NF. This is\nparticularly important given the widespread attention to restoration activities21, 22, 23, such as the UN\nDecade on Ecosystem Restoration (2021–2030). Understanding the resilience differences between\nPF and NF can provide valuable guidance for designing and implementing effective restoration\nstrategies to enhance ecosystem resilience.”\nAdditionally, following the suggestions of Reviewer #4, we have included comparisons involving\nother forest management types. Therefore, the main text focuses primarily on the differences\nbetween PF and NF, while comparisons of other types are included in the Discussion Section (Lines\n514~533) and Supplementary Materials (Supplementary Note 4) with relevant descriptions.\nL607: “widely utilized in ecosystem studies”, then maybe give a few more references instead of\nonly one.\nR: Action taken. We have corrected the text in the revised version of the manuscript according to\nthe reviewer’s suggestion. Here are the modified texts, now reading as (Lines 715~717):\n“XGBT combines gradient learning and tree-boosting steps, making it a powerful ensemble learning\nalgorithm widely utilized in ecosystem studies39, 90, 91. ”\nL612: what does “its” refer to?\nL615” what does “their” refer to?\nR: Action taken. We have corrected the text in the revised version of the manuscript according to\nthe reviewer’s suggestion. Here are the modified texts, now reading as (Lines 728~733):\n“We used XGBT models to predict the λAR1 difference between PF and NF (e.g., △λAR1-LAI\n= λAR1-LAI – λAR1-LAI ), based on eight environmental factors and the differences\n(PF,NF) (PF) (NF)\nin these factors between PF and NF. The AR1 difference between PF and NF was the target variable,\nwhile the eight environmental factors (e.g., CEC value of the grid) and their differences (e.g., △\nCEC = CEC – CEC ) between PF and NF were used as predictors.”\n(PF,NF) (PF) (NF)\nL621: “datasets availability section” should be corrected to “Data and code availability section”.\nR: Action taken. Based on the feedback from yours and Reviewer #3, we have added descriptions\nof the relevant variables in the Method Section. The revised sentence now reads (Lines 735~740):\n“The eight environmental factors include annual temperature (Temp derived from TeerraClimate),\ndryness index (P/PET derived from TeerraClimate), root depth (RD92), above-ground biomass\n(AGB93), below-ground biomass (BGB93), clay fraction (Clay94), sand fraction (Sand94), and cation\nexchange capacity (CEC94). These factors are selected for their known influence on forest resilience,\nincluding both climatic and land-surface variables. ”\nL671: “low probability pixels” what does this refer to, low probability of what?\nR: Action taken. To clarify, this refers to pixels with a low probability of accurate classification\nwithin the Global Forest Management Dataset (GFMD). Specifically, these are pixels where the\nclassification accuracy of forest management types (e.g., planted forests, natural forests) is below a\ncertain threshold (50%, 60%, 70%) based on the GFMD’s internal validation metrics. We have\ncorrected the text in the revised version of the manuscript according to the reviewer’s suggestion,\nnow reading as (Lines 788~791): “By removing low-accuracy pixels (i.e., low accuracy of\nclassification in GFMD), the results of reanalysis displayed a similar pattern, and different accuracy\nthresholds may have a limited impact on our key conclusions.”\nL681-683: “GFMD is the result of a classification algorithm that takes in input various factors\nincluding vegetation data, the same data to calculate resilience. It may undermine the independence\nbetween forest management and resilience”. This sounds like a potentially huge issue. If the input\ndata are not independent, how can you say anything about the results? The following parts were you\ncontrol with independent data should really be mentioned more centrally.\nL702-703: “This is aligned with our main conclusion. The above analysis results can be seen in the\ndiscussion and supplementary section” good, as mentioned above, I think this confirmation should\nbe more central to the results. Also, I could not find it in the discussion, where it should really be,\nrather than hidden at the end of the methods or in the supplementary materials.\nR: Action taken. Thank you for your valuable suggestions. We agree that the discussion on the\nindependence and reliability of the GFMD data should be more critical to the Discussion Sections.\nThe specific modifications, including details of the expanded discussion, are outlined in the above-\nmentioned response (Modified texts can be read in Lines 496~491 of the Discussion Section).\nL707: the data here means the data produced by the paper, not just the data used as input (although\nthe latter is also needed)\nR: Action taken. We updated the core datasets for the figure’s reproduction. We have corrected the\ntext in the revised version of the manuscript according to the reviewer’s suggestion. Here are the\nmodified texts: (Lines 835~838)\n“The core part of datasets related to the figures and maps are now can be found:\nhttps://figshare.com/articles/dataset/_b_Biophysical_factors_determining_the_effects_of_forest_\nmanagement_on_ecosystem_resilience_b_/28286651”\nResponse to Reviewer #2\nReviewer #2 (Remarks to the Author):\nI co-reviewed this manuscript with one of the reviewers who provided the listed reports. This is part\nof the Nature Communications initiative to facilitate training in peer review and to provide\nappropriate recognition for Early Career Researchers who co-review manuscripts.\nR: We greatly appreciate the time and effort you have dedicated to reviewing our manuscript and\nproviding important feedback. We have modified and responded point-by-point accordingly.\nResponse to Reviewer #3\nReviewer #3 (Remarks to the Author):\n1. This study explores the resilience of forests under varying management intensities, ranging from\nunmanaged natural forests to managed, planted forests and plantations. The research examines how\nresilience, as defined by the critical slowing down (CSD) theory, differs among forests with different\nlevels of management. The study concludes that natural, unmanaged forests exhibit greater\nresilience, while planted forests are more vulnerable and thus show lower resilience. Additionally,\nthe study investigates the factors influencing these differences, reporting that wetter regions are less\nvulnerable and, consequently, more suitable for planted forests.\nI find the study both interesting and novel, contributing important findings beyond previous similar\nstudies. However, several points need to be addressed to strengthen the results. The authors\nconducted multiple sensitivity analyses to demonstrate that the findings are consistent, considering\nfactors such as the sliding window width, forest coverage threshold, definitions of planted (PF) and\nnatural forests (NF), reliability of forest management maps, and preprocessing methods of satellite\nvegetation data. However, one crucial aspect that should be considered is forest age: natural and\nmanaged forests are typically older than planted forests and plantations, and resilience is known to\nvary among different age groups. Therefore, I recommend including an additional sensitivity\nanalysis to compare planted forests with younger unmanaged stands, to verify that the results are\nnot merely due to age differences between PF and NF.\nR: We sincerely appreciate your detailed comments and constructive feedback. Below, we address\nyour key suggestions on the forest age. We acknowledge that forest age is an important factor\ninfluencing resilience, and the differences between younger and older forests could play a role in\nthe resilience comparison between planted forests (PF) and natural forests (NF). In this study, we\nhave already accounted for the potential impact of forest age on our findings (Lines 687~691). In\nthe Method Section, we used the paired grids method to compare vegetation resilience differences\namong different forest management types. We have excluded the paired grid with more than 10\nyears of forest age difference. For example, if the age difference between a paired PF and NF\nexceeded this threshold, the grid was excluded from the analysis. Here are the original texts in the\nmain text: “To mitigate the potential effects of varying forest age and elevation on vegetation\nresilience, we ensured that the differences in these variables within each paired grid were restricted\nto within 10 years and 500 meters, respectively16. The global forest age data are obtained from\nBesnard et al.17, while elevation data are derived from the Shuttle Radar Topography Mission\n(SRTM) v4 dataset18.”\nAction taken. We noticed that the critical information about paired grids and forest age difference\nwas mentioned only once in the Method Section, and the manuscript did not adequately introduce\nor emphasize this important methodological approach. To ensure readers do not overlook these key\ndetails, we have mentioned the description of the paired grids method and the age thresholds in the\nResult Section. These revisions include as following: (Lines 159~165)\n“Specifically, grids are paired using the sliding window method, which moves a fixed-size window\nacross the study area. Within each window, paired grids are selected based on similar forest age\n(differences of ≤10 years), elevation (differences of ≤500 meters), and a minimum of 20% forest\ncoverage for each management type (more details in the Method Section). We point out that the\neffects of potential covariates have been ruled out within the sliding window approach.”\nAdditionally, according to the reviewer’s comment, we conducted a sensitivity analysis by applying\ndifferent age difference thresholds (e.g., 5, 7, and 10 years) to filter paired grids. The results\nconsistently showed similar findings: forest management tends to weaken resilience, but under\ncertain conditions, well-managed PF exhibited resilience levels comparable to or even exceeding\nthose of NF (Supplementary Fig. S19). Importantly, the findings remained robust regardless of the\nage threshold applied. Here are modified texts and figures: (Lines 753~757)\n“Specifically, we tested the impact of the sliding window widths (0.25°, 0.35°, and 0.45°),\nthe\nthreshold of forest age differences (5, 7, and10 years), and forest coverage thresholds (10%, 15%,\nand 20%) on the results, and found that they did not affect the key conclusions (Supplementary Figs.\nS17-S19).”\nFigure S19. The impacts of forests' age difference on our main findings. Similar to Figure 2 (10\nyears), but with different thresholds of forests age difference (5 and 7 years).\n2.Furthermore, in the final chapter discussing potential limitations, the authors should expand on\nthe limitations of the CSD theory.\nR: We agree with the importance of elaborating on the limitations of the critical slowing down (CSD)\ntheory, as you suggested. Action taken. We have expanded on this point in the revised Discussion\nSection, highlighting the challenges of applying CSD in real-world ecosystems, such as (a) the\ninfluence of environmental noise, (b) the assumptions of stationarity, (c) the difficulty in detecting\nthe time of resilience shifts, (d) and the impacts of the optical index on the accurate ecosystem\nresilience estimating. Here are the modified texts as follows: (Lines 492~513)\n“It is important to acknowledge the potential limitations of CSD. Based on the CSD theory, changes\nin CSD metrics like λAR1 serve only as early warning indicators for resilience recovery and loss15,\n67. These indicators exert challenges to accurately pinpoint the time of critical transitions in system\nstates. Furthermore, CSD metrics like variance are easily influenced by environmental noise, which\ncan lead to inaccurate detection of slowing-down events43. For example, previous studies reported\nthe critical slowing down in the Amazon Forest20, 32, while the variance metric shows a decrease44.\nCSD theory assumes that systems exhibit stationary behavior. However, in ecosystems, external\ndrivers such as abrupt events and human disturbances often lead to non-stationary conditions, which\ncan weaken the accuracy of early warning signals derived from CSD15. Moreover, CSD indicators\nmay be scale-dependent because ecological processes function across multiple spatial and temporal\nscales45. Ignoring this dependence may result in the uncertainties of resilience estimation.\nAdditionally, optical remote sensing indices may affect the accuracy of resilience estimates. Since\ngreenness indices such as the LAI are susceptible to saturation, assessing resilience using these\nindices in areas with high forest coverage may induce uncertainty45, 68, 69. It is thus recommended to\nuse biomass indices such as GPP or VOD instead. Due to the coarse resolution of VOD indices\n(approximately 0.25°), we use the high-resolution GPP index to validate our findings, which show\nsimilar patterns to the results by optical remote sensing indices (more details in Supplementary\nmaterials). These limitations highlight the importance of careful application and interpretation of\nCSD-based resilience assessments.”\n3.Aside from these points, I have only a few minor comments detailed below:\nAbstract: The abstract lacks specific mention of the aspect of resilience studied. I suggest\nmentioning the CSD theory in the abstract.\nIntroduction: The introduction is well-prepared, effectively outlining the concept of resilience,\nidentifying the research gap, and clearly stating the study's objectives.\nMethods: The authors utilized validated datasets, such as GFMD, the TerraClimate database, and\nMODIS. Although the methods involve some subjective choices, the sensitivity analysis provides\nreassurance about the robustness of the results. The classification of resilience into stable, increasing,\nor decreasing is based on a well-established theory and is therefore justified.\nR: We appreciate your positive comments on the Introduction and Method Section. Action taken.\nIn response to the comment on the Abstract Section, we have made the following revisions to\naddress your concerns: (a) Mention of CSD Theory, (b) Clarification of Resilience Quantification,\nand (c) Clarification of method and datasets in the abstract. (Lines 4~22)\n“Policymakers are concerned that newly planted trees may not provide ecosystem services as\neffectively as native forests. The efficient functioning of ecosystem services relies on the capacity\nof forest ecosystems to recover from external disturbances (i.e., ecosystem resilience). This study\nutilizes the critical slowing down (CSD) theory to quantify resilience differences among forest\nmanagement types, particularly between planted forests (PF) and natural forests (NF). CSD\nresilience indicators such as AR1, λAR1, and λVar are estimated from satellite-based vegetation\nindices. This study further investigates the biophysical drivers for the resilience difference between\nPF and NF by machine learning. Results show that most forest management activities weaken forest\nresilience. Forest resilience experiences greater losses under intensive managements However, some\nforest managements can enhance resilience under specific conditions, such as wet climates. PF\ndisplay higher resilience than NF in P/PET (the ratio of precipitation to potential\nevapotranspiration) > 1.5 (wetter climates), while NF is more resilient in P/PET < 1.5 (drier\nclimates). We further detect such resilience differences in all climate zones between PF and NF\ndecreasing from 2001 to 2015. These differences are influenced by biophysical contexts such as\nlocal climate, vegetation, and soil. Regions characterized by cold, wet conditions, dense vegetation,\nand higher soil fertility minimize differences in resilience, suggesting that rational site selection for\nafforestation can enhance forest resilience compared to native forests.”\n4.Detailed Comments:\nLine 106: The phrase \"In Spatial\" is unclear—please clarify.\nR: Action taken. To clarify, we intended to convey that the analysis results are based on spatial\npatterns. We have revised the text to improve clarity and readability. The modified text now reads\n(Lines 129~131): “Results show that areas characterized by human activity, as reflected in the forest\nmanagement map, generally have lower resilience than natural forest without management (NF\nWOM) (Fig. 1).”\nLine 108: Introduce the abbreviations GHFC and GFMD here (or in the Introduction) as they are\nused for the first time.\nR: Action taken. We have corrected the text in the revised version of the manuscript according to\nthe reviewer’s suggestion. The revised text now reads as (Lines 116~119):\n“We derived the reference distribution of forest managements for the year 2000 (Fig. 1a) by\noverlapping tree cover data available for the year 2000 derived from the Global Hansen Forest\nChange Map (GHFC) and the spatial mapping of forest management types retrieved from the Global\nForest Management Datasets (GFMD) for the year 2015.”\nIn addition, we carefully reviewed the manuscript for similar instances where abbreviations were\nintroduced without prior definition and made the necessary changes to ensure consistency and\nclarity (like AR1-LAI).\nLine 176: Could you provide an explanation for this pattern already here? Why are the results\ndifferent in these regions?\nR: Thank you for your important question regarding the observed differences in resilience between\nPF and NF under varying P/PET (precipitation to potential evapotranspiration) conditions. Our\nfindings suggest that these differences are closely tied to the sensitivity of vegetation to water and\nenergy availability, which vary across dry and wet regions. In dry regions, where forest greenness\nis primarily limited by water availability, PF exhibit higher sensitivity to changes in water resources\n(Fig. 5). This means that PF are more vulnerable to disturbances such as droughts, which reduce\nsoil moisture and water availability, leading to more severe degradation and slower recovery.\nIn contrast, in wet regions, where forest greenness is mainly limited by energy availability, NF show\ngreater responsiveness to variations in energy (Fig. 5). This difference in sensitivity to water and\nenergy explains the observed resilience patterns between PF and NF across wet and dry regions\n(Figs. 3–5).\nWe hypothesize that the observed resilience differences may be related to the sensitivity of\nvegetation to water and energy availability. For example, if forests are more sensitive to water\nresources, disturbances such as droughts will result in more severe degradation due to the reduction\nin soil moisture and water availability. Consequently, the recovery process will be prolonged as the\nforest requires more time to accumulate resources and regenerate, especially when the regeneration\nrate is limited by drought-induced stress. Over time, this will lead to a decrease in the overall\nrecovery rate of the forest ecosystem.\nAdditionally, we already discussed the reasons why are resilience different in these regions (Lines\n427~448). Other factors may also contribute to this pattern: The species planted in PF are often\nselected for fast growth and may lack the physiological adaptations in some dry regions27.\nAdditionally, in humid regions, management interventions in PF (e.g., thinning) might mitigate\nenergy limitations more effectively than NF, enhancing PF resilience28.\nAction taken. We have added the sentence for clarity and to guide the reader to a detailed\ninterpretation of the Discussion Section. The updated sentence now reads as (Lines 228~231): “This\nthreshold likely reflects the point at which water availability becomes sufficient to alleviate stress\nfrom human managements, enabling PF to exhibit superior resilience to NF (more details in the\nDiscussion Section).”\nLine 183: Please check for grammatical accuracy—consider changing to \"showing similar patterns\"\nwithout \"the.\"\nR: Action taken. We have revised the sentence for improved accuracy and clarity. The updated\nsentence now reads as (Lines 139~142):\n“Resilience estimated by using other vegetation indices (KNDVI, FAPAR, GPP_MODIS,\nGPP_PML) and alternative CSD indicators (λAR1, λVar) exhibits similar spatial patterns\n(Supplementary Fig. S2) and confirm the robustness of our assessment.”\nLine 242: Instead of \"cycle,\" did you mean \"circle\"?\nR: Action taken. We have corrected the text in the revised version of the manuscript according to\nthe reviewer’s suggestion. The revised text now reads as (Lines 261~263):\n“In the figure, a cross (×) indicates the mean value for each vegetation index, a circle (○) denotes\nthe median value, and a red inverted triangle (▼) signifies the aggregated average of the means and\nmedians.”\nLine 320: Please provide references to the \"previous studies\" that also documented a global decline\nin resilience.\nR: Action taken. We have corrected the text in the revised version of the manuscript according to\nthe reviewer’s suggestion. The revised text now reads as (Lines 348~349):\n“This analysis reveals a general decline in the resilience of global forests, which is consistent with\nprevious studies18, 48. ”\nLine 425: I suggest adding a sensitivity analysis comparing younger NF and younger PN here.\nAdditionally, consider mentioning the limitations of the CSD theory.\nR: Action taken. We have addressed the related points in responses (above-mentioned), including\na discussion of the limitations of the CSD theory (Lines 492~513).\nLine 538: The reference \"Simon et al.\" is missing from the reference list.\nR: Action taken. We got the rules of referencing wrong and should have referenced the last name.\nWe have corrected the text in the revised version of the manuscript. The revised text now reads as\n(Lines 690~691):\n“The global forest age data are obtained from Besnard et al.88, while elevation data were derived\nfrom the Shuttle Radar Topography Mission (SRTM) v4 dataset89.”\nLines 614-615: Why didn't you use age as a predictor variable? Age could significantly impact\nresilience and potentially affect your findings. Please elaborate on this. Since PF are generally\nyounger than NF, I recommend an additional sensitivity analysis comparing PF and NF of similar\nages.\nR: Thank you for your important suggestion regarding the potential impact of forest age on\nresilience and its role as a predictor variable. Below, we elaborate on our reasoning and the steps\nwe have taken to address this concern:\nControl for Forest Age in Pre-Processing:\nIn our analysis of resilience differences between planted forests (PF) and natural forests (NF), we\ncontrolled for large disparities in forest age by excluding paired grids where the age difference\nexceeded 10 years. This step minimizes the confounding effect of forest age on resilience differences.\nSimilarly, we also controlled for elevation by excluding grids with an elevation difference exceeding\n500 meters. These pre-processing steps ensure that the observed resilience differences are primarily\nattributable to forest management practices and other environmental factors, rather than large\nvariations in forest age or elevation.\nAdditional Sensitivity Analysis:\nIn response to your suggestion, we conducted additional sensitivity analyses by applying varying\nage difference thresholds (e.g., 5, 7, and 10 years) to the paired grids (Supplementary Fig. S19). The\nresults consistently showed that the observed patterns of resilience differences between PF and NF\nremained robust across these thresholds. This reinforces the reliability of our findings and indicates\nthat forest age differences, within the tested thresholds, do not significantly alter the overall\nconclusions.\nRationale for Not Using Age as a Predictor Variable:\nWhile forest age is undoubtedly an important factor, including it as a predictor variable in\nsubsequent modeling would not be appropriate. Forest age was already accounted for in the pre-\nprocessing stage through paired grid filtering. Using it again as a predictor variable in the model\nwould risk introducing redundancy and circular reasoning (i.e., independence), potentially biasing\nthe results. Therefore, to maintain the independence of the model and avoid overfitting, we did not\ninclude forest age, canopy coverage, or elevation as predictor variables in our analysis.\nLines 617-621: Although you mention that \"soil and vegetation variables\" are detailed in the dataset\navailability section, I believe these should be briefly described in the Methods section, as they are\ncritical inputs for your model.\nR: Action taken. We agree that soil and vegetation variables are critical inputs for our model and\nshould be briefly described in the Method Section. In response, we have revised the text to include\nthe following details (Lines 735~740):\n“The eight environmental factors include annual temperature (Temp derived from TeerraClimate),\ndryness index (P/PET derived from TeerraClimate), root depth (RD92), above-ground biomass\n(AGB93), below-ground biomass (BGB93), clay fraction (Clay94), sand fraction (Sand94), and cation\nexchange capacity (CEC94). These factors are selected for their known influence on forest resilience,\nincluding both climatic and land-surface variables.”\nLine 650: I couldn't find the reference to \"Taylor et al.\"—please ensure it's included in the reference\nlist.\nR: Action taken. We have corrected the error by replacing the incorrect citation with the appropriate\nreference. The revised sentence now reads (Lines 767~769):\n“Smith et al. calculated resilience by using all available records from a single satellite45, while\nForzieri et al. used growing season data from a single satellite18.”\nAdditionally, we have carefully reviewed the entire manuscript to ensure that all in-text citations\ncorrespond accurately to the references listed in the bibliography.\nI have not reviewed the code, but I have confirmed that it is available online as stated by the authors.\nR: Thank you for confirming the availability of the code online. As part of our commitment to\ntransparency and reproducibility, we will further organize and document the code after the\nmanuscript is accepted. The finalized version of the code will be permanently archived on Fig share,\nensuring long-term accessibility for the research community.\nResponse to Reviewer #4\nReviewer #4 (Remarks to the Author):\nReview report\n1. The study is focused on the forest management practices (human factors) and the (natural factors)\ndriving forest ecosystem resilience. There are sustainable forest management practices that support\nforest ecosystem resilience while other unsustainable practices undermine forest ecosystem\nresilience. This is a very interesting topic, especially as the resilience of forest ecosystems is critical\nto the supply if ecosystem services benefiting people. The effect of natural factors that also\ninfluence forest ecosystem resilience, and that is a key objective of the manuscript is not capture in the\nmanuscript title? The resilience of forest ecosystems is not only explained by forest management\npractices (type and intensity) but also more so, by biophysical factors, including climate change,\nand as such should be captured in the title.\nBiophysical factors affected the impacts of forest management on ecosystem resilience?\nR: Thank you for your important suggestion and constructive feedback. We agree that forest\necosystem resilience and the interplay between human activities and climate variables (i.e.,\nbiophysics factors). This study has indeed examined the influence of biophysical factors on the\nresilience differences between planted forests (PF) and natural forests (NF), which reflects the\nimpact of forest management on resilience. For example, we found that regions with higher P/PET\n(precipitation to potential evapotranspiration ratio) show reduced resilience differences between PF\nand NF. It indicates that afforestation in wet zones may minimal the resilience gap between planted\nforests and local natural forests.\nAction taken. Below, we address your concerns regarding the manuscript title. To address this, we\npropose revising the manuscript title to better capture the dual focus on forest management practices\nand biophysical factors. By combing your advances and the discussion with the author group, the\nmodified new title now reads as: “Biophysical factors determining the effects of forest management\non ecosystem resilience”. If you have further suggestions, let us know.\nAdditionally, we have made revisions throughout the manuscript to further emphasize the\nbiophysical factors (i.e., we changed the original “context/background factors” to “biophysical\nfactors” or “biophysical background/context”), such as Abstract (Lines 19~20), Introduction (Lines\n98~99 and Lines 109~111), Figures’ caption (Line 254, Line 294), and Discussion (Lines 411~415).\nHere is the related reversion:\nIn Abstract and Introduction Section:\n“These differences are influenced by biophysical contexts such as local climate, vegetation, and\nsoil.” (Lines 19~20)\n“(2) explore the associated biophysical factors (e.g., climate, soil, and vegetation conditions)\naffecting resilience differences between PF and NF by machine learning” (Lines 109~111)\nIn Result Section:\n“Fig. 3 Relationships Between Resilience Differences (PF vs. NF) and Biophysical Factors.” (Line\n254)\n“Fig. 4 Analyzing the Impact of Biophysical Factors on Vegetation Resilience Differences Between\nPF and NF.” (Lines 294~295)\n“The biophysical contexts including climate, vegetation, and soil variables may affect the resilience\ndifference between PF and NF (details and models for other paired types in Discussion). We further\nuse biophysical factors (and factors difference between PF and NF) to predict the resilience\ndifference between PF and NF by the Extreme Gradient Boosting model (XGBT) (Fig. 4).” (Lines\n265~269)\nIn Discussion Section:\n“Additionally, aridity is not the sole biophysical factor influencing the resilience differences\nbetween PF and NF. Lower temperatures, higher aboveground and belowground biomass, and\nhigher cation exchange capacity (e.g., indicating higher fertility) all contribute to reducing the\nresilience gap (Fig. 4).” (Lines 411~415)\n2. Besides, I have some general major concerns.\nComparing planted forests and natural forests seems odd for me, as the difference in their resilience\nis quite obvious from their management practices. The results stated in the abstract are quite obvious:\nFor example, “forest resilience suffers greater losses under intensive management, as seen in human\nplantations like agroforestry, which are much less resilient than unmanaged natural forests”. And\n“Such resilience difference between PF and NF is related to background conditions such as local\nclimate, vegetation, and soil”. Topography also plays a role?\nR: We respond to the reviewer’s comment by addressing separately the two questions.\nQuestion 1: why do we compare the resilience difference between PF and NF?\nResults show that planted forests can surpass the resilience of natural forests. We indicated that\nnatural forests include both managed and unmanaged ones. So, under certain conditions, the\nresilience of planted forests can exceed that of natural forests. Forest management can cause large\ndamage to the resilience of natural forests without management (i.e., intact forests) (Fig .2). In\naddition, PF here are defined as trees with a longer rotation period (over 15 years) and for restoration\npurposes. Considering these reasons, the resilience of planted forests may indeed exceed that of\nnatural forests. We also explain the mechanism from the perspective of climate sensitivity (see\nDiscussion Section, Lines 427~448).\nBased on previous studies and our findings, it seems difficult for planted forests to surpass the\nresilience of intact forests (Fig. 2). Our research indeed shows that unmanaged forests (i.e., natural\nforest without management) possess the highest forest resilience due to their high biodiversity and\ncomplexity (Fig. 2).\nIn the Abstract Section, such a description is just an example of \"Forest resilience experiences\ngreater losses under intensive management.\" Perhaps a more appropriate example here would be\n\"as seen in human plantations like agroforestry, which are much less resilient than planted forests\nfor restoration.\" After a discussion with our authors, we have finally deleted the sentence.\nAdditionally, we further describe the purpose of this study. Our research differs from that of\nForzieri et al.15, as we further detail the differences in resilience across different forest management\ntypes (not only focusing on the intact and non-intact forests). We found that forests subjected to\nintensive human managements such as rotation and agroforestry exhibit significantly reduced resilience.\nHowever, forests planted with the goal of restoration (i.e., PF with longer rotating years) do not\nnecessarily have lower resilience than natural forests; this depends on their biophysical contexts\nsuch as climate, vegetation, and soil factors. Our results also indicated the regions where forest\nmanagement has a lesser impact on resilience, serving as a policy reminder.\nAction taken. We added the reasons why we mainly focus on the resilience between PF and NF\n(Lines 708~713).\n“We first explored the biophysical drivers for resilience differences between PF and NF. This is\nparticularly important given the widespread attention to restoration activities21, 22, 23, such as the UN\nDecade on Ecosystem Restoration (2021–2030). Understanding the resilience differences between\nPF and NF can provide valuable guidance for designing and implementing effective restoration\nstrategies to enhance ecosystem resilience.”\nQuestion 2: Topography also plays a role?\nThanks for highlighting the potential role of topography. The topography is including factors such\nas landforms, elevation, slope, and aspect. In this study, we have taken steps to mitigate the impact\nof elevation on the estimation of resilience differences. Specifically, to control for the potential\nconfounding effects of varying forest age and elevation, we ensured that the differences in these\nvariables between paired forests did not exceed 10 years and 500 meters, respectively. Here are the\nwords from the original texts in Lines 687~691 (Method Section): “To mitigate the potential effects\nof varying forest age and elevation on vegetation resilience, we ensured that the differences in these\nvariables within each paired grid were restricted to within 10 years and 500 meters, respectively59.\nThe global forest age data are obtained from Besnard et al.88, while elevation data were derived\nfrom the Shuttle Radar Topography Mission (SRTM) v4 dataset89.”\nAction taken. In the Result Section, we re-emphasize this point. (Lines 159~165)\n“Specifically, grids are paired using the sliding window method, which moves a fixed-size window\nacross the study area. Within each window, paired grids are selected based on similar forest age\n(differences of ≤10 years), elevation (differences of ≤500 meters), and a minimum of 20% forest\ncoverage for each management type (more details in the Method Section). We point out that the\neffects of potential covariates have been ruled out within the sliding window approach.”\nAdditionally, while we accounted for elevation differences in our analysis, we acknowledge that\nbroader topographical features—such as slope, aspect, and landforms—may also influence\nresilience differences between PF and NF on a global scale. These factors can affect soil water\navailability, microclimate, and vegetation growth, potentially interacting with forest management\npractices. In response to your comment, we modified and added related limitations in the discussion,\nas follows: (Lines 514~517)\n“This study prominently examines the differences in resilience between PF and NF, further focusing\non their response to biophysical factors. Certain factors (e.g., landforms and slope), that influence\nthe vegetation recovery process, may not have been adequately accounted for in the model70.”\n3. These are obvious well documented results and do not provide any novelty, even where a\nquantified difference is provided. Is the resilience of well managed planted forests at the same level\nas or even more than unmanaged natural forests?\nMy suggestion is to compare similar forest ecosystem types (e.g., managed and unmanaged natural\nforests; planted commercial forests and planted forests mainly aimed to restore ecosystem functions\nand service, etc.). Comparing similar forest ecosystem types with different or similar forest\nmanagement practices and biophysical conditions would be more intriguing in understanding\ndifferences in their resilience and inform restoration and conservation decision making. Comparing\nwithin planted forests with primary objective for restoration functions and services and planted\nforests under the commercial forestry approach can enhance our understanding of the resilience of\nthese ecosystems in terms of supply of ecosystem services. The resilience of planted forests\nmonocultures and multiple tree species forests can be compared to further understand why forest\nmonocultures are discouraged.\nR: We agree with your comments that comparing the resilience of similar forest types is essential.\nWe added the additional analysis to compare the resilience difference between similar forest types,\nincluding managed versus unmanaged natural forests and long-rotation versus short-rotation planted\nforests (i.e., planted commercial forests and planted forests mainly aimed to restore ecosystem\nfunctions) (Fig S24-S29). In addition, we investigate the relationship between tree species richness\nand resilience, to further validate why forest monocultures are discouraged (Fig S24). Here, we\nfurther respond to the reviewer’s comment by addressing separately the two questions.\nQuestion 1: “Is the resilience of well-managed planted forests at the same level as or even more\nthan unmanaged natural forests?”\nAlthough we responded to this question briefly in the previous comment, here we will express it in\ndetail. After many pairwise comparisons, our findings indicate that unmanaged natural forests\nexhibit the highest resilience among all forest management types potentially due to their complex\nand integrated ecological process (Fig. 2 and Supplementary Fig. S26). Previous studies have\nsuggested that management practices such as irrigation, rotation, and logging reduce ecosystem\ncomplexity and biodiversity, thereby diminishing the forests' capacity to recover from disturbances.\nOur findings also indicated that greater intensity of forest management corresponds to reduced\nresilience. For example, the resilience difference between plantations and NF is larger than that\nbetween planted forests and NF (Fig. 2). The plantations include the planted forests, planted forests\nwith rotating up to 15 years, agroforestry, and oil palms, which have more intense human activities\nor management. Through this comparison, we can further infer those commercial forests (such as\ntimber forests-PF with short rotating years, agroforestry, and oil palms), which are more intensively\nmanaged, have lower resilience than planted forests aimed at restoration. Thus, increased\nmanagement does not make planted forests more resilient than unmanaged natural forests (even\nmakes planted forests less resilient).\nHowever, under specific biophysical conditions (e.g., wet conditions with P/PET > 1.5), planted\nforests (PF) can exhibit higher resilience than natural forests, which both include the managed and\nunmanaged NF. These differences hinge on various biophysical factors including climate, soil,\nand vegetation, which vary significantly across regions. This comparison and investigations inform\noptimal reforestation activities for reducing the resilience between PF and NF. Furthermore,\nalthough planted forests can surpass natural forests in resilience under certain conditions (Fig. 2-5),\nthe resilience of unmanaged natural forests is still larger than planted forests (Supplementary Fig.\nS26). This further implies that human managements may significantly impair the resilience of\nnatural forests.\nIn summary, it is important to distinguish between three definitions: well-managed forests,\nintensively managed forests, and human-managed forests considering the biophysical contexts.\n\"Well-managed\" represents the result of human intervention that is challenging to quantify precisely\n(it's multi-dimensional including environment, economy, and benefits). Our findings only suggest\nthat human managements can generally lead to a reduction in forest resilience. Forests under higher\nlevels of human pressure (i.e., intensive managements) tend to exhibit lower resilience. However,\nwhen biophysical contexts are considered in managements (e.g., choosing humid or cold regions for\nafforestation), planted forests maybe have higher resilience than natural forests. However, a\nnoticeable resilience gap remains compared to unmanaged natural forests. Natural forests without\nmanagement have much higher resilience than that in natural forests with management (Fig. 2).\nQuestion 2: why not compare similar forest ecosystem types? (compare resilience between planted\ncommercial forests and planted forests mainly aimed at restoring ecosystem functions and services)\nWe agree with your suggestions, and we first added the analysis and discussion of comparing\nresilience between planted commercial forests and planted forests mainly aimed at restoring\necosystem functions and services. In the original analysis (Fig. 2), we conducted two comparisons:\nthe resilience difference between plantations and natural forests, and the resilience difference\nbetween planted forests and natural forests. The plantations conclude the planted forests, planted\nforests with rotating up to 15 years, agroforestry, and oil palms. We found that the resilience\ndifference between plantations and natural forests (NF) is larger than that between planted forests\nand NF. This further partly supports that commercially oriented plantations have lower resilience\ncompared to restoration-oriented planted forests.\nFurthermore, we attempted to pair the grids method to directly compare the resilience between\nwood-oriented plantations (refer to the planted forests up to 15 rotating years) and restoration-\nfunction-oriented plantations (refer to the planted forests with over 15 rotating years). Unfortunately,\nthe number of paired grids was limited (only ten to several dozen), which may compromise the\nrobustness of our analysis. We have to consider the number of paired grids to design paired\nexperiments. In this regard, we reclassified the GFMD to ensure that there were enough paired grids\nfor use in the analysis. In other words, the comparison results of different forest management types\nin Figure 2 already consider the number of paired grids available. Other pairing results may not have\nsufficient grids for further analysis (typically requiring more than 500 grids).\nAction taken.\n(1) In response, we added the limitation of the paired grids method in the Method Section (Lines\n700~706):\n“Therefore, we designed six types of comparative experiments to examine the differences in\nresilience across different vegetation management types, including resilience between NF and PF,\nNF and Plantations, NF WOM and F WM, NF WOM and NF WM, and NF WOM and Plantations.\nFor each type of comparison, over 500 paired grids were analyzed to ensure statistical robustness\nand comprehensive coverage of the different management comparision.”\n(2) To provide a more robust analysis, we utilized original crowdsourced data to examine the\nresilience of wood-oriented plantations and restoration-function-oriented plantations\n(Supplementary Fig. S25). By definition29, planted forests with shorter rotation periods (<15 years)\nare often commercially oriented plantations aimed at timber production, whereas those with longer\nrotation periods are typically restoration-oriented planted forests (> 15 years). Our enhanced\ncomparison revealed that: In regions with P/PET < 3.5, planted forests with low rotation years (i.e.,\nwood plantation) exhibit significantly lower resilience compared to those with high rotation years\n(i.e., planted forest for restoration). In regions with P/PET > 3.5, the difference in resilience between\nlong and short-rotation age plantations becomes minimal. Overall, the resilience of long-rotating\nplanted forests is larger than that with short rotation periods. We further found that the resilience of\noil palm plantations is generally lower than that of long-rotation planted forests, again highlighting\nthat commercially oriented plantations have lower resilience compared to restoration-oriented\nplanted forests.\nIn addition, we further investigated the relation among P/PET, tree species richness, and forest\nresilience. We found that resilience in all plantations increases with higher tree species richness.\nThe increase in tree species richness further enhances resilience under similar climate conditions\n(Fig. S26). These findings highlight the importance of encouraging multi-species plantations as a\npolicy recommendation for the enhancement of forest resilience and ecosystem services.\nFigure S26.The relation among forest resilience, P/PET, and tree species richness across various\nmanaged types of forests. (a-f) The AR1 by P/PET in various forest management trees. (g-l) The\nrelation between AR1-LAI and tree species richness across various management forests. (m-n)\nThe AR1-LAI heatmaps illustrate the relationship between tree species richness and P/PET across\ndifferent forest types. Panel (m) shows the AR1_LAI for all forests, while panel (n) focuses on\nplantations. Here we used the original crowdsourcing datasets, long-term LAI resilience, tree\nspecies richness, and the P/PET for statistics. The grey bar is the mean value, and the error bar is\nthe mean ± 2*standard error. The red cross is the median value. All datasets such as forest\nresilience, P/PET, and tree species richness are resampled to 1km. Global tree species richness has\n0.025° spatial resolution, derived from\nhttps://figshare.com/articles/dataset/The_global_map_of_tree_species_richness/17232491.\nAccording to the reviewer’s suggestion, we added above-mentioned discussion of resilience\ncomparison between planted commercial forests and planted forests mainly aimed at restoring the\necosystem. We also added the description of afforestation by multiple tree species forests. Here is\nthe related discussion, now reading as (Lines 472~482):\n“For example, PF with rotations up to 15 years (i.e., wood plantations/shorter-rotation PF) exhibit\nlower resilience compared to those with rotation periods over 15 years (i.e., planted forests for\nrestoration/longer-rotation PF) (Supplementary Fig. S26). These shorter-rotation PF are\neconomically composed of a single species and are often partially harvested at regular intervals66.\nAs a result, species richness and biodiversity in shorter-rotation PF are generally lower compared\nto longer-rotation PF, reducing the capacity to recover from disturbances52. Our study also supports\nthe indications that the increase in tree species richness can further enhance forest resilience\n(Supplementary Fig. S26). It encourages multi-species plantations as a policy recommendation for\nthe enhancement of forest resilience for afforestation.”\nComparing managed and unmanaged natural forests can provided important insight on the effect of\nforest silvicultural management practices on forest ecosystem resilience. Unmanaged old growth\nforest may sequester less carbon due to older age trees compared to young secondary regrowth\nforests. Do unmanaged natural forests need some sought of management intervention to increase\nresilience or should they be kept intact?\nThe resilience of forest ecosystems is not only explained by forest management practices (type and\nintensity) but also more so, by biophysical factors. But it is clear this was not investigated in this\nstudy, making the scope so one dimensional. This is clearly acknowledged by the authors. For\nexample:\nIn line 446 to 448, the authors clearly stated that “…the narrowing disparity in resilience between\nPF and NF may be attributed to multiple factors, including vegetation physiology, climate change,\nand the intensity of human activities. This complex issue merits further investigation”.\nIn lines 459 to 460, “… but there is less exploration of whether the climate, soil, and vegetation\nbackground can mitigate the effects of other forest managements”.\nR: We agree with comments on the importance of comparing forest ecosystems under similar\nmanagement practices. According to the reviewer’s comment, we built additional four XGBT\nmodels for investigation of the mitigation of the effects of forest management by the climate,\nsoil, and vegetation background (Supplementary Fig. S29). Here, we respond to the reviewer’s comment\nby addressing separately the two questions.\nQuestion 1: Do unmanaged natural forests need some sought of management intervention to\nincrease resilience or should they be kept intact?\nWe have employed a paired study design to compare: (1) Managed Natural Forests vs.\nUnmanaged Natural Forests and (2) Managed Forests vs. Unmanaged Forests. The results indicate\nthat managed natural forests exhibit significantly lower resilience compared to unmanaged natural\nforests, highlighting the potential adverse effects of management practices on natural forest\necosystems (Part of Figure 2). The results of our crowdsourced analysis also show similar patterns\n(Supplementary Fig. S26-S27). Thus, for unmanaged natural forests, maintaining their integrity\n(i.e., being kept intact) can enhance their resilience, aligning with previous study30. Unmanaged\nnatural forests have the highest resilience among all forest management types. We should avoid\nthe disruption of pristine ecological processes by human management such as irrigation changing\nthe root zone or trimming reducing biomass31. This disruption could reduce the resilience of\nunmanaged natural forests. Policymakers are designating protected areas to facilitate natural\nregeneration processes.\nPart of Figure2. Further comparisons of resilience are performed across various forest management\ntypes using the same methodologies. These comparisons include NF vs. PF, NF vs. Plantations,\nNF WOM vs. F WM, NF WOM vs. NF WM, and NF WOM vs. Plantations. The notation “X vs. Y”\nindicates the resilience difference, calculated as the resilience of Y minus that of X. Statistical\nsignificance tests are performed to assess whether these resilience differences are statistically\ndistinguishable from zero. Bars marked with an asterisk (*) denote statistical significance at p <\n0.05, while bars labeled with 'NS' denote non-significant differences. (l) A conceptual figure\nillustrates different forest management types (e.g., NF WOM). Abbreviations used in the figure\ninclude: NF-Natural Forests, PF-Planted Forests, NF WOM-Natural Forests without Managements,\nNF WM-Natural Forests with Managements, and F WM-Forests with Managements.\nQuestion 2: What’s the relationship between resilience differences and biophysical factors?\nWe agree that these biophysical factors play a significant role in resilience differences, and we did\nnot fully investigation in original main texts. We have further developed and integrated relevant\nmodels to offer a more comprehensive explanation of our results.\nAction taken. We further built an XGBT model to explore other comparing resilience between NF\nand Plantations, NF WOM (without management) and F WM (with management), NF WOM and\nNF WM, and NF WOM and Plantations. These models allow us to explore the role of biophysical\nfactors, such as tree root depth, climate, and soil conditions, in shaping forest resilience under\ndifferent management practices. For instance, we found that the depth root of trees is the most\nimportant factor for predicting the resilience differences between managed and unmanaged NF.\nTherefore, we should avoid management activities like irrigation in areas with deep roots, as this\nsignificantly disrupts the development of natural forest root zones, reducing their ability to mitigate\ndisturbances by accessing groundwater (thus decreasing resilience). Previous study has also\nreported that deep groundwater roots can alleviate the extent of damage to the Amazon rainforest\nduring droughts32, leading the vegetation recovery soon. Here are the modified texts and figures:\nIn Method Section: (Lines 745~748)\n“In addition, we also established XGBT models to explore other comparisons of resilience between\nNF and Plantations, NF WOM and F WM, NF WOM and NF WM, and NF WOM and Plantations\n(more details in Supplementary materials) (Supplementary Fig. S29).”\nIn Discussion Section: (Lines 514~533)\n“This study prominently examines the differences in resilience between PF and NF, further focusing\non their response to biophysical factors. Certain factors (e.g., landforms and slope), that influence\nthe vegetation recovery process, may not have been adequately accounted for in the model70.\nFurthermore, we also built XGBT models to investigate the influence of biophysical contexts on the\nresilience differences between NF and Plantations, NF WM and NF WOM, NF WOM and F WM,\nand between NF WOM and Plantations (Supplementary Fig. S29). Results reveal that temperature\nis the most important factor in predicting the resilience differences between NF and Plantations, and\nbetween NF WOM and Plantations (Supplementary Fig. S29a, d). The root depth is the most important\nfactor in driving resilience differences between NF WM and WF WOM (or between F WM and NF\nWOM). Thus, we should avoid management activities like irrigation in regions with deep-rooted\nforests, as such activities disrupt the development of natural forests' root zones, reducing their ability\nto mitigate disturbances by accessing groundwater. A previous study reported that deep-rooted forests\ncan alleviate drought-induced damage in the Amazon rainforest71, leading to faster recovery. More\ndetailed analyses and discussions regarding the XGBT model results can be found in the\nSupplementary Materials.\nUltimately, our comparisons across different forest management types offer crucial guidance for\nsustainable forest development. Policymakers can better choose where to perform managements to\nmitigate its potentially negative effects on resilience.”\nIn Supplementary Materials: (Lines 91~114)\n“Supplementary Note 4: The relation between biophysics factors and the resilience difference\namong various managed forests.\nWe further developed four XGBT models to explore the relationships between resilience differences\nand environmental factors (Supplementary Fig. S29). Specifically, the models analyze resilience\ndifferences between NF and Plantations, NF WOM and F WM, NF WOM and NF WM, and NF\nWOM and Plantations. The overall performance of these models, as measured by the KGE (Kling-\nGupta Efficiency), ranged between 0.6 and 0.7 (Supplementary Fig. S29), indicating good\nperformance. Results reveal that higher temperature, the most important variable, may exacerbate\nthe resilience differences between NF and Plantations, and between NF WOM and Plantations\n(Supplementary Fig. S29a, d). In regions with high temperatures (e.g., tropical regions), Plantations\nare often dominated by monoculture cash crops, such as oil palm, which lack the structural and\nbiological complexity compared to local forests9. Furthermore, human activities like logging\nprominently disrupt the integrity and biodiversity of forests, leading to lower resilience in\nPlantations compared to NF WOM or NF WM. Furthermore, management practices in regions with\ndeeper root depth may increase the resilience differences between NF WM and NF WOM (or F WM\nand NF WOM) (Supplementary Fig. S29b, c). Forests with deep-rooted depth often rely on\ngroundwater, which provides a buffer against drought or water stress14. However, human\nmanagements can disrupt this reliance through irrigation15, making forests more sensitive to drought\nor water variability, thereby significantly reducing resilience. Furthermore, in regions with deep root\ndepth (e.g. mountainous regions), root and vegetation structures in NF WOM are often more\ncomplex and diverse. Human management in such regions tends to simplify root systems, degrade\nsoil quality, and reduce biodiversity/complexity16. Overall, to mitigate the negative impacts of\nhuman forest management on resilience, it is essential to implement managements considering\nclimate conditions and vegetation biophysical characteristics (Fig. 4 and Supplementary Fig. S29).”\nFigure S29. Similar to Figure 4, but focusing on predicting the resilience differences for the\nfollowing comparisons: NF vs. Plantation, NF WOM vs. F WM, NF WOM vs. NF WM, and NF\nWOM vs. Plantation. The R2 and KGE in the lower right corner are derived from the 5-fold cross-\nvalidation.\nReferences:\n1. Dakos, V., Van Nes, E. H., d'Odorico, P. & Scheffer, M. Robustness of variance and\nautocorrelation as indicators of critical slowing down. Ecology 93, 264–271 (2012).\n2. Smith, T., Zotta, R.-M., Boulton, C. A., Lenton, T. M., Dorigo, W. & Boers, N. Reliability of\nresilience estimation based on multi-instrument time series. Earth Syst. Dynam. 14, 173–183\n(2023).\n3. Smith, T., Traxl, D. & Boers, N. Empirical evidence for recent global shifts in vegetation\nresilience. Nat. Clim. Change 12, 477–484 (2022).\n4. Smith, T. & Boers, N. Reliability of vegetation resilience estimates depends on biomass\ndensity. Nat. Ecol. Evol. 7, 1799–1808 (2023).\n5. Zhang, Y. et al. Warming and disturbances affect Arctic-boreal vegetation resilience across\nnorthwestern North America. Nat. Ecol. Evol. 8, 2265–2276 (2024).\n6. Blaschke, L. L. et al. Spatial correlation increase in single-sensor satellite data reveals loss of\nAmazon rainforest resilience. Earth's Future 12, e2023EF004040 (2024).\n7. Boers, N. Observation-based early-warning signals for a collapse of the Atlantic Meridional\nOverturning Circulation. Nat. Clim. Change 11, 680–688 (2021).\n8. Pickering, M. et al. Assessing the relationship between forest structural diversity and\nresilience in a warming climate. EGU General Assembly 2024, Vienna, Austria, 14–19 Apr\n2024, EGU24-16017 (2024).\n9. van de Leemput, I. A. et al. Critical slowing down as early warning for the onset and\ntermination of depression. Proc. Natl. Acad. Sci. USA 111, 87–92 (2014).\n10. Hu, Z., Dakos, V. & Rietkerk, M. Using functional indicators to detect state changes in\nterrestrial ecosystems. Trends Ecol. Evol. 37, 1036–1045 (2022).\n11. Dakos, V., Carpenter, S. R., van Nes, E. H. & Scheffer, M. Resilience indicators: prospects\nand limitations for early warnings of regime shifts. Philos. Trans. R. Soc. B 370, 20130263\n(2015).\n12. Hansen, M. C. et al. High-resolution global maps of 21st-century forest cover change.\nScience 342, 850–853 (2013).\n13. Camps-Valls, G. et al. A unified vegetation index for quantifying the terrestrial biosphere.\nSci. Adv. 7, eabc7447 (2021).\n14. Wang, Q., Moreno-Martínez, Á., Izquierdo-Verdiguier, E., Camps-Valls, G. & Zhu, X.\nEstimation of vegetation traits with kernel NDVI. ISPRS J. Photogramm. Remote Sens. 195,\n408–417 (2023).\n15. Forzieri, G., Dakos, V., McDowell, N. G., Ramdane, A. & Cescatti, A. Emerging signals of\ndeclining forest resilience under climate change. Nature 608, 534–539 (2022).\n16. Yu, Z., Liu, S., Wang, J., Wei, X., Schuler, J., Sun, P., Harper, R. & Zegre, N. Natural forests\nexhibit higher carbon sequestration and lower water consumption than planted forests in\nChina. Glob. Change Biol. 25, 68–77 (2019).\n17. Besnard, S. et al. Mapping global forest age from forest inventories, biomass and climate\ndata. Earth Syst. Sci. Data 13, 4881–4896 (2021).\n18. Jarvis, A., Reuter, H. I., Nelson, A. & Guevara, E. Hole-filled SRTM for the globe Version\n4, available from the CGIAR-CSI SRTM 90m Database. (2008).\n19. Scheffer, M. et al. Early-warning signals for critical transitions. Nature 461, 53–59 (2009).\n20. Liu, Y., Kumar, M., Katul, G. G. & Porporato, A. Reduced resilience as an early warning\nsignal of forest mortality. Nat. Clim. Change 9, 880–885 (2019).\n21. Zhong, Z. et al. Higher sensitivity of planted forests' productivity than natural forests to\ndroughts in China. J. Geophys. Res. Biogeosci. 126, e2021JG006306 (2021).\n22. Fernández, M.E., Gyenge, J. & Schlichter, T. J.T. Water flux and canopy conductance of\nnatural versus planted forests in Patagonia, South America. Tree Physiol. 23, 415–427\n(2009).\n23. Martín-Benito, D. et al. Response of climate-growth relationships and water use efficiency to\nthinning in a Pinus nigra afforestation. For. Ecol. Manage. 259, 967–975 (2010).\n24. Kuuluvainen, T. & Gauthier, S. Young and old forest in the boreal: critical stages of\necosystem dynamics and management under global change. For. Ecosyst. 5, 1–15 (2018).\n25. Bohn, F. J., May, F. & Huth, A. Species composition and forest structure explain the\ntemperature sensitivity patterns of productivity in temperate forests. Biogeosciences 15,\n1795–1813 (2018).\n26. Morin, X., Fahse, L., De Mazancourt, C., Scherer-Lorenzen, M. & Bugmann, H. Temporal\nstability in forest productivity increases with tree diversity due to asynchrony in species\ndynamics. Ecol. Lett. 17, 1526–1535 (2014).\n27. Brockerhoff, E.G. et al. Plantation forests and biodiversity: oxymoron or opportunity?\nBiodiversity and Conservation. 17, 925–951 (2008).\n28. Bauhus, J. et al. Mixed-species forests: ecology and management. Mixed-species forests:\necology and management. Springer; 2017. p. 1–24. https://doi.org/10.1007/978-3-662-54553-\n9_7\n29. Lesiv, M. et al. Global forest management data for 2015 at a 100 m resolution. Sci. Data 9,\n199 (2022).\n30. Liu, X. et al. Enhancing ecosystem productivity and stability with increasing canopy\nstructural complexity in global forests. Sci. Adv. 10, eadl1947 (2024).\n31. Payn, T. et al. Changes in planted forests and future global implications. Forest Ecol. Manag.\n352, 57–67 (2015).\n32. Chen, S. et al. Amazon forest biogeography predicts resilience and vulnerability to drought.\nNature 631, 111–117 (2024)."}, {"round": 3, "reviewer_comments": {"reviewer_4": "The authors have provided a comprehensive explanation of the limitations of the study, in terms of variable input, analysis\nand future directions. They also have responded satisfactorily to my comments, especially presenting the magnitude of\ndifference in the resilience between forest management types. Therefore, I approve of their revision but encourage the\nauthors to consider highlighting the limitations of the study on the remote sensing and disturbance regime effect sides.\nThe limitations on the remote sensing side in influencing the data analysis and results are profound, but can be overlooked\ndue to the global scale nature of the study. The spectral signatures of vegetation types and management conditions undergo\nconstant change in strength and amount, following various natural and human disturbances in forest ecosystems, whether\ncontinuous, abrupt or subtle disturbances. So this in itself creates challenges for remote sensing data processing and\nanalysis, in addition to image distortion challenges. So, I think the results offer a foundation for future studies using more fine\nresolution data such as from the global ecosystem dynamics investigation (GEDI). I encourage the authors to add these\nremote sensing related limitations in the manuscript to inform the readers of their awareness and transparency.\nThe limitations on the effect of specific disturbance type on the resilience of forest management types is also glaring, but\ndwelling on this in this study is unfeasible (it is a topic on its own), as such should be highlighted in the discussion and\nadded as a recommendation for future studies. Different disturbance regimes affect different forest ecosystems differently in\nspatial and temporal scale."}, "author_response": ""}, {"round": 4, "reviewer_comments": {"reviewer_2": "The authors have satisfactorily responded to all my main comments. Before acceptance, I would only suggest the following\nminor changes:\nSupplementary L318-319: \"PF generally possess lower structural complexity and shallower, simplified root systems\", please\nadd source(s) for this statement.\nSupplementary Figures S22, S28...: you added the coefficient of variation, but labelled this as R instead of R2. If I am not\nmistaken, this should be R^2 (R-squared). Please check this and double check that it is kept consistent across all figures.\nSimilarly, the p-values are missing from e.g., Figures S28-S29."}, "author_response": ""}]

Decision letter
[Review submitted as PDF attachment - text not found in document]