Accelerated land surface greening caused by earlier permafrost thawing

Domain: physics
Published: 2025-12-16
Editorial outcome: major_revision

Abstract
Persistent and above-average warming has advanced the start of spring permafrost thawing, intensifying climate warming through carbon feedbacks. However, the extent to which variations in spring permafrost thawing contribute to greening trends in permafrost-affected areas (i.e., increases in vegetation greenness) over time remains unclear, limiting our understanding of the ecological consequences of permafrost degradation. Analyzing 40-year freeze/thaw data and multiple satellite-derived greenness indicators, we identify widespread increases in the sensitivity of spring greenness to spring permafrost thawing based on moving-window analyses, indicating that advances in spring permafrost thawing have played a progressively stronger role in promoting spring greening, particularly in boreal forests and tundra regions underlain by continuous permafrost. In addition to the region-specific climate and permafrost conditions, we uncover biogeophysical pathways accounting for the increase in sensitivity of spring greenness to spring permafrost thawing, including reduced albedo, earlier vegetation phenology, and enhanced soil moisture infiltration. Notably, state-of-the-art Dynamic Global Vegetation Models consistently underestimate both the magnitude and variability of sensitivity of spring greenness to spring permafrost thawing. These findings highlight the temporal changes in vegetation responses to freeze-thaw dynamics, necessitating improved model projections concurrent with climate change.

Introduction
Permafrost regions, which cover ~25% of the Northern Hemisphere, are undergoing substantial greening over recent decades due to warming temperatures, manifested as increased vegetation greenness, biomass, productivity, and growing season length1–3. Greening over permafrost regions profoundly influences carbon cycling and energy fluxes, with far-reaching implications for global climate systems4–6. Concurrently, persistent and above-average warming has accelerated permafrost degradation, contributing to positive climate feedback through enhanced soil carbon release (i.e., permafrost carbon feedback)7. This degradation also affects vegetation dynamics by altering regional energy and water balance8,9, hydrothermal properties of the soil10,11, and nutrient mineralization12,13. Despite these advances, the relationship between permafrost changes and greening remains poorly understood, partly due to a research focus on soil warming and deepening active layers, which often overlooks the timing of freeze-thaw cycles—an essential indicator of permafrost-climate interactions13.

Warming-induced advances in the start of spring permafrost thawing (SOT) have been documented through ground-based measurements14 and satellite observations15,16, revealing consistent trends but varying magnitudes. Here, SOT refers to the timing when the seasonally frozen active layer above permafrost begins to thaw, serving as a key indicator of land surface changes and permafrost feedback mechanisms17. Acting as an aquitard, permafrost retains soil moisture near plant roots, supporting nutrient cycles essential for plant growth and productivity18. Consequently, an earlier SOT may facilitate vertical and horizontal soil water movement, enhancing root-level water and nutrient availability, and promoting spring greening through earlier green-up dates (GUD) and increased leaf area12,16. However, under warming and drying conditions, earlier SOT has also caused surface soil drying and nutrient depletion, irreversibly accelerating desertification in alpine ecosystems such as the Qinghai-Tibet Plateau19. In addition to the spatially non-uniform responses of vegetation growth to SOT changes, the impacts of SOT are likely to be temporally non-stationary due to rapid anthropogenic climate change and vegetation acclimation20,21. Climate warming may destabilize the interactive relationship between permafrost and vegetation by altering landscapes and soil hydrology13,22. Permafrost greening cannot be fully explained by climate warming alone; soil and hydrological changes associated with permafrost dynamics are also critical contributors21,23–26, highlighting the urgent need to investigate the role of SOT shifts in shaping vegetation trends over time.

To address this, we quantify changes in the sensitivity of spring greenness to SOT (\beta_{SOT}) from 1982 to 2021 for different permafrost and vegetation types over permafrost regions in the northern hemisphere (Supplementary Fig. 1) using a moving window approach. By integrating daily freeze/thaw data, satellite-derived greenness metrics, and a suite of environmental variables (Methods; Supplementary Table 1 and 2), we confirm the advancements in SOT and permafrost greening trends in spring (Supplementary Figs. 2 and 3), identify regions with increased and decreased \beta_{SOT}, and explore the mechanisms driving these variations. Employing machine learning techniques, causal inference analyses, and structural equation models, we analyze the temporal variability of \beta_{SOT}, identify the spatial patterns and environmental pathways contributing to these changes, and assess the ability of Dynamic Global Vegetation Models (DGVMs) to replicate the observed trends (Methods; Supplementary Fig. 4).

Results
Increasing sensitivity of spring greenness to SOT
Advances in SOT were significantly associated with increased spring greenness (April to June) across different vegetation types in Northern permafrost regions during 1982–2021. Based on ridge regression, the estimated sensitivities (\beta_{SOT}) were predominantly positive, with spatially averaged values of 14.8, 2.4, and 3.6 ‰ day–1 for kernel normalized difference vegetation index (kNDVI), solar-induced chlorophyll fluorescence (SIF), and gross primary productivity (GPP), respectively (Supplementary Fig. 5), suggesting that earlier thawing tends to correspond with enhanced spring vegetation activity. Some regions also exhibited negative \beta_{SOT}, possibly reflecting locally constrained responses to earlier thawing. We then determined the time series of \beta_{SOT} using 15-year moving windows from 1982 to 2021 at each pixel and found that 38.6% of regions exhibited a significant increase in \beta_{SOT} for kNDVI, while 22.0% showed a decreasing trend (p<0.05) (Fig. 1A). Comparable patterns were found for SIF (38.2% increase; 24.4% decrease, Fig. 1D) and GPP (44.8% increase; 17.3% decrease, Fig. 1G). The strongest increasing trends in \beta_{SOT} were concentrated in continuous permafrost regions, and diminished with decreasing permafrost continuity (Supplementary Fig. 6). Temporal changes in spatially averaged \beta_{SOT} confirmed significant increasing trends, with slopes of 0.70, 0.17, and 0.35 ‰ day–1 window–1 for kNDVI, SIF, and GPP (Fig. 1B, E, H). By vegetation type, Boreal Forests and Tundra regions exhibited the most pronounced increases in \beta_{SOT}, consistent across kNDVI, SIF, and GPP (Fig. 1C, F, I). To test the robustness of these results, we performed additional analyses using alternative moving window sizes and several other greenness indicators from complementary datasets beyond those presented in the main text. The results remained consistent, differing only in magnitude (Methods; Supplementary Fig. 7). In addition, a spatial consistency analysis across six greenness datasets showed that 77.2% of pixels exhibited consistent directional responses, further supporting our findings (Supplementary Fig. 8). As another validation, we applied a log-log regression to account for possible nonlinear SOT-greenness relationships, alongside a standard multiple linear regression for comparison (Methods). Both approaches consistently confirmed the widespread increases in \beta_{SOT} (Supplementary Fig. 9). Taken together, these robustness checks demonstrate an increasing and stimulating impact of SOT on permafrost vegetation greening in spring over the past four decades.

Spatial attribution analysis of \beta_{SOT} trends
We quantified the contributions of permafrost, climate, hydrology, vegetation and soil nutrition factors in explaining the spatial variability of \beta_{SOT} trends, using a random forest (RF) algorithm combined with SHapley Additive exPlanations (SHAP) values (Methods), with \beta_{SOT} trend derived from ridge regression based on kNDVI as a representative indicator. The RF model explained, on average, 81% of the spatial variability in \beta_{SOT} trends (R2=0.81) (Fig. 2G). Among the environmental factors, mean SOT, mean temperature, and temperature trend had strong association in \beta_{SOT} trends, emphasizing the combined influence of baseline thaw timing and climate warming in driving spatial variations in \beta_{SOT} trends (Fig. 2H). Grouping pixels into different vegetation types overall confirm the importance of background permafrost and climate conditions in accounting for the spatial patterns of \beta_{SOT} trends (Supplementary Fig. 10). The six most important environmental factors, ranked by their absolute SHAP values, were selected to explore regional differences in \beta_{SOT} trends. Comparing regions with increased and decreased \beta_{SOT}, we found distinct differences in these environmental characteristics. For permafrost factors, regions with decreased \beta_{SOT} experienced earlier thawing by ~10.7 days (Fig. 2A) and later freezing by about 6.9 days (Fig. 2D) compared to regions with increased \beta_{SOT}. These results suggest that in severely degraded permafrost regions with extended non-frozen periods, vegetation shows weaker sensitivity to thaw timing, as reflected by attenuated \beta_{SOT}. For climatic factors, regions with increased \beta_{SOT} exhibited faster warming rates (0.02°C year−1), despite having mean temperatures that were 2.04°C lower than decreased \beta_{SOT} regions (Fig. 2B, C). Lower shortwave radiation in increased regions (2.2 Wm−2 less on average, Fig. 2F) suggests greater reliance on non-radiative energy inputs enabled by earlier thawing to meet vegetation energy demands. Greening trends also influenced \beta_{SOT} trends, with regions of increased \beta_{SOT} characterized by slower greening rates (Fig. 2E), suggesting that inherent vegetation characteristics may modulate the sensitivity of greenness to thaw timing.

Biogeophysical processes mediating the SOT-greenness relationship
While the spatial patterns of divergent \beta_{SOT} trends are relatively well characterized, understanding the mechanisms driving their temporal evolution at a given location remains more challenging. We hypothesized that the temporal changes in \beta_{SOT} are mediated by (i) energy balance processes through changes in albedo, (ii) phenological dynamics through shifts in GUD, and (iii) water supply through variations in soil moisture and runoff. To test these hypotheses, we applied causal inference and partial correlation analyses to examine how variations in SOT propagate through intermediate processes, i.e., albedo, GUD, and the difference between root-zone and surface soil moisture (\DeltaSM), that mediate its influence on spring greenness. We then used structural equation modeling (SEM) to quantify the individual and combined contributions of these mediating factors during two separate periods (1982–2001 and 2002–2021).

Using a causal inference analysis (i.e., convergent cross mapping, or CCM) for the entire period (1982–2021), we identified significant causal links between \beta_{SOT} and albedo, GUD, and \DeltaSM (Fig. 3A–C). The overall CCM results, represented by the red lines, confirmed causality for all three variables, with Spearman correlation coefficients (rho, Methods) reaching 1 (p=0.00). GUD exhibited the highest proportion of pixels with causality (65.9%), followed by albedo (65.6%) and \DeltaSM (65.3%). Partial correlation analyses for the two periods revealed that the roles of albedo, GUD, and \DeltaSM as intermediaries linking SOT to greenness differed between 1982–2001 and 2002–2021 (Fig. 3D–F). SOT-albedo correlations were consistently positive in both periods (~33%) but showed a shift in albedo-kNDVI relationships for two periods, with significantly negative correlations increasing from 48% to 54.5% (p<0.05), indicating a growing negative mediation by albedo in recent decades. GUD maintained stable positive correlations with SOT (~12%) but exhibited a significant increase in negative correlations with greenness, rising from 13.1% to 25.6%, reflecting its intensified role in modulating spring greenness. \DeltaSM consistently showed negative correlations with SOT (~11%) but experienced a notable rise in positive correlations with kNDVI, from 4.9% to 8.9%, underscoring its increasing importance in supporting vegetation water demands. SEM analyses overall confirmed divergent regulations of SOT on greenness via albedo, GUD, and \DeltaSM for two periods (Fig. 3G–I). For regions with decreasing \beta_{SOT}, we found that shrub-dominated biomes with deeper root systems tended to experience more negative \DeltaSM trends (Supplementary Fig. 11), suggesting that early thaw may limit water accessibility in deeper soil layers, thereby weakening the vegetation response to SOT.

Evaluation of DGVM performance in reproducing \beta_{SOT} trends
We further tested whether current state-of-the-art ecosystem models can reproduce the observed increasing \beta_{SOT}. Using the 16-model ensemble mean vegetation productivity from the “Trends and drivers of the regional scale terrestrial sources and sinks of carbon dioxide” (TRENDY) project, we found that the proportions of regions with significantly increased (32.2%) and decreased (28.8%) \beta_{SOT} were comparable, indicating no clear spatial dominance in trend direction (Fig. 4A). The models predicted an overall increase in \beta_{SOT} based on the moving-window approach (slope = 0.09 ‰ day–1 window–1, p=0.00, Supplementary Fig. 12), but the magnitude of the increase was 87.1% weaker than the observed results of kNDVI. The high standard deviation of \beta_{SOT} among models, particularly in the continuous permafrost zone (Fig. 4B), highlights substantial inconsistencies in model simulations for permafrost regions. Analysis of individual model outputs showed that, with the exception of ORCHIDEE (Organizing Carbon and Hydrology in Dynamic Ecosystems) and YIBs (Yale Interactive Terrestrial Biosphere), most TRENDY models did not fully capture the increase-dominant trends in \beta_{SOT} (Fig. 4C). Spatial consistency analysis at pixel level revealed that approximately half of the simulated trends differed from observations (Fig. 4D), reflecting variability and limitation in model projections across regions.

Methods
Determination of permafrost thawing and frozen dates
We determined the permafrost thawing and frozen dates using three independent datasets (Supplementary Table 1): (1) a microwave-based freeze-thaw Earth system product, (2) a reanalysis-based skin temperature product, and (3) in-situ observations.

The Freeze-Thaw Earth System Data Record (FT-ESDR) provides daily freeze-thaw classifications derived from 37 GHz passive microwave sensors (SMMR, SSM/I, SSMIS) at 25 km resolution43. Each grid cell is assigned one of four surface states: frozen (0), thawed (1), transitional (2), or inverse-transitional (3). Quality control steps were applied to remove erroneous signals near open water bodies or complex terrain.

We used the ERA5-Land skin temperature product from the ECMWF ERA5-Land Daily Aggregated dataset44. This dataset provides daily aggregates of hourly skin temperature at ~9 km spatial resolution. Each day was classified into frozen (<0°C), thawed (>0°C), or transitional (0°C) states, aligned with the FT-ESDR classification scheme.

We compiled daily surface temperature records from 39 active-layer monitoring sites within the Global Terrestrial Network for Permafrost (GTN-P, https://gtnp.arcticportal.org/), spanning 2001–2021 and representing major permafrost regions (Supplementary Table 5). These site-level data were processed using the same classification thresholds as ERA5-Land to ensure comparability.

For each dataset, we applied a threshold-based moving-window algorithm widely used in freeze-thaw studies45,46. Specifically, the start of spring permafrost thawing (SOT) was defined as the first day within a 15-day window in which ≥12 days were classified as thawed or transitional, i.e., indicating the onset of persistent thaw conditions while minimizing the influence of short-term noise. The start of autumn permafrost frozen (SOF) was defined similarly as the first day within a 15-day window in which ≥12 days were frozen, denoting the return of frozen conditions in autumn47. Sensitivity tests confirmed that window size slightly altered absolute dates but did not affect interannual variability16. The analysis was constrained to the January–June period to capture spring thawing dynamics. To improve robustness, we compared SOT estimates from FT-ESDR and ERA5-Land and found consistent spatio-temporal patterns (Supplementary Fig. 2). Therefore, the averaged SOT and SOF from both sources was used in all analyses. Site-based SOTs from GTN-P were used for validation (Supplementary Fig. 2).

Indicators of spring greenness
Satellite-based greenness indicators. To quantify spring vegetation greenness, we used multiple satellite-based indicators for the permafrost region of the Northern Hemisphere during 1982–2021. As spring phenology in permafrost regions typically begins in mid-April16, we defined spring as starting from April to June. The indicators include the kernel normalized difference vegetation index (kNDVI), SIF, and GPP, which are widely used as the proxies of vegetation greenness and photosynthetic activity of vegetation48–50. We used spring-averaged values to represent vegetation activity from different perspectives and scales. It should be noted that greenness indicators differ in their ecological meaning (e.g., greenness vs. productivity), and while collinearity exists, they are not fully synchronous5,51,52. Caution is therefore needed in interpretation.

kNDVI (Eq. 1) is a nonlinear transformation of NDVI that enhances sensitivity to vegetation greenness by reducing soil background noise and mitigating saturation in dense canopies48. The NDVI data were derived from harmonized AVHRR-MODIS observations with sensor calibration, orbital correction, and machine learning-based fusion53. LCSPP (Long-term Continuous SIF-informed Photosynthesis Proxy) is a SIF-based proxy for photosynthesis reconstructed by training neural networks on OCO-2 SIF and bias-corrected AVHRR reflectance. It shows high consistency with both satellite SIF and site-level GPP54,55. GPP was obtained from BESS v2.0 (the Breathing Earth System Simulator), a satellite-driven model integrating coupled carbon, water, and energy processes. BESS shows strong agreement with FLUXNET data and provides reliable global GPP estimates56.

kNDVI = tanh(NDVI^2) (1)

To enhance robustness, we further incorporated additional greenness-related datasets with shorter temporal coverage, derived from different sensors and/or reconstruction approaches: (1) GOSIF is a global 8-day 0.05° SIF product reconstructed from OCO-2 SIF, MODIS reflectance, and reanalysis data, with strong agreement to FLUXNET GPP across ecosystem57. (2) CSIF provides spatially continuous SIF estimates at 4-day 0.05° resolution using neural networks trained on MODIS reflectance and OCO-2 SIF, capturing realistic spatiotemporal photosynthetic dynamics under clear- and all-sky conditions58; (3) VOD (VODCA v2) is a harmonized microwave-derived daily 0.25° vegetation optical depth product59, reflecting canopy water content by integrating multiple sensors and frequencies through the Land Parameter Retrieval Model. All datasets consistently demonstrated an increasing sensitivity of spring greenness to SOT (Supplementary Fig. 7 and 8), confirming that our conclusions are primarily driven by biogeophysical effects rather than potential biases in observational signals. Full details of all datasets are provided in Supplementary Tables 1 and 2.

Process-based model projections. We utilized model-simulated GPP from 16 terrestrial biosphere models provided by the “Trends and drivers of the regional scale terrestrial sources and sinks of carbon dioxide” (TRENDY) project to analyze vegetation sensitivity to SOT (Supplementary Table 4)60,61. TRENDY is an ensemble framework that estimates ecosystem carbon dynamics under varying initial forcing conditions. We used outputs from simulation scenario S3, which incorporates temporal changes in land use, climate, and CO2 concentration, allowing for an integrated assessment of human and climate-driven influences on vegetation. Model outputs, originally at varying spatial resolutions, were re-gridded to 0.1°×0.1° and aggregated to calculate spring GPP. These simulations were systematically compared with satellite-based observational data to evaluate their ability to capture temporal variability of SOT effects, providing insights for improving model representations of permafrost-vegetation-climate interactions.

To compare spatial consistency of \beta_{SOT} trends, we established an observational baseline using a majority voting method with three satellite-derived datasets. For each pixel, if two or more datasets indicated an increase (or decrease) in \beta_{SOT}, the pixel was classified as increasing (or decreasing).

Climatic and ancillary data
To investigate the relationship between SOT and spring greenness, we utilized a comprehensive range of climatic and ancillary datasets spanning 1982–2021. For consistency across analyses, all datasets were re-gridded to a spatial resolution of 0.1° using the bilinear interpolation method.

Climatic data. Monthly climatic data for 1982–2021 were derived from TerraClimate62, a dataset of monthly climate for global terrestrial surfaces at a spatial resolution of 1/24°. We calculated the average 2-m temperature (Temp), 2-m minimum temperature (Tmin), precipitation (Prec), and shortwave radiation (Srad) from this dataset for spring (April to June). Global averaged marine surface CO2 was obtained from the Global Monitoring Laboratory (GML) of the National Oceanic and Atmospheric Administration (NOAA). These variables served as controlling factors to analyze the sensitivity of spring greenness to SOT. Further details on the data are provided in Supplementary Table 1.

Hydrological and surface energy data. To analyze hydrological and energy processes linked to SOT, we used ERA5-Land reanalysis data (0.1°×0.1°)44, including surface albedo, snow depth, surface runoff, sub-surface runoff, and volumetric soil water. Snow depth during the cold season was averaged from November to April. Surface soil moisture was derived from layer 1 (0–7 cm depth), while sub-surface soil moisture was calculated as the weighted mean of layers 2 (7–28 cm) and 3 (28–100 cm)63. The vertical soil moisture difference (\DeltaSM) was defined as sub-surface soil moisture minus surface soil moisture, representing vertical gradients. Similarly, the vertical runoff difference (\DeltaRunoff) was calculated as sub-surface runoff minus surface runoff to capture vertical hydrological dynamics.

Vegetation phenology and structure data. Vegetation phenology, including spring GUD, peak of season (POS), and date of foliar senescence (DFS) during 1982–2021, was derived from the GIMMS NDVI3g+ dataset64. Pixels with annual NDVI values below 0.1 were excluded to eliminate sparsely vegetated areas65, and the NDVI time series was smoothed using a Savitzky-Golay filter66. GUD and DFS were estimated as the average of two widely used methods to reduce uncertainty: (1) a dynamic-threshold method, where GUD (or DFS) was defined as the date when the scaled NDVI ratio crossed 0.5 during the spring (or autumn) transition; (2) a double-logistic fitting method, where GUD and DFS corresponded to the date of maximum (or minimum) rate of NDVI change67,68. POS was defined as the day when NDVI reached its annual maximum during the summer growing season69. The maximum root length at a grid of 0.0083° was obtained from the EartH2Observe70. Canopy height data were obtained from the ETH Global Canopy Height 2020 product71, which integrates Global Ecosystem Dynamics Investigation (GEDI) and Sentinel-2 observations to provide a global canopy height estimate at a 10 m spatial resolution.

Soil attribute data. Active Layer Thickness (ALT) data were obtained from the Permafrost Climate Research Data Package provided by the ESA CCI Permafrost project (1997–2021; 1 km resolution), derived from a thermal model constrained by satellite observations72. Additionally, soil nutrients, including soil organic carbon (SOC) and nitrogen, were sourced from SoilGrids, a global dataset produced by ISRIC—World Soil Information Center, which employs machine learning to predict soil attributes at six standard depths globally (250×250 m)73.

Land cover data. To minimize anthropogenic and hydrological influences, we excluded croplands, urban areas, and water bodies using the MODIS MCD12C1 land cover product (Collection 6.1). A conservative masking strategy was applied: any pixel classified as cropland, urban, or water in any year from 2001 to 2021 was excluded from analysis. To reduce disturbance-related uncertainty, we masked areas with frequent or severe burning using the MOSEV global burn severity product74. We calculated the maximum dNBR for each January–June period (2001–2021), then excluded pixels with a multi-year mean dNBR >660, effectively removing ~10.4% of the study area from greenness sensitivity analysis (Supplementary Fig. 15).

Permafrost regions and types were delineated using the Circum-Arctic Map of Permafrost and Ground-Ice Conditions from the National Snow and Ice Data Center75 (NSIDC), which classifies permafrost based on extent into four categories: continuous (90–100%), discontinuous (50–90%), sporadic (10–50%), and isolated patches (0–10%) (Supplementary Fig. 1).

Vegetation types were defined according to the Terrestrial Ecoregions of the World (TEOW) framework, which classifies global land areas into 14 biomes and 846 ecoregions (Supplementary Table 3). For our analysis, we focused on stable plant functional types (e.g., tundra, shrublands, grasslands) to assess long-term responses at the biome level.

Estimation of \beta_{SOT}
To quantify the sensitivity of spring greenness to the timing of permafrost thaw (\beta_{SOT}) and its temporal evolution, we applied a moving-window regression approach68. Given the 40-year satellite record, a 15-year sliding window was selected to balance sample size and temporal resolution, with sensitivity tests using 13- and 17-year windows confirming robustness (Supplementary Fig. 7). Within each window, \beta_{SOT} was estimated using three complementary approaches (Supplementary Fig. 4):

Ridge regression. To address multicollinearity among driving variables, including SOT and climate variables, we employed ridge regression76,77 as the primary method for all main-text analyses. This regularization technique penalizes coefficient magnitude to reduce variance inflation, ensuring more stable estimates. In each window, the normalized anomalies of spring greenness indicators (GIs: kNDVI, SIF and GPP) were regressed against SOT, spring mean temperature (Temp), precipitation (prec), shortwave radiation (Srad), and atmospheric CO2 concentration (CO2):

GI_{anl} = C_0 \times SOT + C_1 \times Temp + C_2 \times Prec + C_3 \times Srad + C_4 \times CO2 + RES (2)

Here, “anl” indicates the anomaly (i.e., deviation from the long-term mean), which removes the effects of background values while preserving interannual variations influencing vegetation growth78. C0–C4 represent regression coefficients of each driving factor. RES is the residual.

For interpretability, SOT advancement was defined as positive, and thus \beta_{SOT} was computed as the negative of the SOT coefficient, standardized by the mean GI, to enable comparisons across indicators:

\beta_{SOT} = -C_0 / Mean(GI) (3)

Ordinary least squares regression. As a baseline, we also applied standard OLS regression79 using the same explanatory variables as the ridge model. While OLS is more straightforward and interpretable, it does not account for variable collinearity, which may inflate variance or bias coefficients. OLS estimates, therefore, provide a reference point for evaluating the robustness of ridge-derived results (Supplementary Fig. 9).

Log–log regression. To assess potential nonlinear effects and reduce artifacts from linear assumptions80,81, we further applied a log-log regression with SOT as the sole driving variable:

log(GI) = C_0 \times log(SOT) + RES (4)

Accordingly, \beta_{SOT} was calculated as:

\beta_{SOT} = -C_0 (5)

This approach provides an alternative view of the SOT-greenness relationship under multiplicative scaling and complements the multivariateregression analyses.

In summary, ridge regression was selected as the main method due to its robustness to collinearity, OLS served as a transparent baseline for comparison, and log–log regression offered insights into potential nonlinear scaling. Together, these approaches ensured that our estimates of \beta_{SOT} were not artifacts of a single statistical choice but reflected consistent patterns across methods (Supplementary Fig. 9).

Analyses
Trend analysis. To quantify long-term changes in the sensitivity of spring greenness to thaw timing, we calculated the trend of \beta_{SOT}, which is defined as the Theil–Sen slope of \beta_{SOT} estimates across moving windows from 1982 to 2021 (Supplementary Fig. 4), expressed in ‰ day–1 window–1. This method provides a robust, median-based trend estimator less sensitive to outliers. The significance of the trend was assessed using the non-parametric Mann–Kendall test, with a significance level set at P<0.05.

Spatial attribution analysis. To identify the key environmental factors controlling the spatial variation in \beta_{SOT} trends over 1982–2021, we employed an explainable machine learning framework combining Random Forest (RF) modeling and SHAP. The dependent variable was the \beta_{SOT} trend estimated from ridge regression, while environmental factors were grouped into five categories: permafrost, climate, hydrology, vegetation, and soil nutrition. The permafrost category included SOT, SOF, and ALT; climate factors comprised spring temperature, precipitation, shortwave radiation, snow depth, and winter Tmin; hydrological factors included soil moisture and runoff; vegetation factors covered GUD, POS, DFS, kNDVI, and canopy height; and soil nutrition was represented by SOC and soil nitrogen. For multi-year variables, both mean states and linear trends were included. Details of all variables are provided in Supplementary Table 1.

RF models were trained using these factors. As a data-driven machine learning algorithm, RF does not require statistical assumptions about explanatory variables or target variables, and it leverages large datasets to enhance stability82. While collinearity among variables does not affect the overall model performance of RF, it can influence variable importance rankings, as strongly collinear variables may mask each other’s contributions, reducing interpretability83. To mitigate this issue, we calculated variance inflation factors (VIF) within each category and excluded variables with VIF≥3, including soil moisture, spring Tmin, mean kNDVI, and mean GUD84.

To interpret the RF results, we applied SHAP to isolate the marginal contributions of each environmental factor to the target variable. SHAP values, calculated as absolute weighted average marginal contributions, were used to rank variable importance and identify dominant factors influencing \beta_{SOT} trends63,85. The RF models were fitted using the “ranger” package86, and SHAP values were extracted with the “treeshap” package87 in R (version 4.3.1).

Mechanistic pathway analysis. To better understand how the timing of spring thaw modulates vegetation greenness, we developed a mechanistic framework grounded in ecohydrological and climatological theory. Previous studies have suggested that permafrost thawing may influence vegetation through multiple interacting processes. Earlier thawing alters the surface energy budget by reducing albedo, thereby increasing absorbed solar radiation and accelerating soil warming are generally favorable for early-season plant growth26,28. Thaw timing also affects the onset of vegetation activity; by advancing GUD, it may extend the growing season and enhance carbon assimilation in spring16,88. In addition, changes in thaw dynamics influence subsurface hydrology, potentially increasing root-zone soil moisture through enhanced infiltration and reduced surface runoff, thereby improving water availability during the early growth period17,18. Based on this process-level understanding, we hypothesized that the sensitivity of vegetation greenness to thaw timing (\beta_{SOT}) is regulated through three interconnected pathways: energy balance, phenological timing, and water availability. To test this hypothesis, we employed a suite of complementary approaches.

To disentangle the potential mechanisms underlying \beta_{SOT} trends, we adopted a causal inference framework using CCM, which is a nonparametric approach for detecting causality in coupled nonlinear systems89. Unlike correlation-based methods, CCM tests whether the historical states of a putative cause can reconstruct the time series of an effect, thus enabling causal interpretation in complex ecological systems. This method has been increasingly applied in Earth system science to uncover feedbacks and interactions beyond statistical association90. In this study, CCM was implemented at two spatial scales: (1) Regional level: We calculated the spatial mean of \beta_{SOT}, albedo, GUD, and \DeltaSM within each moving window. CCM was then applied to the resulting time series (1982–2021) to detect overall causal linkages across the study region. (2) Pixel level: CCM was further conducted at each grid cell to obtain spatially explicit patterns of causality. For each CCM test, the predictive skill (\rho) was used to assess causality. A variable was deemed to causally influence \beta_{SOT} if \rho increased and converged with larger library size (L, or step length). To statistically confirm convergence, we calculated the Spearman correlation between \rho and L; a significantly positive correlation (p<0.05) indicated a robust causal relationship. All input variables were processed as anomalies to remove long-term trends, and the time lag was set to zero to focus on contemporaneous causal effects.

We conducted pixel-wise partial correlation analysis to assess their intermediary roles in the SOT-kNDVI relationship, while controlling for the influence of climatic variables. This approach isolates the unique contribution of each intermediary variable by removing the effects of confounding factors67,91. To assess temporal shifts, the analysis was performed separately for 1982–2001 and 2002–2021. Statistical significance for partial correlations was determined at a threshold of p<0.05.

We used SEM to quantify direct and indirect effects of SOT on spring greenness92. The analysis was conducted separately for 1982–2001 and 2002–2021 to assess temporal shifts in pathway strengths based on the three hypothesized pathways introduced above. All variables were first standardized to ensure comparability across regions and units. The models were implemented at the pixel level, where each grid cell was treated as an independent unit. For each pixel, a separate SEM was constructed using the time series of gridded variables, and maximum likelihood estimation was used to derive the standardized path coefficients. After pixel-wise model fitting, we computed the direct, indirect, and total effects by multiplying the standardized coefficients along each path. To characterize general patterns across permafrost regions, we then calculated the spatial average of the pixel-level coefficients and reported the aggregated results in Fig. 3G, H. Model validity at each pixel was assessed using standard SEM fit indices: the χ² test (p>0.05), Comparative Fit Index (CFI>0.9), Standardized Root Mean Square Residual (SRMR<0.08), Goodness of Fit Index (GFI>0.95), and Root Mean Square Error of Approximation (RMSEA<0.08). A model was considered acceptable if at least three of five criteria were met93. All SEM analyses were performed using the “lavaan” package in R94.

Uncertainty analysis. To assess the robustness of our estimated \beta_{SOT} trends, we conducted uncertainty analyses from both spatial and temporal perspectives.

Spatial subsampling analysis: To evaluate the robustness of \beta_{SOT} estimates within each time window, we performed a spatial bootstrap procedure across all 15-year moving windows. Specifically, for each window, we randomly sampled 30% (or 70%) of valid pixels from the study region and computed the spatial mean of \beta_{SOT}. This process was repeated 200 times per window to characterize the uncertainty associated with spatial heterogeneity. Distributions of the resulting estimates were summarized with boxplots to quantify spatial uncertainty. The consistent temporal pattern and narrow uncertainty ranges across iterations indicate that the observed \beta_{SOT} dynamics are not driven by specific pixel configurations (Supplementary Fig. 13).

Temporal subsampling analysis: To evaluate the influence of interannual sampling variability within each 15-year moving window, we implemented a temporal bootstrap approach. Specifically, for each window, we randomly selected 10 out of the 15 years (without replacement) and recalculated pixel-wise \beta_{SOT} using ridge regression. This process was repeated 10 times, producing 10 sets of \beta_{SOT} trend estimates across the 26 windows from 1982 to 2021. For each pixel, we computed the Theil–Sen slope of \beta_{SOT} over time in each replicate and summarized the results by calculating the mean and SD across replicates. The SD map quantifies the sensitivity of trend estimates to temporal sampling variability. The high consistency across replicates confirms the stability of \beta_{SOT} trends over time (Supplementary Fig. 14).

Together, these spatial and temporal uncertainty analyses provide a comprehensive evaluation of the robustness of our findings and support the reliability of the detected \beta_{SOT} patterns and trends.

Peer review rounds
[{"round": 1, "reviewer_comments": {"reviewer_1": "The authors presented a study that offers some insights into an interesting topic related to the spring thawing in permafrost under the climate warming background. The impacts of the advanced starting time of spring thawing on vegetation growth were analysed by using remote sensing observations. Although this study is generally well conducted, some parts may need improvements. My comments are listed as follows.\n1) My main concern is that, in Eq. 4, the sensitivity of spring greenness to the start of spring thawing (SOT), i.e., β_SOT, was determined by a linear regression model. The accuracy or the residual of this model should be presented for confirming the robustness of quantifying β_SOT by this way. Additionally, the four factors in the equation may have some intrinsic connections to each other, such as SOT is directly related to temperature and radiation, temperature is correlated with radiation, CO2 is related to temperature, etc. But if we use the coefficient in a linear model to represent the sensitivity of a certain factors to a target variable, it is better to use factors that are independent to each other. How authors consider this problem when using this approach to quantify β_SOT?\n2) Authors used large amount of datasets and machine learning models to show that β_SOT is mainly driven by SOT, which is shown in Figure 2. I am very curious that, whether the result has a large change if adding more comprehensive influencing factors, such as the minimum temperature in spring in the current year and in winter in the last year. Also, it is interesting to check if there are any impacts from the phenological variables (e.g., peak of growing season, end of growing season, length of growing season, and/or freeze time) in the last year?\n3) From the section of Method, SOT was determined by using a moving window calculation based on remote sensing derived surface temperature data. Is it possible to collect some ground-sourced observations to validate the remote-sensing-based estimation of SOT?\n4) In causal analysis, authors constructed the model structure (includes the influencing variables and their effects’ directions) at first and then used the observed data to fit the model parameters and thereby calculated the magnitudes of those direct and indirect effects. This is heuristic approach to constrain the model structure. I would like to see if authors can provide more supports from previous work to confirm the selected key factor in the model, and the model structure, are valid and sufficiently to be used to reflect the physical, chemical, biological processes of vegetation growth enhancement in response to the warming through the advanced spring thawing and other factor changes.\n5) It might be better to show the trend of β_SOT in a continuous colormap in Figure 1, or to classify the values into several intervals, to display the spatial pattern of β_SOT trend.\n6) In Figure 2, can authors show the direction (positive or negative) of these factors’ effects?\n7) In Figure 1 and 4, a space needs to be added between “Frequency” and “(%)”.", "reviewer_2": "The authors derive global measures of the date of start of near-surface soils thaw over 40 years using various remote sensed data in microwave bands (I no little of such products). Using a variety of graphical and statistical methods, they report on the sensitivity of this date to three greenness indicators, and to other factors climatic, edaphic and hydrological. They also report on linear trends in this sensitivity over time (they call these \"the trend\"). Greenness is derived from multi-spectral imagery in the visible and NIR spectral bands (i.e. presumably independent of the microwave-frequency data). Quite an exemplary assembly of disparate global datasets, I would say, in passing.\nTo my mind, the main results are\n1) the coherent spatial patterns in the areas exhibiting significant trend, especially positive trends, indicating an increasing strength of the relationship between SOT and greenness at high latitudes and/or in areas of continuous permafrost;\n2) the attribution of mechanism to these trends by causal inference and by structured equation modelling, the latter indicating the change in SOT acts through increasing albedo, but only (or by far the most strongly) over the past 20 years when albedo drives greenness.\n3) the determination that ecosystem models / DGVMs generally did not reproduced these findings very well, especially at high latitudes and/or in zones of continuous permafrost where the 16 models in the ensemble evaluated are in poor agreement.\nI think this work could have broad significance to the fields of ecosystem modelling and ecological forecasting on the basis of the 3rd point alone. The structural equation models might further inspire developments in ecosystem modelling especially in boreal and tundra systems, e.g. the importance of phenology and the need to include snow-cover in annual models, as affecting albedo. I am not sufficiently familiar with this global literature to make an informed comparison but the work certainly expands upon the regional analyses from North American that I do know.\nI believe the work does support the conclusion and claims made, barring a few places in the manuscript where I found the methods unclear (see my annotated copy). I found no flaws in the analysis or methodology. The statistical methods may even raise the bar somewhat for future analysis, the use of causal inference methods, and SEMs in particular.\nIf any part of the manuscript were less necessary to the story, I think it would be the random forest model, in regard to which I am curious as to why SOT would have selected as a covariate at all, and why it turned out to be far the most important in explaining variation in the trend in sensitivity. I suppose areas where SOT is latest have the most room to move in response to warming, especially given polar amplification?\nIn my experience, written descriptions are scarcely sufficient to reproduce even the simplest model or analysis. However combined with code on the git repo, the text should do, given the clarifications requested in my annotated ms.\nI have one suggested change in notation that might improve readability. The term \\beta_{SOT} is used as an index of sensitivity...it comes from the estimated coefficient of one or other regression model that's calculated for a pixel over some temporal window, and/or for various regional or global means of the pixel level values. It would hard to remove all ambiguity without more complex notation, but I think the *trend* in such coefficients could more clearly be expressed as ∂β, denoting the partial derivative of β with respect to time, which I think is what it pretty much is.", "reviewer_3": "With their manuscript “Accelerated land surface greening with earlier permafrost thawing” the authors present research results for a most relevant topic associated to climate change and increasing temperatures: What is happening to the vast permafrost regions and how are they responding to the changing environment in terms of vegetation activity and land surface greening. Relying on 40 years of data, and being based on causal analyses, this study contains lots of valuable insights. The following contains some comments about certain aspects of the manuscript. At the end, there is a list of smaller remarks and potential modifications.\nKey results\nThe key results of this study are connected to the influence of SOT (Start of Spring Thawing) on vegetation greening in permafrost regions. It was found that changes in temperature, Active layer Thickness, and precipitation are influencing factors in the question whether earlier SOT leads to enhanced greening. In most permafrost regions, earlier SOT can be associated with positive feedback (more greening), while in some other regions, the opposite is the case.\nValidity\nIn this study, a lot of different datasets from varying sources are coming together. That makes it hard to estimate, which potentially negative effect one single dataset might have on the overall results. The input data comprise satellite-based passive microwave data, which is inherent to problems along the coast lines and in the proximity of water/lakes in general, and is problematic in mountain regions. The different sensors throughout the 40-year time series feature different bandwidths, and I’m not sure whether this could be an influencing factor referring to the quality and consistency of the data. The climate data about precipitation can be questionable as well at times, and this aspect has not been discussed enough from my point of view. ERA5 data about albedo and snow depth for example can be similarly uncertain. The dataset of permafrost ALT is only available from 1997 – 2018, but it is not clear to me how the years before and after were accounted for, or were they skipped from the analysis including ALT? The MODIS land cover product is only available after 2000; The NSIDC permafrost map is a compilation of other existing data sources, which may have varying levels of quality and consistency. The problem that arises is that each dataset introduced a certain error, and in combination, these errors can become bigger. An assessment of this uncertainty is missing.\nSignificance\nThis study provides new insights into the complex relationships between spring thawing, greenness, and climate change in boreal forests and tundra regions. The findings have significant implications for our understanding of the impacts of climate change on these ecosystems, which are critical components of the global carbon cycle and are particularly vulnerable to climate change. The study identifies key bio-geophysical pathways that account for the increases in the sensitivity of spring greenness to SOT, providing new insights into the underlying mechanisms driving these changes. It also reveals that state-of-the-art models consistently underestimate the magnitude and variability of the sensitivity of spring greenness to SOT, emphasizing the need for improved model projections concurrent with climate change. It is therefore a significant study with significant results, although the discussion still lacks addressing the data quality and uncertainty part, which needs to be added.\nData and methodology\nAs stated above in the general validity section, the problem I have with the data is the uncertainty connected to the input data, which differs greatly in its magnitude. I would like for the authors to address this issue, including the following steps: Include an estimation for uncertainty propagation; Quantification of the uncertainties for each input dataset and in their combination; Performing a sensitivity analysis to understand how the results change in response to variations in the input datasets; Discuss the issue openly.\nAnalytical approach\nSome of the applied methods are probably only known to experts from the fields of causal analysis, which is not too often used in the context of arctic greening/ greening in permafrost regions. I’m also not a dedicated expert in this field, which makes it difficult to assess the analytical approach of the results. I would prefer more context to the methods in the materials and methods section, outlining why a certain model or method was chosen and what makes it a suitable method for the presented use case.\nSuggested improvements\nMy main suggestions are connected to the input data uncertainty and quality, and that this topic needs to be addressed, quantified, and discussed. It is very important to make clear that a lot of the input data can be ambiguous and unprecise. Combining several of these datasets can multiply the uncertainty, leading to error propagation which is not addressed at all in the manuscript so far.\nClarity and context\nSometimes the authors are using unclear abbreviations and very subject-specific language, which is hard to follow for non-experts. In some of the figures, abbreviations can be found that have never been introduced. It is sometimes tedious to go back and forth to the methods section and back just to find a terminology that is being used in the main text. Would be good to address that.\nReferences\nThe references appear to be well chosen and representative for this field of research, containing state-of-the-art publications and studies.\nSome additional remarks:\nLine 1: The abstract suggests a causal relation between land surface greening and the dynamics of permafrost thawing. Maybe change the wording to emphasise this link, e.g.: “Accelerated land surface greening caused by earlier permafrost thawing”\nLine 16: It is clear from context that the abbreviation is referring to permafrost thawing, but perhaps it should be added here just once to clarify the abbreviation, e.g.: “the start of spring permafrost thawing”\nLine 18: How are permafrost greening trends different from vegetation greening trends? Does this term refer to vegetation greening trends occurring due to changes of permafrost soils or is it a synonym to Arctic greening?\nLine 68: It seems that “SOT” is missing in this abbreviation.\nLines 81-82: General remark: positive βSOT is defined as “advances in SOT lead to enhanced greenness in spring”, therefore negative βSOT can probably be assumed to correspond to “advances in SOT lead to reduced greenness in spring”. Are there any cases of delayed SOT and how would this be treated in relation to βSOT?\nLines 83-84: The formatting of the percentage changes of the three indices is not quite clear: change of sensitivity of spring greenness to start of spring thaw onset in percentage per 10 days for the entire time period?\nLine 98: It is curious that the trend pattern for the GOSIF dataset is different to all other vegetation indices, with smaller trends for the 15-year window compared to the 13-year window. Why might that be?\nFigure 1: This figure is very informative and visually appealing. Minor improvements in readability cloud be made by increasing the font size of the y-axis of the box plots and the colour saturation of the trends. In addition, the sub-figures C, F and I lack an explanation for the grey colour, which probably refers to no significant trends.\nLines 127-129: The trend of βSOT is either referred to “enhanced” or “weakened”, which probably refers to positive and negative trends?\nLines 137-140: The result that larger vegetation types like shrubs exhibit a weakened trend of spring greenness in relation to earlier SOT is surprising, as shrubs are often linked to stronger vegetation greening. Negative βSOT trends were associated with detrimental soil moisture conditions for herbaceous biomes (Lines 181-184), which indicates that shorter vegetation types would be more prone to weakened spring greenness.\nFigure 3: It seems that the headline in sub-figure I is missing. The method section indicates that the library size in sub-figures A, B and C refers to the sample size. So the library size represents the 25 different 15-year windows and convergence is already achieved after roughly 5 windows?\nLines 395-397: How was the start date of the spring aggregated vegetation indices derived? Recent publications on phenology in the Northern Hemisphere report shifts of start of growing season (SOS) of up to 30 days earlier (e.g. Chen, X.; Yang, Y.; Du, J. Distribution and Attribution of Earlier Start of the Growing Season over the Northern Hemisphere from 2001–2018. Remote Sens. 2022, 14, 2964. https://doi.org/10.3390/rs14132964), which could shift the SOS into March in the permafrost region.\nLines 493-494: Is there a reason why the vegetation type was delineated using the static TEOW dataset instead of the yearly MCD12C1 product? The increase of shrubs and the subsequent decrease of lichen during the 40 year study period might be significant, however perhaps not on the spatial scales of this study of 0.1°.\nLine 565: SHAP was already defined in the above paragraph.\nLines 571-572: “Enhanced” and “weakened” trends are not defined, probably refer to positive and negative trends but could also be interpreted as a qualitative indicator of significance.\nLine 575: Sentence is repeated.\nLines 579-580: How were the cause variables chosen? Based on SHAP values?\nLines 601-603: Was the SEM applied to the spatially averaged variables or pixel-wise? Figure 3 G and H imply that the resulting graph stems from the spatially averaged variables, however this section mentions “gridded variables”.", "reviewer_4": "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_5": "This study investigates the long-term changes in the sensitivity of vegetation greenness to the start of spring thawing (SOT) in permafrost regions and its underlying mechanisms, shedding light on the complex interactions between freeze-thaw dynamics and vegetation activity. However, I think this paper needs significant revise based on the following comments:\n(1)Limitations in Data and Methods\nSatellite data (e.g., 25 km resolution) may miss fine-scale heterogeneity (e.g., microtopography, local hydrology). Moreover, The analysis of data uncertainty is not sufficiently comprehensive. Although the article mentions the consistency of results across different remote sensing indicators, it does not thoroughly discuss how inherent uncertainties from different data sources (such as satellite observations and reconstructed datasets) may affect the research findings. It is recommended that the discussion section be expanded to include an assessment of how uncertainties stemming from various data sources might influence the outcomes of the study.\nThe discussion of the applicability of causal analysis methods is insufficient. Although the article applies causal inference techniques such as CCM, it does not adequately address the limitations of these methods across different spatial scales or ecosystem conditions. It is recommended to further discuss the applicability and potential limitations of these methods to enhance the robustness of the study's conclusions. In addition, while croplands and urban areas were excluded, natural disturbances (e.g., wildfires, insect outbreaks) impacting greenness were not addressed, potentially introducing bias.\n(2)Deepening Mechanistic Explanations:\nInteractions among the three pathways (albedo, phenology, moisture) require further exploration. For instance, reduced albedo may accelerate thaw via energy absorption, indirectly affecting moisture infiltration—a cascading effect not explicitly discussed. Coupling mechanisms between soil carbon-nitrogen cycling and vegetation growth are under-addressed, potentially overlooking critical biogeochemical feedbacks.\n(3)Concrete Model Improvement Recommendations\nClarify specific deficiencies in current model parameterizations of freeze-thaw dynamics (e.g., active layer thickness, three-phase water transitions) and propose actionable solutions (e.g., vertical soil module integration). Explore hybrid approaches combining process-based models (e.g., CLM, LPJ-GUESS) with machine learning to better capture nonlinear relationships."}, "author_response": "We sincerely thank the reviewers for their constructive and detailed feedback. We have\ncarefully revised the manuscript to address all concerns, with changes highlighted in red\nfont throughout the revised version. The major revisions we made are summarized as\nfollows:\n1. Uncertainty Quantification and Data Robustness:\nAs noted by both reviewers 3 and 5, the use of multiple satellite, reanalysis, and ecological\ndatasets raises important concerns regarding input data uncertainty. To address this, we\nimplemented a multi-pronged strategy to quantify and mitigate uncertainty across key\nvariables. We quantified potential uncertainties across key variables—such as freeze–thaw\nstate, vegetation greenness, and climate inputs—by integrating long-term data from\nmultiple sources and supplementing them with ground-based site validation. In addition to\nexcluding croplands, urban areas, and water bodies, we also masked wildfire disturbances\nby removing pixels with high dNBR values derived from the MOSEV dataset. To assess\nthe robustness of our findings, we conducted cross-dataset comparisons using\nindependent vegetation indices, as well as spatial and temporal subsampling analyses to\nquantify sensitivity to sampling design. A dedicated paragraph has been added to the\nDiscussion to openly acknowledge uncertainty propagation in multi-source datasets and to\nemphasize the need for improved ground validation, disturbance monitoring, and high-\nresolution observations in future research on permafrost ecosystems.\n2. Clarification of Analytical Framework and Methods:\nIn response to concerns raised by reviewers 1 and 3, We refined the statistical framework\nto ensure robustness and interpretability. For the sensitivity analysis, we conducted\nresidual diagnostics on the multiple linear regression and replaced all reported β values\nSOT\nin the main text with ridge regression results to account for interdependence among\ncovariates. To improve transparency for a broader readership, we expanded the Materials\nand Methods section to introduce the rationale behind ridge regression and Convergent\nCross Mapping (CCM), emphasizing their suitability for quantifying sensitivity and detecting\ncausality in nonlinear, feedback-rich systems. We clarified how CCM was applied to test\nthe causal influence of albedo, GUD, and ΔSM on β , and how SEM was implemented\nSOT\nat the pixel level and summarized across permafrost regions. In addition, we reviewed and\nstandardized abbreviations and technical terms throughout the manuscript and figures to\nimprove accessibility and consistency.\n3. Deepening mechanistic interpretation and model relevance:\nAs suggested by both reviewers 2 and 5, we strengthened both the conceptual framing\nand mechanistic interpretation of how earlier permafrost thaw enhances vegetation\ngreenness through energy (albedo), phenology (GUD), and moisture (ΔSM) pathways.\nUsing SEM and CCM, we clarified the causal directions and pathway structure linking thaw\ntiming to spring productivity. Beyond SOT, we also identified additional controlling factors—\nparticularly the start of freeze (SOF)—that significantly modulate the spatial variation in\nβ , as revealed by machine learning-based spatial attribution. We proposed targeted\nSOT\nimprovements to Dynamic Global Vegetation Models (DGVMs) by addressing key\ndeficiencies in freeze–thaw parameterizations, including missing representations of active\nlayer thickness, three-phase water transitions, and abrupt thaw events. We recommend\nincorporating vertically resolved soil thermal–hydrological modules, vegetation–hydrology\ncoupling, and hybrid frameworks that integrate process-based DGVMs with machine\nlearning to better capture nonlinear and region-specific responses.\nWe hope these revisions address the reviewers’ comments satisfactorily and improve the\nclarity, rigor, and relevance of our study. Detailed point-by-point responses are provided\nbelow."}, {"round": 2, "reviewer_comments": {"reviewer_1": "I would like to thank authors for carefully considering all my comments and suggestions to the original version of the manuscript.The revised version of the manuscript has appropriately addressed my concerns.", "reviewer_5": "I think the authors have addressed all my concerns.", "reviewer_6": "Both reviewers pointed out the issues from 1) the data quality and uncertainty from the large volume of inputs in this work including satellite observations/derivatives, climate data, and reanalysis, and 2) multicollinearity or high correlation of these datasets. The authors made efforts to address these two major issues by expanding “independent datasets” indicating greenness, spatial and temporal resampling of the datasets, and adding in-situ observations, and also included a multilinear model to assess other drivers in greenness change.\nThe solutions provided to the data quality and uncertainty issues are not robust to some degree (but I do know that such uncertainty from data inputs are difficult to track). Based on the descriptions of each dataset:\nkNDVI (originally used): from AVHRR and MODIS\nLCSPP (not spelled out, originally used): ANN modeled derived from OCO-2 SIF and AVHRR\nGPP (originally used): briefly mentioned from satellite-driven BESS (not spelled out), not sure which satellites are used.\nGOSIF (not spelled out, added dataset): reconstructed from OCO-2 SIF, and MODIS\nCSIF (not spelled out, added dataset): trained from OCO-2 SIF, and MODIS\nThese datasets from original ones or added ones are not completely independent of each other because all have used inputs from AVHRR and MODIS. The temporal and spatial resampling of each dataset cannot address the uncertainty and quality of each dataset and thus it is difficult to track the effects of these datasets on the analysis results at such a large spatial scale, especially the average SOT change is about 4 days per 15 years based upon the largest trend of 0.7‰. not speaking the 0.17‰ (less than one day) per 15 years.\nTo address the multiple driver issue raised by the reviewer, the authors added multilinear regression model displayed in Equation 2, which caused confusion with Equation 4 representing SOT only with data transformation. It seems that the model was developed for each pixel representing location x (each one has 40 years matched dataset) to catch the temporal trend of each location, if so, then for the entire study domain, coefficient will be different for each location or across different ecotype and a map can be derived in equation 4 showing a spatial pattern, thus difficult to understand next step to model spatial pattern using RF model again.\nThere are many unclear abbreviations/terms pointed by previous reviewers and the authors said these have been solved such as “Figure captions have been updated to be fully self-contained and clear without requiring repeated cross-referencing”. This is not correct. For example, LCSPP has not been defined; BFT, TUN, TCF, MGS, TGSS appeared in Figure 1 are not defined.\nThere are many datasets used in this work. Some are used for SOT estimation, some are used for trend estimation, some are used for association analysis/mask out non-interested areas. I would suggest the authors made a workflow about how to link these datasets and their role in this work, and explicitly present the spatial and temporal resolution of each dataset and rationale in selecting these datasets."}, "author_response": "We are very grateful for the reviewers’ constructive comments and for the\nopportunity to further improve our manuscript. We have carefully revised the text,\nfigures, and supplementary materials in response to the remaining concerns, with\nthe goal of making our research design and conclusions more transparent, robust,\nand accessible. Changes (mainly in the Methods section) highlighted in red font\nthroughout the revised version.\n First, we refined the presentation of data and terminology to improve clarity and\naccuracy. In particular, we replaced the term ‘independent’ datasets with\n‘complementary’ indicators, explicitly acknowledging shared upstream inputs\n(AVHRR, MODIS, OCO-2). We also expanded Supplementary Table 1 and added\nSupplementary Table 2 to detail each dataset’s source, resolution, generation\nmethod, and role in the workflow as suggested. In addition, we systematically\nreviewed all abbreviations to ensure they are either spelled out at first mention\nor minimized.\n Second, we strengthened the robustness of our conclusions by conducting a\ncross-indicator consistency analysis, in which pixels (grids) were classified as\nconsistent when the sign of β was supported by at least two-thirds of the\nSOT\ndatasets. Using this criterion, 77.2% of pixels showed consistent directional\nresponses across the six vegetation indicators, despite expected differences in\nmagnitude between greenness- and productivity-based metrics. These results\nindicate that the spatial pattern of sensitivity is robust and not driven by any\nsingle dataset.\n Third, we incorporated a workflow schematic (Supplementary Fig. 4) that links\ndata inputs, sensitivity estimation, mechanistic analyses, and DGVMs\nevaluation. Together with the expanded dataset tables, this figure clarifies the\nrole of each dataset and the overall analytical structure, thereby improving the\ntransparency and reproducibility of our study design.\nWe sincerely hope that these revisions have addressed the main points of concern.\nWe have taken great care to respond thoroughly and transparently, and we believe\nthe manuscript has been substantially improved in both clarity and rigor."}, {"round": 3, "reviewer_comments": {"reviewer_6": "The authors made efforts to address the issues I pointed out before, and the response letter is clear and easy to follow. Some sections still need to clarify. I am making comments based upon the response letter and the highlighted manuscript.\nFirst, it says again all abbreviations were spelled out when first used. This is not correct as I pointed out before. For example, kNDVI, first appeared in line 92, not spelled out (spelled out in line 493 after having been used many times before), followed by SIF, GPP in the same part having not spelled out. This causes difficulties going through this paper. Are Results not main text or does the journal start Introduction followed by Results which are not considered main text? Similarly, TRENDY, ORCHIDEE, YIBs, etc. It is very difficult to figure out these. This needs to be fixed and checked again.\nLine 137-138: “derived from” is duplicated.\nAs in the responding letter, the authors said the RF part was not for predictive modeling purpose, instead it was for demonstrating its spatial heterogeneity and corresponding drivers/factors. Thus, in lines 135-160, words like predictors etc. should not be used to avoid misleading. Also note that permafrost-climate-vegetation is a system, and thus these factors are associated with each other.\nLine 145: Six predictors: permafrost, climate, hydrology, vegetation and soil nutrition, these are five not six. Fig 2 says six key environmental factors: mean SOT, mean temperature, temperature trend, mean SOF, kNDVI trend and mean shortwave radiation. These are confusing. Which is corresponding permafrost factor, climate factor, hydrology factor, vegetation factor and soil nutrition factor? I have not seen any of these six environmental factors that are SOC, ALT and hydrological indicators. This part is confusing, as also reflected in Figure 2 caption. Section 5.2, again, as with line 134 section, confusing about five categories and six environmental factors. Clearly, these six environmental factors do not represent hydrology, soil nutrition, etc. while datasets mentioned ALT, soil and hydrologic datasets. Also noted that hydrologic dataset mentioned multiple parameters and it was not clear which one was used, or all were combined. These two sections need to clarify."}, "author_response": "We sincerely thank the editor for informing us that our manuscript has been\naccepted in principle for publication. We are deeply grateful for the editor’s and\nreviewers’ insightful and constructive feedback throughout the review process.\nIn this final revision, we have carefully addressed all editorial checklist items and\nthe comments from Reviewer #6. Changes highlighted in red font throughout the\nrevised version. Specifically, we:\n(1) standardized all abbreviations and corrected minor textual errors;\n(2) clarified the classification and terminology of environmental factors and\ncategories to remove ambiguity;\n(3) strengthened cross-references and figure citations to ensure consistency\nbetween the main text and Supplementary Information.\nWe believe these clarifications have fully resolved all outstanding concerns.\nDetails of each modification are provided below."}]

Decision letter