Summary of the modeling approach

Permafrost-region 21st-century soil carbon emissions are compared between two model types using IPCC RCP4.5 and RCP8.5. Both models represent basic sets of permafrost processes and have multiple soil organic carbon pools, but the models differ in their complexity of how individual processes are described. While CLM4.5BGC7 simulates the full physical interactions between the atmosphere and the soil, including vegetation uptake of CO 2 , AThaw parameterizes soil thaw rates, depending on ground thermal properties, mean annual ground temperatures, active layer depth, and magnitude of the regional warming anomaly which drives permafrost degradation23. In its current model design, CLM4.5BGC simulates gradual thaw in terrestrial uplands. We utilized CLM4.5BGC emission data for years 1950–2100 partitioned according to non-permafrost carbon, originating from present-day active layer horizons, and permafrost carbon, which becomes thawed from the top down as active layer gradually deepens. In contrast, the more simplistic, but therefore more flexible model design in AThaw also allows for abrupt thaw under thermokarst lakes by tuning model parameters to simulated talik growth rates of a physically-based thermokarst-lake model16. Newly-formed thermokarst lakes and laterally-expanding lake margins are a large net source of atmospheric CH 4 and CO 2 , in contrast to mature thermokarst-lake stages, where fluxes are lower (Supplementary Tables 5 and 6). AThaw simulates only the carbon emissions from newly thawed sub-lake sediments comprising contributions from expansion of thaw lakes present since 1850, and from new lake initiation in response to warming after the year 1850. AThaw does not simulate fluxes for older lake areas already present on the landscape at year 1850. Further, permafrost degradation in CLM4.5BGC is driven by spatially resolved climatic forcing fields, while AThaw focuses on large-scale latitudinal climatic gradients. In addition to differences in spatial resolution, CLM4.5BGC is more complete, simulating a large set of climate variables which all determine surface vegetation and soil conditions, while AThaw considers surface air temperatures only as the key driver for the net balance of lake expansion versus drainage.

AThaw model

AThaw is a conceptual model which projects carbon release from abrupt thaw by accounting for the full chain of processes from formation of new thermokarst lakes under global warming and talik deepening in sub-lake sediments to eventual carbon release to the atmosphere following anaerobic microbial degradation of organic matter and CH 4 oxidation. AThaw also accounts for lake drainage; although, carbon fluxes in drained lake basins are not modeled24,63,64,65. AThaw is incorporated into a multi-box permafrost-carbon release model which allocates soil organic matter into latitudinally and vertically gridded boxes of differing conditions regarding soil physics, carbon quantity and quality, and biogeochemistry23.

Briefly, AThaw model resolution consists of (a) 20 latitudinal bands, ranging from 45°N to 85°N with a 2° gridding, and of (b) 27 vertical soil layers corresponding to layer thicknesses of 25 cm for the upper 4 m, and of 1 m thickness for the depth range of 4–15 m. AThaw assumes typical soil organic matter residence times in permafrost soils, which are determined by partitioning into passive, slow, and fast cycling carbon pools with decomposition timescale parameters based on incubation experiments. Q 10 temperature sensitivity is accounted for as well as CH 4 oxidation [for details see Supplementary Table 1 and Schneider von Deimling et al.23.

AThaw assumes that increasing Arctic temperatures will drive expansion of existing lakes and new thermokarst-lake formation by melting of near-surface ground ice and subsequent ground subsidence. This assumption is in line with Community Land Model (CLM4.5) results from Lee et al.66, whereby surface excess ice in permafrost soils of many regions will largely melt by 2100 when subject to intense warming. AThaw assumes that comparatively small 19th and 20th-century warming rates have initiated some formation of new but rather shallow thermokarst lakes with likely winter-refreeze of lake bottom sediments. Start of abrupt thaw (i.e., formation of sub-lake taliks) is only assumed to occur for stronger warming starting in the 21st century and beyond. The evolution of newly-formed thermokarst lakes is parameterized in AThaw by an optimum function which non-linearly scales the latitudinal thermokarst lake area fraction by the surface air temperature anomaly23 (Supplementary Figs 4 and 5). The lake formation scheme models an increase of newly-formed thermokarst-lake areas with temperature until a maximum extent [FTKLmax (~8 to ~40% increase, depending on soil type)], is reached under a temperature optimum dT′TKLmax. The temperature optimum corresponds to high-latitude surface air temperatures 4–6 °C above pre-industrial; warming above this optimum shifts 21st century lake dynamics toward net drainage. Such responses of thermokarst lake formation to temperature increases are not unprecedented, given that the early Holocene Arctic temperature increase of 1.6 °C37 resulted in a 570% increase in thermokarst-lake formation rates24. AThaw-projected increases in lake area are lower for 2010–2100 compared to the early-Holocene observational record, because permafrost soils are now more protected from warming and thawing by thick organic soil surfaces67 and less likely to form large lakes compared to the less dissected early Holocene periglacial landscape52.

Temperature-driven thermokarst lake dynamics

Rather than simulating individual lake life cycles of formation, expansion, drainage, and re-initiation of later generation lakes (e.g., ref. 43), AThaw quantifies the net effect of new lake formation and drainage in modelled lake areas (FTKL). We capture a wide range of uncertainty in net lake formation and drainage trajectories by varying two key model parameters: The maximum net lake area, FTKLmax, and the optimum high latitude surface temperature increase, dT′TKLmax (Supplementary Fig. 5). Model parameters for the maximum lake area extent (FTKLmax) were chosen individually for the four AThaw soil classes to capture expected differences in the potential for future lake development (see Methane and CO 2 in newly-formed thermokarst lakes). However, the parameter that most strongly controls the dynamics of thermokarst lake formation and drainage is the temperature optimum, dT′TKLmax, the temperature at which the maximum lake area occurs. We prescribed a mean estimate of 5 °C for dT′TKLmax (i.e., high latitude surface air warming above pre-industrial) and consider a full range of 4–6 °C in our model ensemble. This parameter choice is based on paleoenvironmental evidence of Early Holocene warming by a few degrees Celsius in Northern Hemisphere land areas37,68,69 which resulted in rapid and intensive thermokarst activity24,29,70.

This range of dT′TKLmax values was also chosen to capture the sensitivity of future surface ice-wedge melt to climate warming, since ice-wedge melt is a critical step in thermokarst lake formation and drainage. Increases in permafrost temperature, which typically mirror increases in air temperature, have been observed in many Arctic regions, with warming of up to 3 °C since the 1970’s in some of the coldest permafrost regions3,38,71. Permafrost warming is often accompanied by active layer deepening and increased ground-ice melt. Widespread surficial degradation of ice wedges and a significant increase in areas of water-filled polygonal troughs has been linked to climatic warming during the past few decades in northern Alaska, the Canadian Arctic Archipelago, and Siberia14,15,41. Additional permafrost warming on the order of 2–3 °C anticipated by 2050 for RCP4.5 and RCP8.549, is expected to intensify permafrost and surface ice wedge degradation66, enhancing thermokarst. For RCP4.5, permafrost warming slows beyond 205049, supporting AThaw ensembles of lower net thermokarst lake formation during the latter part of the century (Supplementary Figs 4 and 5). In contrast, extreme warming after 2050 for RCP8.5 is expected to heavily degrade permafrost, such that near-surface permafrost disappears entirely from many arctic regions49. These conditions lend support to AThaw parameterization, which leads to net lake drainage in the later part of the century for RCP8.5.

In AThaw, warming above the temperature optimum in the 21st century is assumed to lead to a reduction in AThaw lake area due to increasing lake drainage23,43,45, until a prescribed minimum fraction FTKLmin remains. The minimum fraction of the landscape still covered by newly formed lakes, decreases from north to south23. In northern, continuous permafrost regions, lakes in AThaw are prescribed to drain laterally as melting of the ice-wedge network in the surface surrounding lakes can create drainage pathways52. Other mechanisms of lake drainage include bank overtopping, headward gully erosion towards a lake; and tapping of lakes by streams, rivers, or other water bodies36. Higher drainage potential in southern, discontinuous permafrost regions also encompasses internal drainage of lake water through open taliks that penetrate thin permafrost in groundwater recharge settings54.

While the 4–6 °C dT′TKLmax range used in our modeling was prescribed based on paleoenvironmental evidence, historical observations, and modeling of future permafrost dynamics, we can consider the implications of using smaller or larger dT′TKLmax values. Smaller dT′TKLmax values would imply that the maximum area of thermokarst lake coverage would occur earlier in the century. This would result in a relatively larger AThaw increase to the PCF earlier in the 21st century and a lower AThaw contribution to PCF toward the end of the century due to more widespread drainage. If dT′TKLmax aligned more closely with the range of simulated high latitude warming inferred from higher concentration pathways (e.g., dT′TKLmax > 6 °C), we would expect maximum thermokarst lake coverage (i.e., FTKLmax) to occur later in the century, less drainage of AThaw lakes during the 21st century, and a larger relative contribution of AThaw to end-of-the century PCF under stronger future warming.

It is interesting to note that that while thermokarst-lake initiation and drainage rates are linked to climate, thermokarst-lake growth—the long process between initiation and drainage that results in most of the carbon release—has dynamics (e.g., talik growth, shore expansion) that once started become rather decoupled from climate due to strong linkage with local factors such as ground ice content and ice-layer thickness72,73,74,75. Hence, thermokarst lakes are found across the entire range of Arctic climatic zones and permafrost temperatures as long as there is sufficient ground ice36,76,77. A thermokarst lake on the New Siberian Islands has the same potential to release carbon as a thermokarst lake in Central Yakutia; the differences are largely not determined by climate (or RCP conditions) but by local conditions such as permafrost soil carbon and ground ice contents. Hence, either RCP scenario will result in more lakes (earlier or later) and both scenarios will have a similar emission magnitude linked to maximum lake areas, but relative to anthropogenic emissions, the thermokarst lake emissions from RCP4.5 will be more relevant.

Factors other than temperature that are not included in AThaw, but which can also affect thermokarst lake dynamics include natural and anthropogenic surface disturbance, precipitation changes, and local topography43. However, since CMIP5 models consistently predict a moistening trend in the Arctic (i.e., an increase in precipitation relative to evapotranspiration)51, an explicit accounting for predicted 21st-century precipitation and evaporation trends would reinforce rather than weaken the AThaw lake dynamics driven by temperature changes alone.

Methane and CO 2 in newly-formed thermokarst lakes

AThaw simulates pan-Arctic CH 4 (CH 4TKL ) and CO 2 (CO 2TKL ) release from thermokarst lakes proportional to the amount of newly-thawed carbon that becomes vulnerable to microbial decomposition. The volume of this newly-thawed carbon expands vertically by talik growth in sub-lake sediments and horizontally by increases in the extent to which thermokarst lakes cover the landscape. Here we discuss how these two processes, vertical and horizontal growth, are captured in AThaw and how they compare to other modeling studies and observational evidence.

First, AThaw simulates vertical talik growth rates beneath lakes as a function of atmospheric temperature anomalies, which in turn determine lake bottom temperatures and ultimately drive sub-lake sediment warming. Thaw rates are assumed to depend on soil ice contents, mean ground temperature, depth of the thaw front, and are tuned to reproduce talik deepening simulated by Kessler et al.16. For instance, AThaw simulates typical talik depth beneath lakes of 11 m (7.8–13.4 m, 68% uncertainty range; 20 m max) in warm permafrost environments (i.e., mean annual soil temperatures close to 0 °C) with mineral soils by the year 2050 under RCP8.5. Our helicopter-borne electromagnetic measurements of sites with comparable climatic and soil conditions (i.e., relatively warm permafrost temperatures between −3 and −1°C) (Supplementary Information) have inferred lake talik growth of 8–15 m in < 50 years (Supplementary Fig. 7). It should be noted that advective heat transport by groundwater can accelerate vertical thaw; maximum reported thaw rates are up to a meter per year78,79. When accounting for the effect of ground subsidence in a modeling context, Westermann et al.79 suggest that several meters of ground ice can be removed in less than a decade. The authors simulated a talik growth of 15 m by 2100 under RCP8.5 for a site characterized by yedoma ice complex in cold, continuous permafrost. In contrast, AThaw simulates 6 m (4–8 m, 68% range) talik depth under comparable climatic conditions by end of the century. Therefore, given modeling and observational evidence, we consider that the thaw rates in sub-lake sediments simulated by AThaw are relatively conservative.

A second key AThaw model assumption concerns simulated thermokarst-lake expansion in a warmer climate. Central to this aspect is the question to which extent newly-formed thermokarst lakes will cover degrading permafrost landscapes in a warmer climate. AThaw assumes that future thermokarst-lake formation will strongly depend on soil conditions (e.g., ice content) and landscape morphology. Therefore, the model assumes different maximum lake formation extents for four different ice-rich soil type distributions in the permafrost region of the Arctic: mineral soil (Orthels and Turbels), organic soils (Histels), undisturbed yedoma, and refrozen thermokarst deposits in the yedoma region (Supplementary Table 1, Supplementary Figs 4 and 5). Orthels, Turbels and Histels follow the Northern Circumpolar Soil Carbon Database (NCSCD)80, while undisturbed yedoma and refrozen thermokarst deposits in the Pleistocene-aged ice-rich yedoma soil region are distinguished according to Strauss et al.81 and Walter Anthony et al.24.

For mineral soils, AThaw assumes that newly-formed lakes can degrade a maximum fraction FTKLmax of 7% (4–9%, 68% range; Supplementary Fig. 4) of the landscape. Given the large-scale dominance of mineral soils in the permafrost region soil carbon inventory (540 Pg C, NCSCD)80, this soil class contributes significantly to cumulative pan-arctic CH 4TKL (at year 2100, 46% for RCP8.5, 28% for RCP4.5). Organic soils (120 Pg C, NCSCD)80 are typically richer in ground ice than mineral soils and therefore AThaw assumes a factor two larger potential for the formation of new thermokarst lakes [FTKLmax of 14% (10–17%, 68% range; Supplementary Fig. 4)]. This soil type contributes 17% for RCP8.5 and 14% for RCP4.5 to cumulative pan-arctic CH 4TKL . AThaw explicitly accounts for ice- and organic-rich soils in yedoma regions, separated into undisturbed yedoma landscapes (80 Pg C)81 and drained lake basins in the yedoma region that subsequently formed permafrost following lake drainage (refrozen thermokarst; 240 Pg C)24,81. Given the high ice contents of these yedoma-region soils81,82,83, both soil classes are assigned a high potential for the formation of new thermokarst lakes once rising temperatures have resulted in ground subsidence. Given the high volumetric ice contents of refrozen drained thermokarst lake basins (> 50%)81,83 in conjunction with its basin-type geomorphology, which favors water ponding16, AThaw assumes a maximum FTKLmax of 21% (11–27%, 68% range; Supplementary Fig. 4). Refrozen thermokarst basins in the yedoma region cover only about 4% of the permafrost domain of the Arctic, but they contribute 28% under RCP8.5 (38% under RCP4.5) to cumulative pan-arctic CH 4TKL by 2100. The highest potential for new-lake formation is assumed for unaltered ice-rich yedoma soils, where volumetric segregated ice content is typically 47–53% and high wedge-ice volumes (~40%) further increase permafrost ice concentration81. Here AThaw assumes a FTKLmax of 33% (16–42%, 68% range; Supplementary Fig. 4). Considering that thermokarst activity in ice-rich regions had degraded about 80% of the landscape in specific regions during the Holocene24,81,84, we consider a factor two reduced lake-forming potential plausible, especially under the assumption of strong future warming. While undisturbed yedoma landscapes constitute our carbon pool with the smallest areal extent (covering 0.41 million km2, 2% of the permafrost domain)81, they contribute 9% under RCP8.5 (20% under RCP4.5) to cumulated pan-arctic CH 4TKL by 2100, revealing the importance of carbon release from deep deposits beneath abruptly-formed yedoma thermokarst lakes.

We derived estimates of present-day abrupt-thaw emissions from AThaw and independently from upscaling observations. The estimated emission range from AThaw [20 (7–49) Tg C-CO 2 e yr−1] represents the median and 68% uncertainty range for years 2011–2017 (Supplementary Data 1). The observation-based estimate (19–58 Tg C-CO 2 e yr−1) was derived by upscaling observed lake carbon fluxes from abrupt-thaw features (Fig. 2b) to the extent of observed gross lake expansion areas (5–15%) among pan-arctic regions during the past 60 years21. This range of gross lake expansion is comparable to our observation of gross lake-area increase in northern and western Alaska (i.e., 1999–2014, 1.1–1.7% gross lake area increase upscaled to 5–7% assuming similar expansion/formation rates over the past 60 years; Supplementary Table 3, Fig. 3). We acknowledge that present-day emissions associated with abrupt thaw may be conservative since gross lake-area increases in some regions are higher (e.g., > 50% in Quebec85 and > 23% in West Siberia86), potentially due to hydrological changes in addition to permafrost thaw.

The formation/expansion of new thermokarst-lake areas and subsequent sub-lake talik growth is not the only mechanism of abrupt permafrost degradation making newly thawed permafrost carbon available for microbial decomposition. Further contribution comes from talik growth of present-day lakes which have not yet formed a deep talik and still store large amounts of labile carbon in sub-lake sediments. For instance, on the Alaska North Slope, Arp et al.17 demonstrated that a rapid decrease in lake-ice thickness and duration has already led to many shallow lakes transitioning from bed-ice fast lakes underlain by permafrost to floating ice lakes that have started to develop taliks. This talik formation takes place about 70 years before talik formation is projected for the adjacent terrestrial environment by top-down permafrost models in this cold continuous permafrost zone, a feedback process also not accounted for in AThaw or other models of permafrost degradation. These and other non-lake modes of abrupt thaw13, including coastal and river erosion, thermoerosional gully formation, thaw slumps, collapse of permafrost peat plateaus, and talik formation in upland terrestrial environments, which may occur sooner than predicted by large-scale models when finer resolution soil and vegetation properties are taken into account49, are not explicitly accounted for in the AThaw model description; hence, we consider our assumed FTKLmax values a conservative estimate of the abrupt thaw extent in the 21st century.

Uncertainty

Uncertainty in AThaw modeled carbon fluxes is based on independent sampling of a set of 18 key model parameters which are subject to either observational or to model description uncertainty23 (Supplementary Table 1). For each warming scenario (RCP4.5 and RCP8.5), 500 ensemble runs were performed by applying a statistical Monte Carlo sampling and by assuming uniformity and independence in the distributions of model parameters and initial values. AThaw results are presented as the median and 68% uncertainty range. We acknowledge that our ability to accurately quantify uncertainties is limited given the use of this single model and its assumptions in a highly complex and large system.

Double counting

To avoid double counting CLM4.5BGC emissions from land areas that become thermokarst lakes in AThaw, and to account only for the increase in emissions and CPCRE caused by abrupt thaw, we have subtracted from the AThaw emissions and CPCRE, in our total permafrost landscape calculations, the quantity of carbon emissions and associated radiative forcing already assumed to be emitted from those land surfaces from CLM4.5BGC. We used the AThaw thermokarst-lake area fraction at each time step for each of the four soil classes (Supplementary Fig. 4) and weighted those fractions by the areal extent of the soil classes according to Hugelius et al.80 and Strauss et al.81 based on the implicit assumption of homogenous soil carbon distribution within each of the soil classes. Our calculations consider that this fraction of the land surface, subject to gradual thaw in CLM4.5BGC, undergoes abrupt thaw instead.

Additional methods

We calculated the radiative forcing due to atmospheric perturbations in CH 4 and CO 2 concentration for CLM and AThaw permafrost-soil-carbon flux trajectories following Frolking & Roulet87 (Supplementary Methods). Methodology underlying our remote-sensing based quantification of abrupt thaw in Alaska and field-based estimates of carbon emissions from abruptly-formed thermokarst areas of lakes in Alaska and Siberia are also provided in Supplementary Methods. Field and lab measurements include bubble-trap observations of ebullition fluxes, aerial and ground-based ebullition seep-mapping, and quantification of CH 4 and CO 2 concentrations and radiocarbon dating.

Statistics

To test differences between gradual thaw and abrupt thaw net ecosystem exchange (Fig. 2b) we used the two-sided Kolmogorov-Smirnov test. We used this test also to compare radiocarbon ages of permafrost soil carbon respiration in gradual thaw versus abrupt thaw environments (Fig. 2c) Radiocarbon statistical analysis was performed on percent modern carbon data. All statistical analyses were performed in R88.

Data availability

CLM data are publicly available at the National Energy Research Scientific Computing Center archive (https://www.portal.nersc.gov/archive/home/c/cdkoven/www/clm45_permafrostsims/permafrostRCN_modeldata). AThaw data and calculated radiative forcing for CLM and AThaw fluxes presented in this study are available within the article’s supplementary information file (Supplementary Data 1). Radiocarbon data are provided in Supplementary Data 2. All other data that support the findings of this study are available from the corresponding author (K.M.W.A.) upon request.