Vegetation recovery following retrogressive thaw slumps across northern tundra regions

Recent observations across Arctic and alpine regions reveal substantial greening and browning trends in response to disturbances, such as retrogressive thaw slumps (RTS)^1,2 in ice-rich permafrost regions³. Over recent decades, RTS activity has increased dramatically^4,5,6, with a 60-fold rise observed on Banks Island in the Canadian High Arctic from 1984 to 2015⁷ and a threefold expansion on the Qinghai–Tibet Plateau between 2016 and 2022⁸. RTS occur throughout the Northern Hemisphere permafrost regions^9,10, eroding plant structures (including roots and aboveground biomass), mobilizing previously thawed sediments and soil organic carbon, thereby reducing ecosystem carbon sequestration and increasing microbial decomposition^{11,12,13,14,15,16}. Nonetheless, RTS can stabilize once ground ice is depleted or insulated by sediments^17,18, allowing vegetation succession (gradual process of change in the composition of plant communities over time) to recolonize disturbed areas¹⁹, increasing carbon sequestration. However, the stabilized RTS can re-activate within their disturbed boundaries when a deeper thaw occurs or when ground ice is re-exposed. Consequently, these RTS are considered chronic mass-wasting sites with long-term polycyclic activity (decades to centuries), resulting in complex morphologies and heterogeneous plant compositions^20,21,22. Indeed, the close relationship between plant community composition and the stage of stabilization is well-known and recorded as chronosequences for several RTS field sites^18,23,24,25, yet to date there is limited information for the Northern Hemisphere as a whole about any regional differences of vegetation recovery in terms of timing, composition and succession.

Vegetation greenness can be assessed across large regions using multidecadal satellite-derived normalized difference vegetation index (NDVI) data^3,26,27,28, which serve as an effective proxy for plant productivity. For instance, in Siberian boreal forests, MODIS-derived NDVI indicated greenness recovery (that is, same levels of landscape greenness) 13 years after wildfire disturbance²⁹. Similarly, in Alaskan tundra dominated by lichens, mosses and dwarf shrubs, NDVI values of burned areas returned to undisturbed levels within 10 years, closely corresponding with the recovery time of functional plant diversity³⁰. The growing availability of high-resolution satellite imagery has enhanced NDVI-based assessments, enabling precise detection and monitoring of regional vegetation change trajectories.

Vegetation recovery patterns following RTS disturbances vary widely across tundra regions. For example, in western Siberia, disturbed RTS terrain typically recovers to spectral NDVI levels comparable to undisturbed areas within 20 years³¹. In contrast, on the Qinghai–Tibet Plateau, disturbed areas often remain largely barren even after 20 years^32,33. Recolonization rates in permafrost-affected terrain are influenced by a range of factors, including regional climate (for example, temperature and precipitation)³¹, soil properties (for example, nutrient and moisture availability and soil texture)³⁴, species-specific traits (reproductive effort and dispersal mechanisms)^35,36 and ecological processes (species interactions and succession)³⁷. Despite these insights, vegetation recovery times in RTS remain poorly quantified. The high spatial and temporal variability, driven by dynamic biogeomorphic processes, presents a critical knowledge gap for understanding the resilience of aboveground carbon sequestration in response to widespread permafrost disturbances.

We evaluate the postdisturbance trajectories of vegetation recovery and succession across 12 slump-affected regions in the northern permafrost region. Using time-series image analysis and space-for-time substitution, we incorporate data from repeat satellite, airborne and uncrewed aerial system (UAS) surveys. Vegetation recovery time is defined as the number of years required for NDVI in RTS to return to undisturbed levels, serving as a proxy for recovery in carbon sequestration^28,38. To further investigate the ecological processes underlying spectral greening, we characterize vegetation succession, which represents the change in plant composition over time in response to disturbance. These vegetation change metrics are estimated for eight tundra permafrost ecosystems³⁹ across the Northern Hemisphere, including sites in Alaska, Canada, Siberia and the Qinghai–Tibet Plateau (Fig. 1). The study regions represent diverse plant compositions, nutrient availability and climate conditions, which include (1) Noatak River Watershed (NT) in northwestern Alaska, (2) Brooks Foothills near Toolik Lake (TL) in northern Alaska, (3) Mackenzie Delta (MD) in northwestern Canada, (4) western Banks Island (WB) and (5) eastern Banks Island (EB) in northern Canada, (6) Yamal Peninsula (YP) in western Siberia, (7) central Qinghai–Tibet Plateau (CQ) and (8) western Qinghai–Tibet Plateau (WQ). Detailed site descriptions are provided in Methods. We use PlanetScope imagery⁴⁰ (2016–2024) to construct contemporary greening trajectories based on RTS chronosequences (Extended Data Fig. 1). Piecewise regression⁴¹ identifies breakpoints in NDVI trends, segmenting the time series into distinct stages: undisturbed, initial disturbance, recovery and stabilization. Patterns of vegetation change in RTS are assessed using plant functional types (PFTs) derived from ground surveys, very-high-resolution (<60 cm) satellite, airborne or UAS multispectral imagery from each tundra region, providing insights into associated greenness trends. This study represents the first large-scale assessments of vegetation recovery and succession following RTS disturbance, revealing recovery to vary between years (low-Arctic) to decades (high-Arctic and high-elevation). Notably, the recovery time can be approximated by a model based on a single metric: regional gross primary productivity (GPP) derived from solar-induced chlorophyll fluorescence observations by the Orbiting Carbon Observatory-2 (GOSIF). Model validation using four independent sites, with locations shown in Fig. 1, further confirmed these findings.

Greenness recovery following RTS initiation was rapid in low-Arctic sites but considerably slower in high-Arctic and high-elevation sites, where stable RTS often remained partly barren for decades. In the low-Arctic regions of TL, NT, MD and YP, the NDVI values in undisturbed areas are 0.62 ± 0.03 0.80 ± 0.03, 0.77 ± 0.04 and 0.68 ± 0.03, respectively. Following disturbance, NDVI recovered rapidly to undisturbed levels (0.65, 0.80, 0.75 and 0.69, respectively) within 5–10 years (7.6 ± 0.7 years, 6.1 ± 0.9 years, 5.4 ± 0.9 years and 9.3 ± 0.6 years, respectively), after which vegetation greenness remained stable with only slight increases (Fig. 1). In contrast, recovery in the CQ, WQ, WB and EB regions of the high-Arctic was markedly slower, exhibiting a gradual increasing trend over several decades. These trends indicate that it requires ~49.4 ± 4.4 years, 103.9 ± 90.1 years, 34.6 ± 5.0 years and 105.0 ± 35.8 years for NDVI to recover, respectively (Fig. 1). The greater uncertainty in WQ and EB estimates probably reflects sparse vegetation cover, slower recovery rates and the limited observation period (that is, only 17 years of observations for WQ).

Patterns in vegetation succession varied across tundra regions (Figs. 1 and 2), with rates of change in low-stature vegetation and erect shrubs diverging between low- and high-Arctic and high-elevation sites. In low-Arctic sites (MD, NT, TL and YP), low-stature vegetation rapidly colonized RTS sites within 5–10 years postdisturbance, driving rapid greening (Figs. 1 and 2). Over the following decades, low-stature vegetation was progressively replaced by erect shrubs, the coverage of which is even higher than that of the undisturbed areas (Fig. 2). These estimates align well with local geobotanical surveys in YP⁴². Conversely, vegetation recovery was substantially slower in high-Arctic (WB and EB) and high-elevation (WQ and CQ) sites. Low-stature vegetation was re-established after ~30 years in WB and 50 years in CQ. However, the shorter observation periods of PFT in EB (31 years) and WQ (14 years) showed only minimal increases in vegetation cover, suggesting that full recovery in these regions may take substantially longer.

The variation in NDVI recovery times may be attributed to local and regional factors such as nutrient availability, soil moisture, plant diversity and climate—variables that are broadly represented by regional GPP⁴³. To explore this relationship, we analysed RTS greenness recovery time in relation to GPP data derived from GOSIF⁴⁴, averaged over the period from 2001 to 2024 (Methods). Excluding the WQ region, where GOSIF-estimated GPP is zero, we identified a power-law relationship between recovery time ((tau), years) and GPP (kgC m⁻² yr⁻¹) as τ = 1.35 × (GPP)^−1.68, P < 0.05, R² = 0.97 (Fig. 3). The parameter errors for the scaling factor and the exponent are 0.41 and 0.25, respectively. This relationship suggests that recovery times following thaw-driven disturbance across northern permafrost regions can be readily estimated using regional-scale tundra productivity datasets. Ecosystems with higher productivity exhibited faster recovery. However, notable uncertainty in recovery estimates was observed at extremely low-productivity sites, such as EB, where sparse vegetation cover complicates NDVI trend detection and contributes to increased uncertainty.

Additionally, we assessed the influence of climate by comparing RTS recovery times with regional temperature and precipitation data from ERA5-land, averaged over the same period (2001–2024). No strong relationships were found between recovery time and either temperature or precipitation (Extended Data Fig. 2), suggesting that the climate factors alone do not dominate the vegetation recovery.

To further validate the model, we applied it to four independent RTS sites from different permafrost regions, which included the Peel Plateau (PP), Chukotskiy Rayon (CR), Taymyr Peninsula (TP) and Herschel Island (HI), with mean GPP values of 0.56, 0.36, 0.26 and 0.18 kgC m⁻² yr⁻¹, respectively. Estimated NDVI recovery times based on individual RTS chronosequence for these sites were 4 , 6 , 9 and 25 years, respectively (Extended Data Fig. 3), which correspond well with the model predictions of 3.6 ± 0.5, 7.4 ± 0.5, 13.0 ± 1.0 and 24.1 ± 3.0, respectively (Fig. 3).

Vegetation greenness recovery times in RTS varied widely across permafrost regions. In low-Arctic regions, recovery occurred within 10 years, while recovery required several decades to over a century in high-elevation and high-Arctic regions. Faster recovery in low-Arctic sites such as NT, TL, MD and YP is probably driven by more favourable conditions, including higher nutrient availability, reduced dispersal limitation and greater moisture availability from precipitation and snow melt. Vegetation recovery on the Qinghai–Tibet Plateau and Banks Island was markedly slower. On the forb- and graminoid-dominated Tibetan Plateau^45,46, NDVI recovery spans decades. Recovery in CQ takes ~50 years; whereas in WQ, where initial GPP is near zero, recovery may require up to a century. The harsher climate in WQ, with lower mean annual temperatures (−11 °C in WQ versus −6 °C in CQ) and drier conditions (annual precipitation 339 mm in WQ versus 579 mm in CQ), further limits vegetation regrowth. In the prostrate-shrub tundra of WB and EB⁴⁷, NDVI values reached undisturbed levels after ~30 years and nearly a century. Recovery in EB was substantially slower than in WB, despite similar climate conditions (annual precipitation 233 mm in EB and 227 mm in WB; mean air temperature −13 °C in EB and −12 °C in WB; Supplementary Table 1). Instead, slower EB recovery rates probably reflect drier, highly calcareous gravelly tills that dominate the region⁴⁸. Across these regions, slow recovery thus appears linked to constraints such as low mean annual air temperatures (EB, WB and WQ), limited annual precipitation (EB, WB and WQ), strong solar radiation (CQ and WQ)⁴⁹ and nutrient-poor soils (CQ, WQ and WB)⁵⁰.

During recovery, NDVI values within RTS scars returned to undisturbed levels as early colonizers (for example, forbs, graminoids and shrubs) began to dominate the terrain (Fig. 4a). Rapid NDVI recovery in low-Arctic sites was probably driven by the early colonization of low-stature vegetation, including mosses, graminoids and short dwarf shrubs (Fig. 4b–d,f). Over the following 10–50 years, erect shrubs gradually established and stabilized the terrain, potentially acting as catalysts for further vegetation change⁵¹ and enhancing carbon sequestration⁵², thereby contributing to a negative feedback that enhances carbon uptake potential. Although NDVI values may be at a similar level, the vegetation composition in these disturbed areas differs. Spectral recovery is therefore not synonymous with ecological recovery, as observed in other disturbance works⁵³. Comparisons between disturbed and undisturbed terrain reveal a clear state change in vegetation composition, such as tussock-shrub tundra transitioning to shrub tundra (Fig. 4b–d,f). Under conditions of increased nutrient and moisture availability, shrubs not only expand in coverage but also exhibit trait changes, including increased height relative to undisturbed areas (Fig. 4b,c). The classification of ‘undisturbed’ terrains in this study is based on the best-available knowledge and data. However, it is important to note that some areas categorized as undisturbed may have experienced past disturbances that are no longer detectable in modern imagery or records.

Recovery and succession of vegetation may be influenced by nutrient availability, initial plant communities and climatic conditions. Nutrient limitation is a key constraint on vegetation growth in permafrost regions^54,55, particularly on the Qinghai–Tibet Plateau, where organic layers are shallow and leaching is common. Even after several decades, graminoids and forbs continue to dominate the recovered scars in CQ (Figs. 2 and 4e), probably due to the lack of nearby shrub seed sources and harsh environmental conditions that inhibit the establishment of more productive species. In high-Arctic WB and EB, low-stature vegetation, including graminoids, forbs and dwarf shrubs⁴⁷, has dominated disturbed areas for decades, with minimal shifts in community composition, also probably due to the harsh climate that limits plant growth. In contrast, more southerly ecosystems in northern Alaska, northwest Canada and western Siberia benefit from nutrient-rich soil, suitable summer temperature and higher soil water content due to precipitation and snow melt, which supports seed germination, nutrient availability and root development, ultimately accelerating vegetation regrowth.

The rates of vegetation recovery following RTS disturbances in tundra regions occur more slowly than in drained thermokarst lake basins but are generally faster than recovery following wildfires across the low-Arctic regions. For example, postfire vegetation in boreal Siberia forests typically recovers within ~13 years²⁹, while Alaskan tundra dominated by dwarf shrubs recovers in ~10 years³⁰. In contrast, vegetation in drained thermokarst lake basins can re-establish within just 2 years⁵⁶. Postfire successional patterns in Alaskan tundra also indicate rapid resprouting of vascular plants and Sphagnum moss, followed by a shift towards shrub-dominated communities as nutrient availability increases, resulting in greenness levels that can exceed those of undisturbed areas³⁰, which is similar to that post-RTS disturbance.

The variability in vegetation recovery times across tundra ecosystems can be effectively estimated using regional-scale GPP through a power-law relationship. Validated by four independent regions, we found that ecosystems with higher GPP generally exhibited faster recovery following RTS disturbances. This predictive capability probably stems from the strong associations between high GPP and factors such as suitable climate, greater plant species diversity^57,58, adequate nutrient availability and high soil moisture—all of which enhance vegetation growth and ecosystem resilience³⁸. Diverse ecosystems often host a broader array of plant species with varied strategies for growth, survival, seed dispersal and reproduction, as well as expanded plant reproduction and dispersal⁵⁹, collectively promoting faster recovery and the restoration of ecological functions. For instance, ecosystems with a higher diversity of graminoid and shrub species exhibited notably shorter recovery times. Importantly, although we successfully predicted the recovery of surface greenness over several years to decades, plant species diversity did not fully return within the time frame of this study¹⁹ (Fig. 4b,d,f).

Although RTS are spatially limited—occurring only in areas with ice-rich permafrost, gentle slopes and specific soil characteristics, the ongoing warming of high-latitude and high-elevation regions is making areas with ice-rich permafrost increasingly vulnerable to these disturbances. RTS accelerate soil carbon mobilization and vegetation removal, but in high-productivity tundra regions, disturbed microsites can foster warmer, nutrient-rich and moist conditions that support rapid vegetation regrowth⁵⁵. In such environments, surface greenness often recovers within a decade and may later exceed that of the surrounding undisturbed areas, contributing positively to carbon sequestration. However, in high-Arctic and high-elevation regions, unfavourable environmental conditions leave ground surfaces exposed for decades. As RTS continue to expand across permafrost landscapes^4,6,7,8,60, this study provides critical insights into the spatiotemporal variability in vegetation succession, surface greening and browning, which reflect carbon sequestration following permafrost disturbances across northern tundra regions.

Study areas and datasets

We selected eight study areas (Supplementary Table 1) located in various tundra ecosystems, including NT and TL, MD, WB and EB, YP and CQ and WQ. The ecosystem of NT and TL is characterized by Arctic tussock-shrub and low-shrub tundra. The MD study area is located in northwestern Canada, on the tundra uplands east of the MD that are characterized by rolling hills and small lakes⁶² and where tall shrub tundra and erect dwarf shrub tundra dominate. WB and EB are located in the high-Arctic tundra, where sedge, dwarf shrubs, prostrate dwarf shrubs and moss tundra are dominant⁴⁷. YP is located within an Arctic tundra, a dissected lowland with long slopes composed of sandy clay deposits in western Siberia^63,64, with erect dwarf shrubs, dwarf shrubs, sedge and moss tundra distributed⁴⁷. Meanwhile, the CQ and WQ host high-elevation alpine meadows and grasslands with specialized species that are resilient to extreme high elevations. Four additional RTS sites in CR, TP, HI and PP served as validation locations. The CR and TP are located in eastern and northern Siberia, characterized by sedge and dwarf shrubs⁴⁷. HI, located in the northernmost part of Canada, is covered by fine-grained marine sediments⁶⁵ and characterized by tussock sedge, dwarf shrub and moss tundra⁴⁷. The PP in northwestern Canada is dominated by woodlands, shrubs and dwarf shrubs⁶⁶.

The core datasets (Supplementary Table 2) used to delineate the RTS include PlanetScope images from 2016 to 2024, Landsat-5 and Landsat-7 images (1984–2010)⁶⁷. High-resolution images include Geoeye-1 (2009–2024), QuickBird-2 (2002–2014) and WorldView-1–3 (2007–2024), Keyhole-4 and Keyhole-9 images since 1965^68,69, Aerial Photo Single Frames (1937–2014)⁷⁰, Aerial Images (1952, 1970), airborne imagery in 2023 and UAS images since 2017. The PlanetScope images from 2016 to 2024 were applied to calculate the NDVI maps. We also applied Harmonized Landsat and Sentinel-2 (HLS) products, which provide the bidirectional reflectance distribution function (BRDF) adjusted Landsat and Sentinel-2 images, to calibrate the PlanetScope images. The optical images from UAS and airborne (<0.1 m), as well as high-resolution satellite imagery (around 0.5 m) since 2016, were used to calculate the percentage cover of PFTs. Then, we averaged the Orbiting Carbon Observatory-2 based, solar-induced chlorophyll fluorescence (GOSIF) derived GPP⁴⁴ in 0.05° × 0.05° resolution since 2001 to measure the growth of the terrestrial vegetation. The monthly mean air temperature and precipitation were collected from the ERA5-land reanalysis data⁷¹.

To further analyse the site distributions, we used topographic, vegetation, permafrost distribution and soil datasets, including (1) vegetation type (data source—the circumpolar Arctic vegetation map, description of western Canada landscapes and landforms⁶⁶ and Wang et al.⁷²), (2) digital elevation model from the shuttle radar topography mission⁷³, (3) soil textures⁷⁴ and (4) permafrost map⁷⁵.

Generating NDVI time series

We used space-for-time substitution (thaw slump chronosequence)⁷⁶ to examine changes in greenness before, during and after RTS disturbance (Extended Data Fig. 1). The chronosequence includes RTS initiated in different years and undisturbed land. First, for each regional ecosystem, we selected at least 15 representative RTS (Supplementary Table 1), identified as either active from headwall expansion in recent years or stabilized for decades. We modified and delineated the disturbed area for every RTS based on the published inventories^7,8,77,78 and the high-resolution images. Then, we divided it into subparts according to the respective initiation years. Reactivated RTS were treated as new perturbation phases in our analyses, as debris deposition resets portions of the thaw slump scar by removing or covering vegetation and can alter soil moisture and nutrient status. We identified several undisturbed areas near the RTS to serve as controls for modelling undisturbed conditions. The NDVI time series of RTS, obtained from PlanetScope images (2016–2024), were used to evaluate the vegetation recovery times. To generate NDVI maps (Extended Data Fig. 4), we first downloaded the PlanetScope images using bounding boxes (at least 2 × 2 km²) covering RTS-disturbed areas to provide sufficient valid pixels after quality control. The orthorectified PlanetScope images, including red, green, blue and near-infrared bands, as well as the quality-control layer, can be obtained through https://www.planet.com/. Unusable data mask layers were applied to the PlanetScope images to remove low-quality pixels. As the PlanetScope images are subject to the BRDF effect, which results from varying illumination and sensor viewing angles during image acquisition, we cross-calibrated PlanetScope surface reflectance based on the BRDF-adjusted HLS products. We filtered the HLS based on its band-specific quality control (QA layer) to retain high-quality pixels that are not affected by heavy cloud cover. Then, we selected one HLS image that contained the most valid pixels during every peak-growing season. The peak-growing season was determined for each subregion using long-term (2000–2024) local phenological patterns of NDVI from Landsat-7 and Landsat-8 (ref. ⁷⁹). For each year, we selected one PlanetScope image acquired within the nearest 10 days relative to the HLS image.

To generate the consistent NDVI time series, we first coregistered the yearly PlanetScope images based on HLS images using AROSICS⁸⁰. Then, we cross-calibrated PlanetScope surface reflectance based on the corresponding HLS images using histogram matching⁸¹. The method first upscales the PlanetScope images to the HLS resolution and calculates the band-specific histograms. A linear transformation was then applied to downsampled PlanetScope images to match the mean and standard deviation with those of the HLS data. The recorded band-specific adjustment coefficients were applied to the original resolution of PlanetScope images. Finally, we calculated the NDVI as (near-infrared − red)/(near-infrared + red).

We then averaged the NDVI values on the basis of the RTS delineation-derived subparts and created NDVI time series according to the initiation years of these subparts. To understand how vegetation changes since initiation, we transferred the timescale from the calendar year of the images to the number of years that have occurred since RTS initiation (Extended Data Fig. 1). Assuming RTS in the surrounding areas with similar climate conditions and vegetation cover would exhibit similar vegetation recovery trajectories, we generated the NDVI time series using all the selected RTS within each regional ecosystem.

Estimating NDVI recovery time

The pronounced difference in vegetation cover before and after RTS initiation separates the NDVI time series into several segments. To best fit the NDVI trend, we applied a piecewise regression method⁴¹ to separate the NDVI time series into distinct periods: undisturbed, disturbance, recovery and stabilization (where applicable, as not all regional ecosystems exhibited a stabilization period). For each period, we applied linear regression to estimate the trend. Finally, we estimated the recovery time of NDVI as the number of years required for NDVI values to return to predisturbance levels or become stable. For the ecosystems with a stabilization period, the recovery time was determined as the temporal separation between the breaking points of the recovery and stabilization periods. For ecosystems without a stable trend in the recorded NDVI timescales, we hypothesized that the NDVI value would continue increasing at the same rate until reaching undisturbed levels, then turn stable. The recovery time was estimated from the difference between the average NDVI value before the RTS disturbance and the NDVI value at the time of disturbance, divided by the subsequent recovery trend. The uncertainty of the recovery time was estimated by accounting for the trend uncertainty.

Linking PFT cover fraction with NDVI

We generated PFT maps in the same year as PlanetScope images using UAS, Geoeye-1, QuickBird-2 and WorldView-1–3 images. To enhance the spatial resolution of the composited red–green–blue images for visual interpretation, we applied pan-sharpening to the multispectral bands using the panchromatic band for GeoEye-1, QuickBird-2 and WorldView-2 and -3. To match the resolution of the PlanetScope images, we divided the high-resolution images into 3 × 3 m² grids. We classified the vegetation cover into three PFTs: bare ground, low-stature vegetation and erect shrub and manually assigned the type to every grid based on the textures and colours observed in the images (Extended Data Fig. 5). Although the resolution of satellite-based images is relatively low (~50 cm), the greener colour and coarser texture of erect shrubs presented in the images help us differentiate them from low-stature vegetation. The percentage of PFT cover was calculated for each RTS subpart. Using these percentage cover values, we generated composite colours (red for bare ground, blue for erect shrubs and yellow for low-stature vegetation) and assigned these colours to the corresponding NDVI time series (for example, Fig. 1).

Fitting relationship between GPP and recovery time

We used orthogonal distance regression to fit a power-law model between the NDVI recovery time and the mean GPP across the regional ecosystem while accounting for the errors in both variables. To mitigate interannual variability of GPP, which was observed in both tower-measured and estimated data⁴⁴, we averaged GOSIF GPP over 2001–2024 and calculated annual standard deviations. Regional mean GPP was then computed within 2.5-km buffer zones around RTS sites, with pixel-level errors propagated throughout the analysis. In addition, when fitting the power-law model, we added the associated error bars from the site-level variability in recovery time (vertical error bars) and the regional variation of GPP (horizontal error bars). On the basis of the covariance matrix of the best-fit model, we calculated the confidence interval. We also set lower (0.05 kgC m⁻² yr⁻¹) and upper (0.6 kgC m⁻² yr⁻¹) bounds of the function based on the distribution of GPP in permafrost areas covered with graminoid and shrubs, respectively. To validate the fitted model, we chose four independent RTS sites in CR, PP, HI and TP (Fig. 1). The selected RTS (Extended Data Fig. 3) with subparts initiated at different times, with the colours in the high-resolution images indicating when the greenness recovers. We compared the recovery time with the model-inferred ones to validate the performance of the model.