Northeast Siberian Arctic Lowlands (NESAL) permafrost landscapes are shaped by past climatic, geomorphic, and ecological processes. Climate-driven ground ice accumulation and melting are key factors, creating history-dependent landscapes with varying ground ice distributions. The present-day ground ice distribution significantly influences permafrost thaw pathways and future landscape evolution. Thawing of ice-rich permafrost, particularly the melting of massive ground ice, causes thermokarst—landscape change resulting in characteristic landforms. In continuous permafrost, thermokarst is seen in transitions between polygon types or the formation of thaw lakes and gullies, creating landscape-scale feedbacks on hydrology and carbon decomposition. Unlike gradual thawing in ice-poor terrain, thermokarst can lead to rapid permafrost degradation within years or decades. The NESAL, with its abundance of ice- and organic-rich permafrost deposits, is highly susceptible to thermokarst. These deposits, particularly within the Yedoma domain, are estimated to store around 100 GtC, equivalent to a significant atmospheric CO2 increase if released. Despite recent warming trends, NESAL permafrost temperatures remain cold. Earth System Models (ESMs) generally project NESAL permafrost as stable beyond 2100, even under RCP8.5, but these models typically neglect thermokarst processes, which are likely to be significant due to abundant ice wedges. This study uses the CryoGrid 3 model to investigate how present-day ground ice distribution in the NESAL affects future landscape evolution, permafrost degradation, and the release of carbon.

Existing literature highlights the influence of past climate and geological processes on current permafrost landscapes. Studies demonstrate the importance of ground ice distribution in shaping permafrost thaw pathways (e.g., Burn, 1997; Grosse et al., 2007; Schirrmeister et al., 2008; Kanevskiy et al., 2014; Gilbert et al., 2017). The role of thermokarst in accelerating permafrost degradation is well-established (e.g., Jorgenson et al., 2006; Kokelj & Jorgenson, 2013; Olefeldt et al., 2016; Walter Anthony et al., 2018). Research using simpler models provides initial estimates of abrupt thaw's global significance, but lacks detailed physical processes (e.g., Lee et al., 2014; Schuur et al., 2015; Liljedahl et al., 2016). The limitations of current ESMs in representing thermokarst and the consequent underestimation of future permafrost degradation are widely acknowledged (e.g., Turetsky et al., 2019, 2020). Previous studies on the NESAL have focused on its carbon storage potential (e.g., Strauss et al., 2013, 2017) and its susceptibility to thermokarst (e.g., Morgenstern et al., 2011). While ESMs generally project stability for the NESAL under future climate scenarios, this study aims to address this knowledge gap by incorporating detailed representations of thermokarst.

This study employs an extended version of the CryoGrid 3 permafrost model, incorporating processes that induce thermokarst. The model represents spatial heterogeneity in surface topography and subsurface stratigraphy using laterally coupled tiles (polygon centers, rims, and troughs). Three major landscape types are considered, reflecting variations in ice-wedge volumes and thickness: Yedoma deposits (YD), drained lake basins (LB), and Holocene deposits (HD). The model simulates ground subsidence due to excess ice melt, small-scale lateral fluxes of heat, water, and snow, and lateral sediment transport, which is a key process for stabilizing ice wedges after initial degradation. Simulations were conducted under different hydrological conditions (water-logged versus well-drained) and future climate scenarios (RCP2.6, RCP4.5, and RCP8.5). Reference runs without thermokarst-inducing processes are included for comparison. The model classifies landscape states based on the relative positions of soil surface altitudes (relict polygons (RP), low-centered polygons (LCP), intermediate-centered polygons (ICP), high-centered polygons (HCP), and water bodies (WB)). The model uses data on ground ice content, ice wedge dimensions and organic carbon concentrations derived from an extensive dataset of soil samples from the NESAL to represent the sub-surface stratigraphy. The study area is defined as all landmass north of 66°N within the political borders of Yakutia, focusing on ice-rich lowlands covering approximately 493,000 km². Simulations were run for the period from 2000 to 2100. The model evaluates the amount of carbon currently locked in permafrost that becomes susceptible to thaw due to model-simulated active layer deepening and ground subsidence. This process is then scaled to calculate the estimated total amount of carbon that could be affected by thaw across the NESAL, under various climate change scenarios. A sensitivity analysis was performed to test the robustness of the results.

The model simulations demonstrate a substantial increase in permafrost degradation under RCP4.5 and RCP8.5 scenarios compared to reference runs without thermokarst processes. Maximum thaw depths increased significantly, with relative increases ranging from 1.3 to 8.0. Ground subsidence due to excess ice melt contributed significantly to this degradation (0.2m to 4.7m by 2100). Under RCP4.5, stabilization of permafrost was projected towards the end of the 21st century, due to the accumulation of an ice-poor layer preventing the thaw front from reaching deeper layers containing excess ice. However, under RCP8.5, permafrost continued to degrade beyond 2100. The abundance of excess ground ice strongly controlled thaw rate and magnitude; landscapes with high excess ice content showed more rapid and severe degradation. Hydrological regime also played a crucial role, with water-logged conditions leading to deeper thaw and the formation of thaw lakes, whereas well-drained conditions resulted in predominantly unsaturated conditions and a high-centred microtopography. Scaling the simulation results to the NESAL's area, the study estimated the amount of organic carbon becoming susceptible to thaw. Under RCP4.5, 3.2-9.3 GtC could be affected by 2100, substantially more than the 2.7 GtC projected by the reference run. Under RCP8.5, 12.5-64.4 GtC could be affected, significantly exceeding the 5.3 GtC projected by the reference run. Under RCP2.6, the amounts of carbon affected by thaw were comparatively small. The difference between simulations with and without excess ice increases with the severity of the warming scenario. The hydrological conditions influenced not only the amount of thawed carbon but also its decomposition pathways (anaerobic vs. aerobic).

These findings demonstrate that the inclusion of thermokarst processes is crucial for accurately projecting permafrost thaw and carbon release. The significant difference between the model results with and without thermokarst highlights the substantial underestimation inherent in simpler models used in ESMs. The results support the idea of a critical threshold in climate warming; under RCP4.5, stabilizing feedbacks slow degradation and establish a new equilibrium, whereas under RCP8.5, positive feedbacks dominate, leading to continuous and rapid thawing. The significant increase in thaw-affected carbon under RCP4.5 and RCP8.5 compared to the reference simulations emphasizes the necessity of incorporating thermokarst processes into ESMs for a more reliable assessment of the permafrost carbon-climate feedback. The strong influence of both excess ice content and hydrology emphasizes the need for detailed spatial and temporal analysis of these factors for accurate projections. The results are particularly relevant for cold, ice-rich permafrost regions, where thermokarst may significantly increase the amount of carbon released compared to gradual thaw alone.

This study demonstrates that thermokarst significantly accelerates permafrost thaw in cold, ice-rich regions like the NESAL. Simulations show a much larger potential for carbon release than indicated by ESMs that ignore these processes. The results underscore the necessity of including thermokarst in ESMs for accurate projections of permafrost carbon-climate feedback. Mitigation of climate change is crucial to limit the impacts on permafrost ecosystems and prevent the widespread landscape collapse projected under high-emission scenarios. Further research should focus on incorporating detailed spatial variability in topography and hydrology, as well as meso-scale hydrological and geomorphological processes, to improve the accuracy of future projections.

The model simplifies some aspects of landscape complexity. The model does not explicitly resolve all processes occurring at meso-scale (e.g., thermo-erosional valleys and their lateral interactions), and the impacts of extreme weather events are not fully accounted for. Although the use of multiple soil samples and area estimations increases the generality of the findings, a certain degree of uncertainty remains due to inherent variability in the NESAL’s geological and hydrological parameters. The lack of a biogeochemical model prevents a complete assessment of carbon decomposition and greenhouse gas fluxes. Further refinement of these processes would enhance the model's predictive capacity.