Canopy-mediated climate feedbacks in the boreal continuous permafrost zone

Approximately a quarter of the continuous permafrost zone is forested (Fig. 1), storing substantial carbon, 95% of which is below ground in the active layer and permafrost^1,2,3. These pools are increasingly vulnerable as the boreal biome warms nearly four times faster than the global average, with mean annual air temperatures projected to rise by 6–11 °C by 2100 under high-emissions scenarios⁴. Although some regions show increased carbon uptake through greening, these gains are heterogeneous and potentially offset by wildfire and wetland emissions^5,6,7. Beyond the net ecosystem carbon balance, forest and active layer properties are important factors in preserving near-surface permafrost and the soil organic carbon (SOC) stored therein in a state of thermal disequilibrium with climate conditions. This makes the system highly sensitive to disturbances, underscoring the need to quantify vulnerabilities of the predominantly soil-based carbon stock^8,9.

SOC stocks in boreal permafrost regions vastly exceed live biomass, storing approximately 1,000 ± 200 Pg SOC in the top 3 m (refs. ^10,11), representing 33% of global SOC despite covering only 15% of global soil area¹² (Fig. 2). The vulnerability of this SOC upon thaw is not uniform but is determined by prefreezing decomposition, carbon origins and age, influenced by formation mechanisms, burial rates and environmental conditions^13,14, as well as topographic, edaphic and depositional processes controlling contemporary carbon distribution^15,16. The thermal protection of this heterogeneous reservoir is critically dependent on the thermal buffering capacity of the forest canopy, which maintains the frozen state of permafrost through its control over surface energy fluxes.

The canopy regulates ground thermal regimes through multiple competing mechanisms¹⁷. Summer shading reduces incoming solar radiation by up to 90%, cooling the ground and limiting thaw depth¹⁸. However, the magnitude of this is largely controlled by canopy structure and canopies also exert warming effects, such as lower surface albedo, and longwave radiation retention by dense canopies¹⁹. The net effect is a dampened seasonal temperature amplitude. The reduced summer warming drives shallower thaw depths than in open landscapes, while the relative dominance of summer versus winter warming effects depends on climate, canopy structure and snow cover regime^20,21,22. This thermal buffering is reinforced by synergistic mechanisms including reduced soil thermal conductivity from dry soils, modulated snowmelt timing, and insulating litter and organic layer accumulation^9,23,24.

Forest loss triggers active layer deepening and can accelerate thermokarst formation in ice-rich permafrost²⁵. However, Earth system models poorly represent these forest–ground thermal interactions^26,27, particularly in permafrost regions^{28,29,30,31,32}. The stability of the permafrost carbon pool, thermally regulated by the canopy, is jeopardized by forest loss from compound disturbances. These include logging, wildfires, drought and storms^33,34,35, often amplified by climate-intensified biotic threats such as insect outbreaks^36,37,38. This confluence underscores the need to quantify the size of the carbon pool currently kept in a frozen state by canopy insulation, as it is increasingly vulnerable to destabilization.

Boreal forest disturbances exert complex and opposing thermal effects on climate, complicating forest-based climate interventions^39,40,41: while forest loss causes regional atmospheric cooling through increased snow-season albedo^{42,43,44,45,46,47,48}, it simultaneously diminishes ground insulation, increasing heat flux to subsurface layers and accelerating permafrost degradation^22,49,50,51. These atmospheric and subsurface thermal responses are decoupled: increased winter albedo reduces net radiation at the surface, cooling the atmosphere and ground, while in summer, increased shortwave radiation at the ground surface and altered turbulent heat fluxes drive subsurface warming and active layer deepening. Consequently, the perceived climatic benefit of albedo-induced atmospheric cooling may be offset locally by heightened subsurface warming and soil carbon loss risks. While boreal afforestation can alter permafrost stability through radiative and snow-mediated effects⁴, forest removal increases ALT despite subsequent regrowth, due to persistent changes in shading, heat fluxes and shallow rooting depths that limit evapotranspiration cooling. However, the magnitude of forest thermal buffering and consequently the vulnerability of permafrost-protected SOC remains unquantified across the climatic and biogeographic gradients of the continuous permafrost zone^4,52.

We present a process-based modelling framework coupled with a comprehensive meta-analysis to assess how boreal forest canopies regulate soil carbon storage by controlling permafrost thermal dynamics across climatic and biogeographic gradients. We quantify the thermal buffering capacity of forests by comparing seasonal thaw depths in forested versus unforested scenarios using paired simulations within seven climate zones of the forested continuous permafrost zone⁵³ (Extended Data Fig. 1 and Extended Data Table 1). Our analysis includes three dominant plant functional types (PFTs) and multiple canopy density classes to identify regions where forest cover most effectively insulates SOC from seasonal thaw (Extended Data Fig. 2 and Extended Data Table 2). We complement this with a first assessment of the vulnerabilities of different SOC stocks to forest loss across biogeographic gradients in forested continuous permafrost regions.

Soil C stocks exceed biomass C in forested permafrost

Biomass carbon in the forested continuous permafrost zone is approximately 7.2–18.6 Pg (Fig. 3c). While aboveground live biomass represents the largest pool with 4.5 Pg (ref. ⁵⁴), the other components (root^55,56,57,58, below-canopy^59,60 and litter biomass^61,62,63,64) show substantial regional variability between Eurasian and North American boreal forests (Supplementary Section 6). Eurasian forests, adapted to nutrient-poor soils and shallow active layers, allocate less carbon to roots and accumulate more litter, while North American ENF stands invest more heavily in belowground biomass⁶⁵. The Northern Hemisphere (NH) permafrost region stores over 1,000 Pg of SOC in the top 3 m, with 188 Pg located in forested continuous permafrost^11,12,66. Of this, 64.5% is currently in a frozen state (118 Pg)⁶⁷, clearly exceeding biomass stocks, which amount to only 3–8% of the SOC reservoir.

Canopy thermal regulation sustains shallow active layers

Simulations with the CryoGrid permafrost model (Extended Data Fig. 3) find that a forest canopy substantially lowers active layer thickness (ALT) (Supplementary Fig. 2). Comparing forested and non-forested scenarios under current climatic conditions, the simulated 2018–2023 average ALT is 50.2–62.3% larger under bare-ground scenarios. In absolute terms, the canopy reduces the ALT by 0.87–1.43 m compared with bare-ground conditions (interquartile range 0.98–1.18 m). We find regional variation in the insulating effect. In North America, ALT is reduced by 1.0–1.4 m (mean ± s.d.: 1.2 ± 0.13 m), indicating a stronger insulating effect compared with Eurasia (1.02 ± 0.1 m).

The strength of the canopy’s insulating effect varies by PFT and leaf area index (LAI) (Supplementary Fig. 2). Deciduous needleleaf forest (DNF), covering close to three-quarters of the forested continuous permafrost area (2 million km²) and composed predominantly of Larix spp.⁶⁸, exhibits a modest, regionally consistent effect (1.03 ± 0.07 m; North America: 1.08 ± 0.04 m; Eurasia: 1.02 ± 0.08 m). Evergreen needleleaf forest (ENF), covering one-tenth of the area (0.3 million km²) primarily as stands of spruce, pine and fir in northern Canada, whose persistent canopies intercept snowfall and retain longwave radiation⁶⁹, show a moderate effect (1.10 ± 0.15 m), with strong regional differences: ENF in North America displays substantial ALT reduction (1.18 ± 0.1 m) while the effect is weaker in Eurasia (0.93 ± 0.05 m). Broadleaf deciduous forest (DBF), occupying one-fifth of the zone (0.4 million km²) with dominant species including birch, aspen and poplar⁷⁰, exhibits the largest reduction (1.21 ± 0.18 m), predominantly driven by stands along the North American treeline (1.4 ± 0.03 m), with Eurasia showing moderate reductions (1.11 ± 0.15 m). Sparse canopies (LAI = 1 m² m⁻²) provide the strongest ALT reduction, especially in North America (1.3 ± 0.11 m) but also in Eurasia (1.1 ± 0.11 m). Denser canopies produce a smaller ALT reduction compared with sparse canopies in both North America (1.1 ± 0.13 m) and Eurasia (1.0 ± 0.08 m).

The nonlinear relationship between canopy density and ALT reduction (Supplementary Figs. 3 and 4a) reflects competing radiative and hydrological controls (Supplementary Figs. 4 and 5). Sparse canopies substantially reduce summer shortwave radiation and ground heat flux while maintaining relatively high soil moisture, which sustains thermal conductivity and enhances heat extraction during freeze-back. By contrast, dense canopies, despite stronger shading, promote longwave radiative trapping and higher evapotranspiration rates that dry the soil, lowering thermal conductivity and limiting winter heat loss. This moisture-mediated feedback partly offsets the increased shading, resulting in slightly deeper active layers under high compared with intermediate canopy densities (Supplementary Fig. 4a).

Simulation results for 2023 for the Eurasian deciduous needleleaf-dominated cluster in the continental subarctic (Dfc) climate zone across the four canopy density scenarios (bare, sparse, mid and dense) support this. During the snow-free season (June–September), mean incoming shortwave radiation at the ground decreased by 91% from bare (194.5 W m⁻²) to dense (16.7 W m⁻²), and by 62% from sparse (43.7 W m⁻²) to dense, reflecting progressive canopy interception. By contrast, mean outgoing longwave radiation during the snow season (October–May) increased only modestly (3.9% from bare to dense; 3.6% from sparse to dense). Mean summer canopy transpiration increased by 66% from sparse to dense (0.81−1.34 mm per day), consistent with enhanced soil drying at higher LAI (Supplementary Figs. 4c,d and 5).

Forest canopy attenuates gradual and rapid SOC thaw

SOC in permafrost regions is vertically stratified but retains substantial carbon at depth due to cryoturbation and frozen paleosoils. Because the majority of the SOC stock is in the upper 0–3 m, maintaining a shallow ALT via canopy insulation notably reduces the amount of SOC exposed to seasonal thaw (Fig. 3b). Depending on the LAI, canopy presence limits thaw exposure by 37.9–41.1 Pg, equivalent to 3.2% of the entire NH permafrost SOC pool and 33.5% of the frozen SOC in forested continuous permafrost (Fig. 4b and Supplementary Table 3).

The largest amount of frozen SOC preserved by canopy buffering is found in Eurasian DNF stands within the continental subarctic (Dfc) climate zone, with storage densities from 16.7 to 18.3 kg m⁻² (Fig. 4b). However, the highest SOC densities are located in North American DBF stands (in Dfc), with 42.9–44.5 kg m⁻² (Fig. 4a). ENF in subarctic climate with dry summer (Dsc) and Dfc zones of North America preserve 22.2–25.2 kg m⁻² and 36.5–40.9 kg m⁻², respectively. Overall, SOC densities span a wider range in North America (12.7–44.5 kg m⁻²) than in Eurasia (15.3–20.7 kg m⁻²; Fig. 4b).

To account for abrupt thaw in high excess-ice areas (Supplementary Fig. 6), we quantify the thaw exposure of the full 0–3 m SOC profile. In ice-rich permafrost, total degradation combines active layer deepening and ground subsidence from excess ice melt⁷¹. Thermokarst processes mobilize carbon 10–100 times faster than gradual top-down thaw alone⁷², with consequences varying across thermokarst modes depending on slope position, soil texture and ground ice content. Thermokarst extent across interior Alaska has nearly doubled since 1949, with increasingly transformative expansion projected through 2100⁷³. In addition, a forest canopy delays thermokarst onset by 3–18 years and slows excess ice melt by up to 7 years compared with bare-ground simulations, and post-disturbance subsidence progresses within 2–35 years depending on cryostratigraphy and climate²⁵, with fire further capable of triggering complete permafrost degradation within an 18-m soil column within 120 years in upland forests⁷⁴. The potential for rapid and deep thaw highlights that an additional SOC pool of 17.7–20.9 Pg is currently protected by canopy thermal insulation in high excess-ice areas (Fig. 3d).

Our analysis reveals that forest canopies act as a critical thermal buffer, preserving 59 Pg of SOC in a frozen state. This pool represents 30% of the total forested permafrost SOC stock, far exceeding the 7–19 Pg stored in boreal biomass pools, which underscores the canopy’s disproportionate role in stabilization. However, this protection is fragile, and canopy loss through disturbance could mobilize ancient carbon stocks, with the climate impact hinging on the SOC’s biogeochemical vulnerability and the processes driving forest loss (Fig. 5 and Supplementary Table 2).

SOC vulnerability across boreal permafrost landscapes reflects substantial heterogeneity in soil genesis and depositional history. These differences in carbon origin, formation history and edaphic factors give rise to pronounced contrasts in thaw sensitivity^13,14,15,16. One of the most labile pools is found in the Yedoma deposits, accumulated during the cold, arid late Pleistocene and under limited decomposition (Supplementary Fig. 6). Yedoma is characterized by low carbon-to-nitrogen (C:N) ratios, preserved labile compounds and high ice content (>80%)⁷⁵. Yedoma deposits are partly forested, with the canopy protecting 1.7 PgC. This ancient C is highly vulnerable to rapid decomposition upon thaw, representing a pathway for long-term SOC losses that are effectively irreversible on policy-relevant timescales, complicating temperature overshoot scenarios where temporary exceedance of warming thresholds could trigger permanent carbon release^14,72,76,77.

However, the majority (39 Pg) of thermally protected SOC in forested continuous permafrost resides in other soil classes with distinct formation histories, carbon stabilization mechanisms and thaw vulnerabilities (Extended Data Fig. 4 and Supplementary Table 3). Turbels, formed through cryoturbation processes that physically mix organic matter into mineral horizons, store 14.7 Pg across 32.6% of the forested continuous permafrost area and exhibit high vulnerability to thaw. Orthels, the most extensive class (48.6%), contain the largest carbon pool at 17.2 Pg with medium vulnerability^78,79. Organic-rich Histels and Histosols, representing peatland systems with generally lower vulnerability due to prefreezing humification, although labile fractions remain^80,81,82, store 6.7 Pg combined across 9.4% of the area. Soil class composition differs across PFTs (Supplementary Fig. 7). DNF forests dominate the continuous permafrost zone and are primarily underlain by Orthels (59.9%) and Turbels (28.0%), with the highest Yedoma occurrence (6.0%) of the three PFTs, thereby concentrating medium- to high-vulnerability soils. ENF and DBF forests are instead dominated by Turbels and Histels, reflecting wetter, more cryoturbated conditions and low Yedoma occurrence. Understanding the distinct characteristics and vulnerabilities of these soil classes is critical for projecting carbon dynamics under forest disturbance and climate change.

Empirical evidence confirms that vulnerable C can be rapidly mineralized. Sampling across 18 farms spanning permafrost and non-permafrost soils revealed that conversion of permafrost-affected forests to cropland resulted in SOC losses of 15.6% ± 21.3%, while conversion to grassland resulted in losses of 23.0% ± 13.0% (ref. ⁸³). Notably, the deforestation led to permafrost disappearance, indicating that the land-use change accelerated warming and subsequent thawing. Conversely, soils without permafrost showed no substantial SOC loss and, in some cases, slight SOC gains under grassland^83,84, demonstrating that the thermally protected SOC is not only theoretically vulnerable but is already being lost (Supplementary Section 4.1).

The shallower active layer relative to open ground reflects the net balance of competing radiative and hydrological canopy effects (Supplementary Figs. 4 and 5). The thermal buffer provided by the canopy is increasingly compromised by intensifying disturbance regimes, which are amplified by climate warming, drying and earlier snowmelt^34,85,86. Wildfire represents the most immediate and widespread disturbance. Beyond combustion of biomass and near-surface SOC, directly transferring carbon and nitrogen to the atmosphere⁸⁷, fire can disturb the insulating vegetation layer and lead to accelerated active-layer deepening and soil warming²². The consequence is a reduction in carbon stocks. In burned and subsided soils, carbon pools in the upper 1.25 m are reduced by more than half (from 7.3 kg m⁻² to 3.1 kg m⁻²), with only approximately 10% lost during combustion and the remainder through subsequent thaw-driven microbial decomposition^88,89,90 (Supplementary Section 4.2).

Disturbance regimes can exhibit strong positive feedbacks that promote further disturbances. Wind damage and insect outbreaks generate deadwood that fuels subsequent fires and pathogen attacks, with cascading interactions triggering compounding impacts on soil thermal stability and carbon release^36,91,92. Changes in species composition, such as shifts towards broadleaf deciduous taxa that are less fire-prone, may partially counteract some feedbacks, although resilience is limited under rapidly changing disturbance regimes^70,93. In ice-rich Yedoma regions, any such disturbance can destabilize ice wedges, triggering thermokarst collapse within months to years^25,68,94.

Our analysis reveals a critical paradigm: for the global carbon budget, the primary significance of boreal forests in continuous permafrost zones lies less in their biomass and carbon sequestration capacity, but in their capacity as a thermal insulator for the substantially larger soil carbon pool. The live biomass and litter pools store an estimated 7.2–18.6 Pg (Supplementary Section 6). By contrast, the canopy’s cooling effect maintains approximately 59 Pg SOC frozen, a pool three times larger than the total biomass-related carbon. Spatially, our simulations show that boreal canopies most effectively insulate permafrost and stabilize frozen SOC in zones with high carbon storage (Fig. 3b), with biomass-related carbon playing a smaller overall role (Fig. 3c). Particularly in northeastern Canada, canopy loss would disproportionately impact permafrost-stabilized SOC stocks, as these areas combine high carbon densities with strong canopy-mediated thermal buffering. In addition, disturbances also alter nitrogen cycling, as shown by substantial losses in total and labile nitrogen, along with depletion in foliar nitrogen⁹⁰.

The nonlinear dependence of thermal protection on canopy density (Supplementary Figs. 2, 4 and 5) means that partial canopy loss from drought stress, insect outbreaks or selective logging may not immediately compromise permafrost stability in proportion to canopy loss^20,22,95, and the continuous presence of a canopy is more critical for permafrost stability than the specific forest composition or canopy density.

Enhanced vegetation growth from CO₂ fertilization, warming or increased nitrogen availability could partially offset SOC losses^96,97 but increased plant productivity does not directly translate into higher ecosystem carbon storage^{98,99,100,101,102,103,104} (Supplementary Section 4.3). Birch forest expansion into tundra may stimulate SOC decomposition through priming effects, potentially reducing net carbon storage despite higher productivity¹⁰⁵, while initial carbon gains from boreal forest densification¹⁰⁶ eventually plateau due to growth limitations and climate-induced species mortality¹⁰⁷. The SOC pool will probably release more carbon through thaw and decomposition than can be sequestered by vegetation expansion, with losses persisting even under net-zero emissions scenarios due to ongoing soil respiration¹⁰⁸.

Our study reveals that the potential costs of boreal forest loss extend far beyond the immediate reduction in tree respiration and biomass carbon storage (Supplementary Section 4.4). Our analysis highlights that the forest-affected continuous permafrost region contains approximately 183 Pg of carbon in the top 3 m of soil, which represents about 15% of the total NH permafrost SOC pool. Currently, 32% of this SOC stock is stabilized by the thermal buffering effect of the forest canopy, while the canopy’s own carbon represent a comparatively smaller contribution.

These large SOC pools demonstrate that boreal forests provide critical ecosystem services on regional to global scales. Our findings have direct implications for forest management and climate policy. The canopy-mediated preservation of 59 Pg frozen SOC represents avoided emissions rather than active sequestration, a distinction with consequences for how boreal forests are valued in carbon accounting frameworks. Forest management in continuous permafrost regions must prioritize maintaining canopy cover for thermal protection, even through sparse stands with increased drought resistance and reduced fire-proneness^109,110. While fire is a natural disturbance in boreal ecosystems¹¹¹, targeted management in high-carbon, ice-rich permafrost areas, where canopy loss could trigger irreversible carbon mobilization, may offer climate benefits. A pilot project by the US Fish and Wildlife Service in the Yukon Flats demonstrates such prioritization¹¹², offering a framework compatible with Indigenous-led stewardship¹¹³. Current policy frameworks inadequately account for permafrost thaw and SOC vulnerabilities, revealing a critical gap that our findings can help address.

Conclusion

In boreal landscapes, the net ecosystem C balance is driven primarily by SOC dynamics, not aboveground biomass¹¹⁴. We found a strong linkage between SOC stabilization and forest canopy cover. Consequently, forest loss could disproportionately destabilize the vast boreal continuous permafrost SOC pool, impacting regional and global carbon cycles. The 58.8 PgC (4.8% of the total permafrost SOC pool) thermally protected from seasonal thaw by forests, underscores the critical importance of maintaining canopy cover to mitigate large-scale emissions from thawing boreal permafrost soils.

Our simulations demonstrate that forest canopies critically modulate ground thermal and hydrological regimes, with the vulnerability of permafrost to canopy loss varying regionally based on climate and PFT. Rather than aiming to project future permafrost trajectories, we quantified the present-day thermal and hydrological buffering effects of existing forest canopies on permafrost stability. Our findings underscore the necessity of incorporating detailed canopy–permafrost interactions into larger-scale permafrost models, as the insulating effect of boreal forests is not uniform but contingent on forest density, composition and regional climate conditions. The observed spatial heterogeneity in permafrost responses further indicates that accurate vulnerability assessments must account for these ecosystem-specific controls.

Complete boreal forest removal in the continuous permafrost region would expose 40–59 Pg of permafrost SOC to thaw, equivalent to 3.2–4.8% of permafrost SOC of the NH and 4.3–6.4% of its frozen carbon stock. This SOC pool at risk dwarfs potential carbon losses from aboveground biomass (4.5 Pg) and other biomass-related stocks (subcanopy pool, roots and litter; 7.2–18.6 Pg).

The disproportionate impact of canopy removal on boreal SOC stability would result in a substantially negative net C balance in boreal ecosystems, revealing a key risk for future greenhouse gas emissions from boreal permafrost. This protective role of the canopy far exceeds its contribution to standing biomass carbon stocks, underscoring that canopy conservation is critical for mitigating the permafrost carbon feedback. Our findings can guide the planning management policies for protecting large natural carbon sinks as well as preventing additional carbon losses.

We investigated forest canopy effects on continuous permafrost by integrating satellite-derived vegetation products, climate reanalysis data and numerical modelling. A supervised classification approach grouped the domain into 16 bioclimatic clusters based on:

1. (1)

   climate zones: seven subarctic and semi-arid Köppen–Geiger climate classes defined by temperature, seasonality and precipitation⁵³;

2. (2)

   PFT: Moderate Resolution Imaging Spectroradiometer (MODIS)-derived boreal forest PFTs (ENF/DNF, DBF)¹¹⁵;

3. (3)

   continuous permafrost extent⁶⁶;

4. (4)

   clusters established separately per continent to account for regional climate–vegetation interactions and climate history.

The classification framework was implemented through three sequential steps. First, climate zones were delineated using the Köppen-Geiger system⁵³, which categorizes regions on the basis of critical temperature thresholds and seasonal patterns, emphasizing the coldest month of the year as a determining factor (Extended Data Fig. 1 and Extended Data Table 1).

Second, forest PFT were extracted from MODIS MCD12Q1 at 500 m resolution¹¹⁵, focusing on the PFTs that represent boreal forests. The deciduous needleleaf PFT covers 17.2% or 1.97 × 10⁶ km², particularly in Eurasian Dwc and Dfd zones (Extended Data Fig. 2).

Third, these data were intersected with continuous permafrost boundaries, excluding discontinuous, sporadic and isolated permafrost zones⁶⁶. Continental-scale differentiation was maintained throughout the analysis, with 5 clusters identified in North America and 11 in Eurasia (Extended Data Table 2 and Extended Data Fig. 2). This separation accounts for fundamental bioclimatic differences: Eurasian forests experience stronger continentality with colder winters, while clusters in North America show greater precipitation variability. Within each climate zone, clusters are only established for PFTs that are present.

Processing of remote sensing and climate forcing data

We used Google Earth Engine¹¹⁸ to analyse the spatial distribution of tree PFTs within continuous permafrost regions and their associated Köppen climate classifications. The approach included three steps: (1) data acquisition and preprocessing, (2) masking and classification, and (3) spatial analysis and export.

We obtained three primary datasets: (1) the Köppen–Geiger climate classification⁵³, (2) a continuous dataset on the extent of permafrost⁶⁶ and (3) the MODIS land cover product (MCD12Q1, Collection 6.1, ref. ¹¹⁵). The MODIS dataset provides annual land cover classifications at a spatial resolution of 500 m. We used the 2020 land cover map and extracted the ‘LC_Type5’ band to identify boreal PFTs. The classification identifies areas with ENF, DNF and DBF trees >2.0 m in height and with tree cover >10%. The use of MODIS PFTs to map forested areas is limited by multiple factors that warrant careful consideration. Foremost is the coarse spatial resolution, which obscures fine-scale forest patterns and is inadequate for detecting forest clearing occurring at subpixel scales or mixed pixels in transition zones. In addition, class-specific accuracies vary widely, particularly for spectrally similar vegetation types, resulting in potential misclassification and omission of forest patches due to residual noise in input time series, missing data, and variability in the training database. Alternative land cover classifications and datasets would yield somewhat different estimates of boreal forest extent on continuous permafrost, depending on the classification approach and spatial resolution (Supplementary Section 4.4).

To define forested areas within continuous permafrost zones, we applied a binary mask to the PFT data to retain pixels classified as ENF (class 1), DNF (class 3) and DBF (class 4). The continuous permafrost dataset (1,000 m spatial resolution) was rasterized to match the MODIS grid and used to create a binary mask (1 = continuous permafrost), which was then applied. We note that use of the Obu et al. (2019) permafrost dataset introduces uncertainty due to its model-based representation of permafrost extent, limited ground validation and static snapshot nature, which oversimplifies complex spatial and temporal dynamics⁶⁶.

A combined mask was subsequently applied to the Köppen climate classification raster, filtering climate data to include only areas overlapping with both continuous permafrost and tree PFTs. The processed rasters were clipped to predefined bounding geometries for North America and Eurasia before export. The final output consisted of a two-band image: one band representing the continuous permafrost filtered tree PFT distribution and the other the corresponding Köppen climate classification. Resulting combinations with fewer than 10 pixels were filtered out to reduce noise, resulting in a total of 16 clusters (Extended Data Fig. 2). ERA5 hourly forcing (1990–2023) was spatially averaged per cluster (Supplementary Section 7).

Canopy density was classified into three classes based on mean summer (June, July, August) LAI derived from MODIS MOD15A2H (Collection 6.1, ref. ¹¹⁹): sparse canopy (LAI = 1, representing LAI <1.2; 53% of area), moderate canopy (LAI = 2, representing 1.2 ≤ LAI ≤ 3; 31.1%) and dense canopy (LAI = 4, representing LAI >3; 15.9%). Full details of LAI processing and quality filtering are provided in Supplementary Section 8.

To characterize ground ice conditions, we utilized the Northern Hemisphere Permafrost Ground Ice Volume dataset (GGD318, version 2)¹¹⁶, which provides estimated excess ice content at 12.5 km spatial resolution (Supplementary Fig. 6). The dataset classifies excess ice into five categories based on volumetric ice content and distribution patterns: negligible (<10% by volume), low (10–20%), medium (20–30%), high (>30%) and variable (heterogeneous distribution with localized high-ice-content areas). Excess ice distribution across forested continuous permafrost regions showed considerable spatial heterogeneity. We filtered for high excess ice conditions, indicative of ice-rich permafrost susceptible to ground subsidence upon thaw.

Process representation in the coupled permafrost–vegetation model

To simulate the thermal and hydrological processes within the clusters, we applied the CryoGrid community model, a numerical land surface model designed for permafrost environments¹²⁰. CryoGrid integrates modules for excess ground ice, snow physics, lateral water flow and a multilayer forest canopy, making it well suited for capturing ecosystem-specific interactions between forest and permafrost. The model set-up follows previous implementations^25,121 with updated adaptations to account for PFT-specific canopy characteristics (Extended Data Fig. 3 and Supplementary Tables 4 and 5).

Canopy radiative transfer

The multilayer canopy module resolves shortwave and longwave radiation transfer between atmosphere, canopy layers and ground surface¹⁷. Shortwave radiation is calculated using a two-stream approximation that partitions solar radiation into direct beam, scattered direct beam and diffuse components for sunlit and shaded leaves based on cumulative plant area index and leaf optical properties (reflectance and transmittance). This represents the summer shading mechanism that reduces solar radiation reaching the ground by up to 90% (ref. ¹⁸). Longwave radiation exchange among canopy layers, the ground surface and the atmosphere is represented using a multilayer radiative transfer scheme that accounts for both emission and scattering. This framework resolves canopy-scale longwave radiative trapping, whereby radiation emitted from the surface and vegetation is absorbed and re-emitted within the canopy, enhancing downward longwave flux at the forest floor relative to open conditions¹⁷. The resulting increase in subcanopy longwave irradiance elevates ground surface temperatures, particularly during winter months¹⁸.

Turbulent heat fluxes

The multilayer canopy model computes sensible and latent heat fluxes at multiple within-canopy layers using a diffusive flux approach based on scalar concentration gradients and eddy diffusivity profiles. Above-canopy fluxes are calculated using Monin–Obukhov similarity theory with stability corrections modified for the roughness sublayer. Within the canopy, wind speed and turbulent diffusivity both decrease exponentially with depth below the canopy top, with the rate of decrease determined by a mixing length scale that relates to canopy density¹⁷. This vertical structure creates distinct microclimates within the canopy, with temperature, humidity, and wind speed varying by layer. The ground surface energy balance is solved separately from canopy layers, with ground-level turbulent fluxes computed from the temperature and humidity differences between the ground surface and the lowest atmospheric layer within or below the canopy, mediated by aerodynamic resistance that depends on near-surface wind speed. The canopy effectively buffers ground-atmosphere coupling by reducing wind speeds and vapour pressure gradients at ground level. During cold periods, this suppression reduces turbulent cooling of the ground (warming effect), while during warm periods it limits turbulent heating from warm air above (cooling effect). This multilayer approach allows the model to represent how canopy density (LAI) modulates the degree of ground thermal buffering.

Evapotranspiration and soil moisture coupling

Canopy transpiration is calculated from stomatal conductance optimized for water-use efficiency within constraints imposed by plant hydraulics and soil-to-leaf water transport^17,122. Transpiration fluxes are subtracted from soil layers within the rooting depth, and evaporation/sublimation fluxes from the ground or snow surface are calculated based on surface energy balance¹⁸. Soil thermal conductivity is calculated as a weighted power mean of the conductivities and volumetric fractions of water, ice, air, mineral and organic constituents¹²³, following its implementation in CryoGrid^124,125. Increases in soil moisture enhance thermal conductivity, strengthening conductive heat exchange between the surface and deeper soil layers and thereby altering ground heat flux. Soil moisture dynamics are represented in a one-dimensional vertical scheme that accounts for canopy interception, throughfall, evapotranspiration and precipitation infiltration^17,18.

Snow dynamics

Snow processes are represented using the Crocus multilayer snow scheme¹²⁶ as implemented in CryoGrid¹²⁷, which simulates snow accumulation, compaction, densification, microstructure metamorphism and melt. Snow microstructure is characterized by grain size, sphericity and dendricity, which evolve on the basis of temperature gradients and liquid water content¹²⁶. The canopy module includes snow interception on canopy elements, handled similarly to liquid precipitation interception, which reduces ground snowpack accumulation. Sublimation is calculated from the snow surface based on surface energy balance. Canopy shading reduces incoming shortwave radiation to the snowpack, modulating the timing and rate of snowmelt. The snowpack provides thermal insulation between the atmosphere and ground, with notable implications for winter ground temperatures and subsequent summer thaw depths¹⁸.

Subsurface heat and water transfer

The subsurface extends to 100 m depth where the geothermal heat flux is set to 0.05 W m⁻² (ref. ¹²⁸). The ground is divided into layers with thickness of 0.05 m for the top 8 m, 0.1 m for the next 20 m (8–28 m depth), 0.5 m from 28 to 50 m, and 1 m thereafter²². Heat conduction is solved with thermal conductivity and heat capacity that respond dynamically to changes in moisture content, ice content and temperature following ref. ¹²³ and ref. ¹²⁵. Flow in freezing soil is solved using a modified non-isothermal Richards equation. Phase partitioning among liquid water, water vapour and ice follows ref. ¹²⁹, which eliminates jump discontinuity at freezing and improves performance for unsaturated conditions typical of boreal eastern Siberia. ALT is diagnosed as the maximum depth of seasonal thaw where temperatures exceed 0 °C during the thaw season.

The model set-up follows previous implementations^22,120 with updated adaptations to account for PFT-specific canopy characteristics (see Extended Data Fig. 3 and Supplementary Table 4 for parameters and constants) based on MODIS-derived PFT classifications and ref. ¹⁷. PFT is specified as a parameter at initialization which automatically assigns canopy parameters, including leaf and photosynthetic properties. The canopy is described by LAI, stem area index and leaf area density function. Static canopy configurations are prescribed with summer LAIs of 1.0, 2.0 and 4.0 m² m⁻², stem area index of 0.05 m² m⁻² and tree height of 12.0 m (ref. ¹²¹). The vertical distribution of leaf area is represented by the beta distribution probability density function with parameters p = 3.5 and q = 2 for needleleaf trees, which approximates a cone-like tree shape¹³⁰. Deciduous phenology (DNF and DBF) is prescribed through seasonal LAI variations, with leaf-off periods from October to May when LAI is reduced to minimum values representing stems and branches only, while ENFs maintain constant LAI year-round. The lower atmospheric boundary layer is simulated by 4 m of atmospheric layers above the canopy.

Model initialization and spin-up

We simulate the surface and ground fluxes in forested and bare-ground set-ups for each cluster. Simulations are performed with a 5-min time step to resolve subdaily variations in energy fluxes and soil moisture. Snowpack evolution follows the Crocus snow scheme¹²⁷, ensuring representation of snow accumulation, interception and melt processes¹²⁶. Our modelling framework isolates the thermal effect of canopy cover by comparing forested versus bare-ground scenarios within identical climate-soil-permafrost combinations. For each of the 16 bioclimatic clusters, we apply identical soil parameterizations (mineral layer properties) and identical initial soil temperature and moisture profiles to both forested and unforested model runs, varying only the presence and absence of canopy and the organic layer thickness. For forested simulations, the model prescribes a thin organic layer (0.16 m), typical for upper soil layers below canopies, whereas for bare simulations, a combined organic–mineral soil profile is applied, with a thin organic layer of 0.04 m, following previous model set-ups (Supplementary Table 6). The bare-ground simulation serves as the reference run for unforested conditions. Simulations are run for 1991–2023 and analysed for the 2018–2023 analysis period. This paired design allows us to attribute differences in ALT specifically to canopy-mediated thermal buffering. The model set-up follows previous implementations^22,120 with updated adaptations to account for PFT-specific canopy characteristics (Supplementary Table 4) based on the MODIS-derived PFT classifications¹⁷. Future studies might also consider clustering based on soil composition and organic layer thickness.

Target variables and spatial analyses

To evaluate the physical impact of boreal forest cover on permafrost, we quantified the thermal buffering effect of canopies across varying densities on ALT and the resulting impact on SOC exposure to seasonal thaw.

ALT analysis

We employed a delta approach that combines observational baseline data with model-derived canopy effects. For each of the 16 bioclimatic clusters, we simulated four scenarios using CryoGrid: three forested scenarios with varying canopy density (sparse, LAI = 1; moderate, LAI = 2; dense, LAI = 4) and one bare-ground scenario (LAI = 0). From these simulations, we calculated the cluster-specific change in ALT for each canopy density:

where ({{rm{ALT}}}_{{rm{bare,sim}}}) is the simulated ALT under bare-ground conditions, ({{rm{ALT}}}_{{rm{LAI}},{rm{sim}}}) is the simulated ALT for each forested scenario (LAI = 1, 2 or 4) and the subscript ‘cluster’ denotes that each value is specific to a given climate–PFT combination. Positive values of ΔALT_LAI,cluster indicate that the canopy reduces ALT relative to bare ground, with the magnitude representing the thermal buffering capacity of each canopy density.

To generate spatially explicit estimates of ALT under hypothetical canopy removal scenarios, we applied these cluster-specific ΔALT values to the observational baseline ALT from the Climate Change Initiative (CCI) permafrost product⁶⁷: ALT_bare,est(x, y) = ALT_CCI(x, y) + ΔALT_LAI,cluster where ALT_CCI(x, y) represents the gridded observed ALT under current forested conditions at location (x, y), and ALT_bare,est(x, y) is the estimated ALT if the current canopy were replaced with the corresponding LAI scenario. Each LAI scenario (1, 2 or 4) produces a distinct ΔALT_LAI,cluster value and thus a different estimated bare-ground ALT for each pixel. The ΔALT_LAI,cluster value is applied uniformly to all pixels within each cluster. This approach leverages the spatial detail of the CCI observational product while incorporating model-derived estimates of how canopy density modulates thermal buffering effects.

We evaluated CryoGrid’s capability to reproduce observed ALT magnitudes using field measurements from the Circumpolar Active Layer Monitoring (CALM) programme¹³¹. The CALM dataset comprises 17 sites within four of the simulated clusters (Supplementary Fig. 1). Model validation demonstrates that CryoGrid reproduces the observed range and spatial patterns of ALT in forested continuous permafrost regions.

SOC analysis

The canopy-protected SOC pool was quantified by comparing SOC exposed to seasonal thaw under current forested conditions versus hypothetical bare-ground conditions for each LAI scenario. For each pixel, we calculated: ΔSOC_LAI(x, y) = SOC_thawed[ALT_bare,est(x, y)] − SOC_thawed[ALT_CCI(x, y)], where SOC_thawed[ALT] represents the vertically integrated SOC from the surface to ALT, calculated from the Northern Circumpolar Soil Carbon Database (NCSCDv2)¹⁵. Because each LAI scenario produces different ALT_bare,est values, this yields three estimates of canopy-protected SOC corresponding to sparse (LAI = 1), moderate (LAI = 2) and dense (LAI = 4) canopy scenarios.

The integration accounted for vertical SOC heterogeneity by applying depth-weighted interpolation across the four standardized, available depth intervals (0–0.3 m, 0.3–1 m, 1–2 m and 2–3 m). For each pixel, the absolute ALT depth (ALT_bare,est or ALT_CCI) determined which depth intervals contributed to the seasonally thawed SOC pool. The resulting ΔSOC_LAI(x, y) represents the additional SOC that would be exposed to seasonal thaw under each LAI scenario compared with current conditions, corresponding to the canopy-protected carbon pool at each canopy density. Finally, the total net change in exposed SOC for each scenario was determined by spatially summing the pixel-level ΔSOC_LAI(x, y) values across the forested continuous permafrost domain, yielding regional estimates of canopy-protected carbon stocks in petagrams.

Raster calculations were performed in two stages. Initial computations of SOC exposure to thaw were conducted in a geographic coordinate system (EPSG:4326) at 0.0449° × 0.0449° resolution matching the CCI and NCSCDv2 products. For mass summation, the SOC rasters were reprojected to an equal-area polar stereographic system (EPSG:3413) at 500 m resolution to ensure accurate area-based calculations of carbon stocks in petagrams. Raster operations were executed using GDAL (Geospatial Data Abstraction Library, v3.8.0)¹³² with array-based computations in NumPy (v1.24.1)¹³³ within a Python 3.9.5 environment. Maps were generated using QGIS version 3.34.1-Prizren¹³⁴. For large-scale analyses exceeding available RAM, we employed block-wise processing with optimized tile sizes (1,024 × 1,024 pixels), sequentially processing data subsets while maintaining spatial consistency through edge handling and memory flushing between iterations.

This assumed only gradual thaw, where only the carbon within the active layer is considered at risk of mobilization. In regions characterized by high excess ice based on ref. ¹¹⁶, permafrost thaw is likely to occur rapidly due to thermokarst and thaw lake formation processes (Supplementary Fig. 6). Therefore, in a second analysis step and for high excess ice areas only, we considered the entire SOC column from the surface down to 3 m to be exposed to seasonal thaw upon potential canopy removal (Supplementary Fig. 4).

The spatial distribution of the exposed SOC was analysed across dominant soil classifications (Extended Data Fig. 4) to identify vulnerable pedological units. For soil classification, we integrated the four primary permafrost-affected soil types from the NCSCDv2 (Histels, Orthels, Histosols and Turbels)¹⁵ with Yedoma deposits (high confidence only) from the IRYP v2 database¹¹⁷ (Supplementary Fig. 6). The NCSCDv2 dataset acknowledges substantial uncertainties (Supplementary Section 4.4), particularly regarding spatial coverage of deep pedons and distribution of massive ground ice¹⁵, while the IRYP v2 applies a three-level confidence classification (confirmed, likely, uncertain) to account for variability in Yedoma identification¹¹⁷. By restricting our analysis to high-confidence Yedoma areas only, we minimize classification uncertainty in this particularly vulnerable soil type. Percentage-cover rasters for each NCSCDv2 soil class were used to derive a dominant soil type per grid cell, defined as the class with the highest fractional coverage. This dominant-class approach introduces classification uncertainty in mixed-soil pixels, where secondary soil types may occupy substantial fractions but are not represented in our analysis. Yedoma deposits were treated as a priority class and assigned dominance wherever present, regardless of underlying soil composition. The Yedoma raster was resampled to match the spatial extent and resolution of the NCSCDv2 grid using nearest-neighbour interpolation before classification. The resulting dominant soil classification raster was reprojected from geographic coordinates (EPSG:4326) to the equal-area polar stereographic system (EPSG:3413) at 500 m resolution to align with the SOC data. A masking procedure was applied to retain only pixels with valid SOC values and to restrict the analysis to continuous, forested permafrost regions. For quantitative analysis, SOC was aggregated by soil class using zonal statistics. Total carbon stocks were calculated by summing pixel-level SOC density (kg m⁻²) multiplied by pixel area (m²), with results reported in petagrams. Area calculations and percentage distributions were computed to assess the relative contribution of each soil class to total exposed SOC.

SOC vulnerability classifications were assigned to each dominant soil type. These classifications represent pedological characteristics rather than quantitative thaw rates, and inherit the spatial and classificatory uncertainties of the underlying soil datasets. Yedoma deposits were classified as very high vulnerability due to deep carbon storage, high lability of well-preserved organic matter, and high excess ice content that facilitates rapid thermokarst-mediated carbon mobilization^71,75. Turbels were classified as high vulnerability because cryoturbation distributes fresh organic matter throughout soil profiles, making carbon accessible upon thaw under variable drainage conditions¹⁵. Orthels were classified as medium vulnerability because, despite lower total SOC content, post-thaw drainage favours rapid aerobic decomposition⁷². Histels and Histosols were classified as low vulnerability because waterlogged conditions persist after thaw, constraining decomposition to slower anaerobic pathways regardless of depth¹⁵. These classifications reflect the rate and pathway of potential SOC mobilization following canopy loss and permafrost thaw, rather than absolute carbon stock size.

Soil class composition differs markedly across PFTs (Supplementary Fig. 7). DNF forests, which dominate the continuous permafrost zone, overlie predominantly Orthels (59.9%) with a substantial Turbel fraction (28.0%) and the highest Yedoma occurrence of any PFT (6.0%), concentrating both medium- and high-vulnerability soils. ENF and DBF forests show a strikingly different pattern, with Turbels as the largest class (53.8% and 48.8%, respectively) and Histels as the second most abundant (29.8% and 40.7%), reflecting wetter, more cryoturbated soil conditions and very low Yedoma occurrence.

To estimate the biomass carbon pools we combine the different above- and below-ground live biomass, litter biomass and below-canopy biomass pools from literature, using low and high regional estimates for North America and Eurasia each (Supplementary Section 6).