Fast response of cold ice-rich permafrost in northeast Siberia to a warming climate

Permafrost landscapes of the northeast Siberian Arctic lowlands (NESAL) have been shaped by climatic, geomorphic, and ecological processes over the Late Quaternary, leading to heterogeneous ground-ice distributions that critically influence present-day thaw pathways. Ice-rich permafrost, particularly wedge ice in Yedoma and Holocene deposits, is highly susceptible to thermokarst, which can trigger rapid, abrupt thaw and landscape change compared with gradual thaw in ice-poor terrain. Despite very cold present permafrost temperatures (−8 to −12 °C) and only recently observed warming in boreholes (~0.9 °C per decade), Earth system models (ESMs) have projected NESAL permafrost to remain largely stable beyond 2100, even under RCP8.5. However, these models lack key thermokarst-inducing processes (excess ice melt, ground subsidence, lateral heat/water/snow fluxes), likely underestimating future degradation and carbon feedbacks. This study extends a process-based model (CryoGrid 3) to include these processes and evaluates how present-day ground-ice distributions, inherited from landscape history, control future landscape evolution, permafrost degradation pathways, and exposure of frozen organic carbon under RCP2.6, RCP4.5, and RCP8.5 and contrasting hydrological conditions. The aim is to quantify the magnitude and pace of degradation, assess stabilizing versus destabilizing feedbacks, and estimate thaw-affected carbon stocks to inform carbon-climate feedback assessments and mitigation importance.

Prior work shows that ground-ice distribution governs thermokarst pathways and future evolution of ice-rich permafrost landscapes, with ice-wedge polygon degradation and thaw-lake formation driving rapid (abrupt) thaw. Conceptual models suggest global relevance of abrupt thaw but rely on simplified assumptions. Current ESMs lack representation of excess ice, ground subsidence, and lateral mass/energy fluxes, likely underestimating degradation and carbon release. The NESAL overlaps extensively with the Yedoma domain, which contains deep, ice- and organic-rich permafrost with potential tipping behavior under warming. Observations report pan-Arctic ice-wedge degradation and hydrological impacts, rapid thermokarst development even in very cold permafrost, and large permafrost carbon pools that may be mobilized. Previous modeling advances introduced subgrid heterogeneity, excess-ice melt, and coupled tiles to capture microtopography and lateral fluxes, motivating the process-based approach adopted here.

Model: An extended CryoGrid 3 land-surface model represents polygonal tundra with three laterally coupled tiles (polygon centers, rims, troughs), each with 1-D vertical thermohydrological columns (heat conduction with phase change, unfrozen-ground hydrology). The snow scheme simulates accumulation, melt, infiltration, and refreezing. Lateral fluxes of heat, water, and snow occur between tiles assuming circular symmetry. Troughs can drain to an external reservoir with specified elevation to impose hydrological conditions (water-logged vs well-drained). Geomorphic processes: (1) Excess ground ice scheme simulates ground subsidence from melting excess ice; excess water moves upward, sediments are routed downward, representing ice-wedge degradation. (2) A nonlinear hillslope diffusion–based lateral sediment transport scheme (adapted for periglacial settings) is implemented as slope-dependent advective sediment flux between tiles, activated by steep gradients after subsidence. Parameters include a critical slope angle αcrit = 45°, subaerial and subaqueous transport coefficients Kland = 3e−10 m² s−1 and Kwater = 3e−8 m² s−1. Landscape types and stratigraphy: Three representative landscape types are simulated with distinct surface microtopography and cryostratigraphy: drained lake basins (LB; low-centred polygons), Holocene deposits (HD; low-centred polygons), and Yedoma deposits (YD; initially flat relict polygons). Ice-wedge geometry and ground-ice content differ by type (Table values used in model): wedge depths ~3.8 m (LB), 10 m (HD), 20 m (YD); wedge-ice volumes 10%, 30%, 50%; excess ice contents 13%, 19%, 25%; ground-ice contents 68%, 74%, 80%. An ice-rich intermediate layer 0.2 m (LB, HD) to 0.4 m (YD) overlays wedges. Organic/mineral/ice contents are derived from 984 soil samples across NESAL sites. Forcing: Meteorological forcing (1901–2100) for the central Lena River delta: downscaled CRU-NCEP v5.3 (1901–2014) plus CCSM4 anomalies for RCP4.5 and RCP8.5 (2015–2100). RCP2.6 anomalies are applied by repeating 2000–2014 RCP4.5 period with RCP2.6 anomalies. Climate across NESAL is assumed comparable to the Lena delta for this purpose. Simulations: For each landscape type, runs are performed under RCP2.6, RCP4.5, and RCP8.5 with two hydrological conditions (water-logged vs well-drained, set by reservoir elevation). Initial profiles come from boreholes and a 50-year spin-up (10/1949–12/1999). Analysis covers 01/2000–12/2099. Reference runs emulate typical ESM land schemes: 1-D columns without excess ice and without lateral fluxes. Diagnostics: Microtopographic states are classified by relative tile elevations and reservoir level into relict polygons (RP), low-centred polygons (LCP), intermediate-centred polygons (ICP), high-centred polygons (HCP), or water bodies (WB). Maximum annual thaw depth is the area-weighted mean across tiles; ground subsidence and fractions of saturated/unsaturated thawed cells are tracked annually. Evaluation and sensitivity: Previous work validated CryoGrid 3 tile set-up against observations at Samoylov Island (soil temperatures, thaw depth, moisture, water table, snow, energy fluxes). A sensitivity analysis varied fresh snow density, sediment transport coefficients, field capacity, and natural porosity; results are provided in Supplementary material. Scaling to region: NESAL area delineation combines Yedoma mapping and DEM-based criteria for lowland polygonal tundra (altitude 5–100 m a.s.l., slope ≤4°), yielding areas: LB ~140,000 km², HD ~290,000 km², YD ~63,000 km²; total NESAL ~493,000 km². Thawed organic carbon per unit area from simulations is scaled by landscape areas to estimate regional thaw-affected carbon stocks.

- Thermokarst processes substantially increase projected permafrost degradation relative to models without excess ice and lateral processes. Maximum thaw depths increased during 2000–2100 by factors from 1.3 (well-drained LB, RCP4.5) to 8.0 (water-logged YD, RCP8.5). Reference runs (no excess ice) showed lower increases (factors ~1.7 under RCP4.5 and 2.3 under RCP8.5). - Ground subsidence from excess-ice melt adds 0.2 m (well-drained LB, RCP4.5) to 4.7 m (water-logged YD, RCP8.5) of additional degradation by 2100, a process absent in simplistic models. - Landscape evolution depends strongly on climate scenario and hydrology. Under RCP2.6, landscapes largely remain stable, with limited changes (shallow water bodies in water-logged YD). Under RCP4.5, initial ice-wedge degradation occurs but stabilizing feedbacks (sediment redistribution and formation of an ice-poor layer) tend to establish a new equilibrium by late century, especially under well-drained conditions. Under RCP8.5, widespread landscape collapse is simulated: polygon rims/subsidence produce HCPs and form surface water bodies under water-logged conditions. Water-body depths by late century reach ~1 m (LB), 2 m (HD), and 4 m (YD), with taliks ~3–5 m thick. - Yedoma deposits show the strongest response due to high excess-ice content: water bodies form within 2–3 decades under water-logged conditions; the landscape turns into a water body by 2100 under both RCP4.5 and RCP8.5, with talik formation under RCP8.5. Under well-drained conditions, HCPs emerge after 6–7 decades; under RCP8.5, massive subsidence produces conical thermokarst mounds (baidzharakhs). - Hydrological regime tips with ice-wedge degradation: under water-logged settings, thaw subsidence leads to pervasive inundation and saturated conditions; under well-drained settings, subsidence of rims enhances drainage, producing predominantly unsaturated thawed soils, especially where HCP microtopography develops. - Thaw-affected organic carbon stocks scale strongly with inclusion of thermokarst processes and warming level. By 2050, additional thaw-affected carbon is similar between approaches due to active-layer deepening: ~1.7 GtC (RCP2.6), 1.8 GtC (RCP4.5), 2.5 GtC (RCP8.5). By 2100: • RCP2.6: 0.8–2.0 GtC (with excess ice) vs 1.3 GtC (reference). • RCP4.5: 3.2–9.3 GtC (with excess ice) vs 2.7 GtC (reference) — up to about three-fold. • RCP8.5: 12.5–64.4 GtC (with excess ice) vs 5.3 GtC (reference) — up to about twelve-fold; potentially up to ~two-thirds of NESAL’s ~100 GtC pool becomes thaw-affected. - The deviation between thermokarst-enabled simulations and reference runs increases with warming strength and depends on hydrology, with water-logged cases generally producing deeper thaw and larger thaw-affected carbon. - Results imply that ESMs likely underestimate thaw-affected carbon, especially in cold regions like NESAL that appear stable when only gradual thaw is represented (reference runs show only ~5.3% of NESAL carbon thawed under RCP8.5).

The study demonstrates that including thermokarst-inducing processes (excess ice melt, ground subsidence, lateral fluxes, microtopography evolution, and sediment transport) alters projections of permafrost stability and carbon exposure in ice-rich lowlands. The findings address the central question by showing that present-day ground-ice distributions inherited from landscape history set distinct thaw pathways and rates: landscapes with greater excess ice (e.g., Yedoma) experience more rapid and severe degradation and transition to water bodies compared with less ice-rich terrains. A qualitative regime shift appears between RCP4.5 and RCP8.5: under moderate warming, stabilizing feedbacks (sediment infill forming an ice-poor protective layer) can slow and eventually balance thaw by late century; under strong warming, positive feedbacks (snow accumulation in depressions, higher thermal conductivity in wetter active layers, persistent surface water) outpace stabilization, leading to sustained degradation and landscape collapse. Hydrological connectivity changes with rim subsidence reorganize surface and subsurface water regimes, producing saturated, lake-prone states under water-logged conditions versus enhanced drainage and predominantly unsaturated soils under well-drained settings. These hydrological states have implications for decomposition pathways (anaerobic vs aerobic), though biogeochemical fluxes are not quantified here. Scaling results indicate that models lacking thermokarst processes likely underestimate the fraction of deep permafrost carbon that becomes thawed, especially in cold regions where gradual active-layer deepening alone is modest. Incorporating these processes into ESMs is therefore crucial for realistic assessments of the permafrost carbon–climate feedback, and the results underscore the benefits of climate mitigation (RCP2.6–RCP4.5) in limiting degradation and carbon exposure.

A process-based representation of thermokarst in a permafrost model reveals that cold, ice- and organic-rich NESAL landscapes can respond rapidly to warming, with widespread ice-wedge degradation, significant ground subsidence, and formation of water bodies, especially under RCP8.5. Under RCP4.5, stabilizing sedimentation feedbacks can moderate thaw by late century, whereas under RCP8.5 degradation continues and landscape collapse is widespread. Accounting for thermokarst processes increases projected thaw-affected carbon by up to three-fold (RCP4.5) and up to twelve-fold (RCP8.5) relative to simplistic representations, potentially exposing up to two-thirds of NESAL’s carbon to thaw under high emissions. These results emphasize the need to include thermokarst-inducing processes and microtopography–hydrology feedbacks in ESMs to improve projections of permafrost stability and carbon–climate feedbacks. Future research should couple physical models with biogeochemical schemes to quantify greenhouse gas production and fluxes, refine hydrological and geomorphic interactions at meso-scales (e.g., lake expansion and drainage, snow redistribution), and incorporate spatial variability of ground-ice distributions and landscape history to reduce uncertainty in regional scaling.

- Existing, old thermokarst lakes and large water bodies were excluded; the model only simulates new water bodies formed via ice-wedge thermokarst. - Meso-scale geomorphic and hydrologic processes (e.g., thermo-erosional valleys, lake drainage/expansion, blowing snow redistribution) are not explicitly represented; hydrology is simplified via fixed reservoir levels. - The approach targets decadal-to-centennial scales; long-term ground-ice accumulation during colder periods is not modeled and could offset degradation on millennial scales. - Biogeochemical processes (microbial decomposition pathways, greenhouse gas production and emission) are not simulated; thus, thaw-affected carbon does not directly translate to net carbon fluxes. - Climate forcing is taken from the Lena River delta and assumed representative across NESAL; spatial climate variability is not resolved. - Parameter and structural uncertainties remain (e.g., sediment transport coefficients, snow density, soil hydraulic properties), and sensitivity analyses indicate potential variability in outcomes. - Regional scaling assumes stratigraphies and carbon contents from sampled sites are representative of each mapped landscape class; mapping criteria may omit or include areas with differing properties.