Panarctic lakes exerted a small positive feedback on early Holocene warming due to deglacial release of methane

Climate-driven permafrost thaw can release ancient carbon to the atmosphere, begetting further warming in a positive feedback loop. Polar ice core data and young radiocarbon ages of dissolved methane in thermokarst lakes have challenged the importance of this feedback, but ﬁ eld studies did not adequately account for older methane released from permafrost through bubbling. We synthesized panarctic isotope and emissions datasets to derive integrated ages of panarctic lake methane ﬂ uxes. Methane age in modern thermokarst lakes (3132 ± 731 years before present) re ﬂ ects remobilization of ancient carbon. Thermokarst-lake methane emissions ﬁ t within the constraints imposed by polar ice core data. Younger, albeit ultimately larger sources of methane from glacial lakes, estimated here, lagged those from thermokarst lakes. Our results imply that panarctic lake methane release was a small positive feedback to climate warming, comprising up to 17% of total northern hemisphere sources during the deglacial period.


T
he Arctic is warming four times faster than the entire globe 1 , and predicted warming will continue this trend.A projected global temperature increase of 4 °C by 2100 (equivalent to 8 °C in the Arctic, SSP5-8.5) 2 threatens permafrost region soils, which cover 24% of the northern hemisphere (NH) land surface and store twice as much carbon as the atmosphere 3 .As permafrost warms and thaws, soil organic carbon previously sequestered from the atmosphere for thousands of years becomes available for microbial decomposition, resulting in enhanced release of carbon dioxide and methane (CH 4 ).The projected release of these greenhouse gases from thawing permafrost (92 ± 17 Pg C by 2100, RCP8.5) 4 poses a positive feedback that could accelerate warming and impede ongoing climate change mitigation and adaptation efforts 5 .
The likelihood of permafrost carbon feedback to climate warming has recently been challenged by new radiocarbon analysis of methane from Antarctic glacial ice 6 .The last deglaciation (18 to 10 ka) was an 8000-year period over which global air temperature increased 6 °C7,8 and atmospheric methane nearly doubled 9,10 .Within that period, two distinct events of abrupt warming at 14.7 ka (Bölling-Alleröd warming) and 11.6 ka (Preboreal warming) were associated with abrupt (decade-tocentury scale) increases in atmospheric methane concentration (AMC).During these periods of abrupt AMC rise, Dyonisius et al. 6 found little change in the radiocarbon age of global atmospheric methane.They concluded that since methane release from 14 C-depleted carbon sources, such as old permafrost and hydrates, responded minimally to deglacial warming, these sources will be unimportant in scenarios of future warming.This conclusion was corroborated by other recent studies reporting contemporary radiocarbon ages of methane dissolved in Arctic thermokarst lakes 11,12 ; however, these studies did not account for the relatively older 14 CH 4 age of lake bubbles, particularly those formed in deeper thawing permafrost layers of carbon-rich Pleistocene-aged yedoma soils 13 , the domain of which extended across ~2.6 million square kilometers of north Siberia, Alaska, and northwest Canada during the late Pleistocene 14 .
Here we examine the contribution of permafrost carbon sources to the AMC during the last deglaciation using an independent approach that is more direct than inferring northern permafrost sources from analyses of Antarctic ice.We synthesized 154 published radiocarbon ages of diffusive and ebullitive methane emissions from present-day panarctic lakes 15 and applied these ages to modeled lake methane emissions based on an extensive new paleoenvironmental dataset of 1207 northern lake initiation ages 16 .This yielded a bottom-up reconstruction of lake methane emissions and their radiocarbon values from 18 ka to the present (Methods), which we then compared to top-down estimates based on ice core models.
Thermokarst lakes constitute the most widespread form of abrupt thaw in the panarctic 17 , and the most efficient means of rapidly (over decades to centuries) mobilizing old carbon by deep permafrost thaw [18][19][20][21] .However, in addition to thermokarst lake emissions, we also calculated emissions from glacial lakes since they (1) emit 14 C-depleted methane 19 ; (2) dominate panarctic lake area (i.e., glacial lakes comprise >75% of all lakes north of 45°b y both number and area 22 ); and (3) represent a cryosphere response to warming 23 .Our new dataset increases the number and type of lake initiation records previously used to reconstruct deglacial lake methane emission dynamics 24,25 by >10-fold (Fig. S1).In addition to 162 thermokarst lake initiation dates, it includes 921 glacial lake records, and 124 dates from lakes of other or unknown origin.Altogether, these lake initiation data span the entire spatial domain characterized by past permafrost and glaciations during the Last Glacial/Permafrost Maximum in the northern hemisphere 26 (Fig. 1).This domain expands well beyond the present-day permafrost extent and represents extensive lake-rich regions that were not included in past thermokarst lake methane reconstructions for Beringia 24,25 .We hereafter refer to the lakes within this study domain as "panarctic."To our knowledge, this is the first study to specifically estimate methane contributions from glacial lakes distinctly from boreal peatlands, thermokarst lakes, and other northern (>30°N) sources.Aggregated, all northern sources comprised up to 30% of the global atmospheric methane source during the deglacial period 27,28 .
Despite inescapable uncertainties associated with reconstructing past phenomena, we here present estimates of deglacial to present-day panarctic lake methane emissions based on the best available data.We have attempted to address probable uncertainties where possible, and discuss additional sources of error that we were unable to quantify due to their complexity or a lack of data (Supplementary Information Section 1).We acknowledge uncertainty with respect to the absolute magnitude and timing of estimated emission values, but maintain that our data support important conclusions regarding the function of thermokarst lakes as a positive feedback to deglacial climate warming.We furthermore highlight factors we suspect may have had a potentially significant influence on past lake emission strengths, and introduce opportunities to resolve data gaps that currently prevent better model accuracy.

Results
Patterns of lake formation.Across the entire deglacial period, climate fluctuations impacted landscape-level deglaciation, permafrost degradation, and the associated lake-forming processes 23 .Panarctic lake formation pre-dated rapid peatland expansion 29 by ~2000 years and followed an increase in NH summer insolation 30 and temperature inferred from Greenland ice δ 18 O 31 (Fig. 2).Two peaks in panarctic lake formation (13.2 ka and 10.4 ka, Fig. 2b) were separated by a brief decline synchronous with cooling during the Younger Dryas (12.9-11.6 ka).Thermokarst Lake formation, specifically, appears to have been driven by increases in temperature and changes in water balance towards the wetting of permafrost landscapes, which accelerated and amplified ground ice melt 23 .Our dataset reveals that thermokarst lake formation increased quickly and strongly following the onset of both warming periods (though more so with sustained Preboreal warming) (Fig. 2c), suggesting that permafrost responded rapidly to changing climate conditions.Glacial lake formation, dependent on regressive melt of glaciers and ice sheets with high thermal inertia, responded more gradually to climate fluctuations in both directions (Fig. 2c).Widespread glacial lake formation was likely driven by ice sheet recession that exposed large areas of the land surface to associated lake-forming processes such as inundation of glacially carved depressions, water impoundment by depositional features, melt of buried remnant glacial ice, and basin isolation by isostatic adjustment of the Earth's crust 23 .
Subsequently, total panarctic lake formation decreased dramatically between 10 ka and 8 ka despite sustained early Holocene warmth (Fig. 2b).In permafrost regions, cumulative fragmentation of original, ice-rich uplands by thermokarst lake formation, modification, and drainage, likely increased overall landscape relief and drainage capacity throughout the deglacial period 23 .Development of a highly incised landscape 32 would have strongly controlled lake size and persistence over time, since forming or expanding lakes would have been increasingly likely to drain via pre-existing slopes or channels 18,33 .As lake formation decreased from 10 to 8 ka, major ice sheets also approached their minimum extents 34,35 .The deceleration of glacial lake formation during the early Holocene is therefore attributed to the evolution of landscapes from high-energy proglacial environments to more Fig. 1 Panarctic map of the cryosphere and lake basal ages.The panarctic domain (poleward of the solid black line) is defined as the combined regions affected by the Last Permafrost Maximum permafrost extent 26 (unhatched), the Last Glacial Maximum (LGM) ice sheets 70 (hatched), and exposed continental shelves during LGM 82 (green).Within this domain, the modern continuous permafrost extent (aqua) and non-continuous permafrost extent (blue) 69 .Locations of basal age records (circles) are color-coded by age.
Fig. 2 Panarctic lake and peatland formation relative to climate indicators from 18 ka to present.Average summer insolation 30 (gray lines) at high northern latitudes and Greenland ice core δ 18 O H2O 31 (purple) (a).Binned and cumulative distributions of lakes 16 (dark gray, n = 1112) and boreal peatland 29 (n = 3612, light gray) basal ages (b).Formation rates for glacial (blue) and thermokarst (red) lake types 23 (c).Methane emissions estimates with temporal uncertainty shown by dual curves A and B, (this study; scenarios explained in Methods) 81 and corresponding 1σ error (SD) envelopes for glacial (blue) and thermokarst (pink) lake types (d).
stable periglacial environments 23 .Rates of panarctic lake establishment continued to decrease steadily from 8 to 4 ka, and were much lower from 4 ka to present (Fig. 2b).
Lake methane emission estimates.Combining our lake initiation dataset with lake-type specific areal methane flux rates scaled to reflect changing temperatures 36 , and constraining total allowable permafrost carbon mineralization in thermokarst lakes following a mass-balance approach 37 (Methods), we found that peak thermokarst lake emissions of 15 ± 6 Tg CH 4 yr −1 preceded peak emissions from glacial and other non-thermokarst lake types (hereafter simplified to "glacial lakes"; 15 ± 5 Tg CH 4 yr −1 at 7 ka) (Fig. 2d).Thermokarst lake emissions comprised the majority of the panarctic lake flux during portions of the deglacial period.However, glacial lakes, which dominate the panarctic lake area today, and which were not included in previous estimates of northern lake emissions since the LGM 24,25 , ultimately exceeded thermokarst lake emissions by between 10 and 9 ka.Emissions from glacial lakes were more than double those from thermokarst lakes over the Holocene (Fig. 2d).Combined, the peak in panarctic lake emissions (23-28 Tg CH 4 yr −1 between 10.5 and 9 ka) preceded peak northern peatland emissions 29 by ~8000 years; however, maximum emissions from panarctic lakes appear to have lagged abrupt increases in atmospheric methane at 14.7 ka and 11.6 ka by at least 500 years (Fig. 3).While methane release resulting from a rapid expansion of lake area following the Younger Dryas (Fig. 2) accounted for up to 55% of the increase in AMC attributed to northern sources over several millennia (12-9.5 ka) 27,28 , partitioning of latitudinal sources across periods of abrupt AMC rise is still poorly resolved.Sharp increases in the AMC preceding gradual increases in both northern lake and peatland emissions, imply that other methane sources may have driven these abrupt changes.
Our estimate of present-day panarctic lake emissions (16-18 Tg CH 4 yr −1 ; Fig. 3f) is similar to total northern lake methane emissions independently estimated from modern lake flux measurements applied to present-day lake areas (7-25 Tg CH 4 yr −1 ) 38 and estimates utilizing satellite data for parameterizing ice-free season length (13.8-17.7 Tg CH 4 yr −1 ) 39 .
Composite age of modern lake emissions.The radiocarbon age of methane emitted from lakes is a function of emission mode and lake type.In open-water lake areas, major emission modes include ebullition and diffusion, as well as the time-delayed release of both ice-trapped bubbles and dissolved gas at ice-off (Fig. 4).Summertime diffusive fluxes and background ebullition originate largely from surface lake sediments, where predominately contemporary organic carbon substrates fuel methanogenesis when surface sediments heat up in summer.In contrast, point-source and hotspot ebullition seeps represent freephase gas emissions channelized by bubble tubes that plumb methane formed in sequentially deeper (older) sediments 40 .In winter, surface sediments cool, slowing methanogenesis; however, the thermal lag of summer heat propagation to deeper sediments Fig. 3 Lake methane emissions in relation to polar ice core records.Northern high latitude average summer insolation 30 (gray) and Greenland ice core δ 18 O H2O 31 (purple) (a).Atmospheric methane concentrations from GISP2 (green) 9 , WAIS Divide 10 (purple), GISP (gray) 83 , and NEEM 84,85 (blue) ice cores (b).Methane emissions from northern peatlands 29 (gray) (c).Northern (>30°N) source contributions to AMC inferred from the interpolar gradient 28 (green) (d).Discrete northern (>30°N) source AMC contributions (mean ± SD) inferred from the interpolar gradient by ref. 54 (red), ref. 27 (blue), and ref. 28 (yellow) (e).Total methane emissions estimated from all panarctic lakes, with temporal uncertainty shown by dual scenarios A (dotted line) and B (dashed line) with corresponding 1σ (SD) error envelopes (light gray and dark gray, respectively) in this study 81 ; scenarios explained in Methods), (f).Vertical light-green shaded bars highlight short periods of rapid atmospheric methane increase from (b).
allows methanogenesis to continue through winter in talks (thaw bulbs), allowing lakes with abundant point-source and hotspot emission of 14 C-depleted methane to persist throughout the year 15,41,42 (Fig. 4a).In winter, bubbles escaping sediments may be temporarily impeded by seasonal ice cover.Some 14 C-depleted methane diffuses from bubbles resting under the ice into the water column, where its fate is methanotrophy or evasion to the atmosphere during spring ice melt 43 .The rest of the ice-trapped bubbles escape periodically during winter through ice cracks or upon degradation of the ice cover in spring 43 .These modes of emission have been summarized for ~140 panarctic lakes 38,41,44,45 and are presented according to lake type (yedoma vs. non-yedoma thermokarst lakes, glacial lakes) in Table S1.
The radiocarbon age of methane in bubbles has been shown to have a 1:1 relationship with permafrost soil carbon ages 19 and is dependent on the type of geologic substrates in which lakes form and the depth of sediments in which methane is produced.Approximately 64% of total annual emissions from lakes formed in thick (average thickness 25 m) late Pleistocene yedoma permafrost deposits occur via point source and hotspot ebullition.Relatively old 14 C-CH 4 ages observed in yedoma lake pointsource bubbles (13,000 to 42,900 yr BP) reflect methane formed deep within the anaerobic taliks 40 .Since carbon-rich soils surrounding non-yedoma thermokarst lakes and glacial lakes are thinner (usually <3 m) and younger (Holocene-aged) 19 , point-source ebullition is less common, diffusive emissions are relatively more important, and the radiocarbon ages of emitted methane in both bubbles and dissolved gases is substantially younger than in yedoma lakes.
Our synthesis of published methane radiocarbon ages and emission modes (diffusion, storage, ebullition) for panarctic lakes 15 revealed that ebullition-derived fluxes with a mean age of 9030 ± 5559 yr BP (mean ± SD) dominated emissions in thermokarst lakes formed in Pleistocene-aged yedoma permafrost soils (Table 1).Summertime diffusive fluxes, and those at ice-out, were 9 and 2% of emissions from these lakes, respectively, but the 14 C age of dissolved gases in surface lake water was 1713 ± 1407 yr BP in summer, and 1713 ± 3120 yr BP in spring under ice, yielding a net annual radiocarbon age of total emissions of 7828 yr BP (Table 1).Comparing seasonal ages of methane dissolved in the surface water of interior Alaska yedoma lakes, summer ages were younger (1490 yr BP) than spring ages (2226 yr BP) because bubbles containing older methane contribute to the dissolved gas pool below ice in winter 43,46 and summertime heating of surface sediments enhances methanogenesis from younger substrates 40 .
Diffusive fluxes comprise proportionately more of annual emissions from glacial (20%) and non-yedoma (30%) thermokarst lakes compared to yedoma lakes (9%).In these lakes, the 14 CH 4 ages of both bubbles and dissolved gas were younger than in yedoma lakes (Table 1), yielding net 14 CH 4 ages of 1849 yr BP and 1984 yr BP for glacial and non-yedoma thermokarst lakes, respectively.As with yedoma lakes, 14 CH 4 was younger in the surface water in summer than in spring, under the lake ice.Weighting emissions by lake area and emission modes, the net ages of total panarctic lake emissions and thermokarst lake emissions were 1978 ± 1226 yr BP and 3132 ± 731 yr BP, respectively (Table 1).These ages represent a mixture of methane formed from contemporary and old substrates, but their net age being thousands of years older than modern carbon, indicates the importance of ancient carbon mobilization and release to the atmosphere.
Radiocarbon mixing model results.To compare our bottom-up estimates with top-down constraints from atmospheric Δ 14 CH 4 measurements in ice cores, we re-ran the atmospheric CH 4 box model of Dyonisius et al. 6 in a forward model scheme (Methods) using our reconstructed thermokarst lake emissions as input.We assumed that thermokarst lakes constituted the majority of emissions from permafrost systems because they are the most efficient means for the deep thaw of 14 C-depleted permafrost substrates 20 .For the 14 C activity of the thermokarst lake emissions, we used both the composite age of all thermokarst lake emissions (3132 yr BP) and yedoma thermokarst lake emissions (7828 yr BP) reported in this study.Our use of the two end members (3132 yr BP and 7828 yr BP) accounts for the impact of any spatiotemporal fluctuation in emitted methane radiocarbon ages over time.Figure 5 shows the expected atmospheric Δ 14 CH 4 during the last deglaciation, given our bottom-up thermokarst lake emissions.Though contemporary emissions appear most compatible, our bottom-up emissions are within the uncertainty Fig. 4 Schematic of methane evasion from three different lake types.Thermokarst Lake situated in ice-rich yedoma permafrost shows a higher ratio of ebullitive (white bubbles) to diffusive (near-surface vertical dashes) fluxes, increasing age of permafrost carbon substrate with thaw bulb (tan area) depth, seasonal propagation of heat into deeper sediments stimulating microbial activity (red curve), and higher rates of methane production near the active thermokarst lake margin (a).Thermokarst Lake is situated in thinner, younger non-yedoma permafrost soil showing a lower ratio of ebullitive to diffusive fluxes, and shallower thaw bulb due to less ground ice content (b).Glacial lake surrounded by non-yedoma permafrost or unfrozen ground showing diffusive and ebullitive fluxes and laminated inorganic and Holocene-aged organic sediments (dark brown) fueling methanogenesis (c).
range of ice core Δ 14 CH 4 data 6,47 when the modern composite thermokarst lake emission age of 3132 years is applied, but nearer the lower limit of the uncertainty range when the older age, representing yedoma-region thermokarst lake emissions, is employed.

Discussion
Lake response to deglacial climate warming.Our study presents the first bottom-up estimate of comprehensive panarctic lake methane emissions since the Last Glacial Maximum.It includes not only thermokarst lakes, which dominated panarctic lake emissions during deglaciation, but also glacial lakes, which subsequently became more abundant and surpassed thermokarst lakes in emissions during the Holocene.Today, the entirety of glacial lakes emits more than twice as much methane as all thermokarst lakes.The northern lake methane source is also clearly distinct from the later-initiating boreal peatland source 29 ; still, it does not account for the northern fraction of abrupt AMC  increases during the Last Deglaciation, and could not unless lakes consistently initiated 500-1000 years earlier than indicated by records in our dataset 16 .This suggests that processes other than the northern lake and peatland formation were responsible for the abrupt AMC rise, and supports the hypothesis that strong tropical sources (30°S-30°N), possibly in combination with smaller, yet-unidentified northern sources, drove abrupt increases in AMC that occurred during the last deglacial period 10,48,49 .Methane produced by inundation and thaw of extralimital permafrost deposits in association with sea level rise 50 , thinning and northward contraction of permafrost areas in response to warming 51,52 , and formation of vast proglacial lakes by ice sheet retreat 34,53 may have comprised additional northern sources.We excluded large ice-dammed, proglacial lakes from our dataset and subsequent lake area and methane emission estimates due to their high spatiotemporal complexity.This exclusion likely makes our estimates of methane emissions during phases of expanding large proglacial lakes conservative; and further study of this potential source is needed.
Variations in the panarctic lake methane flux broadly correspond with changing northern (>30°) source contributions to the AMC implied by the interpolar gradient.A ~25 Tg CH 4 yr −1 increase in the northern source between 12 and 9.5 ka 27,28 , and an ~11 Tg CH 4 yr −1 decrease culminating at 6 ka 28,54 , coincide with changes in panarctic lake emissions.An increase of ~14 Tg CH 4 yr −1 , and a subsequent decrease of ~8 Tg CH 4 yr −1 (Fig. 3) in lake emissions account for 55% and 73% of the changes in AMC attributed to northern (>30°) sources over these multimillennial periods.Later methane emissions resulting from the expansion of northern peatlands likely dominated northern sources for the majority of the Holocene 29 , while panarctic lake emissions continued to stabilize or decrease from 8 ka until present (Fig. 3f).The discrepancy between stable to slightly increasing lake methane emissions and decreasing northern source estimates during the Younger Dryas (Fig. 3) could have been caused by a simultaneous reduction in the strength of a yetunconstrained northern source, persistence of lake formation under cooler climate conditions, or an inadequacy of our methods to account for the full effect of regional climate change on methane production in lakes during the Younger Dryas.Another possible explanation is the effect of dust contamination on the Greenland ice core record, which compromises the reliability of the interpolar gradient during colder climate regimes 55 .
Potential variation of thermokarst lake Δ 14 CH 4 .Lacking actual thermokarst lake Δ 14 CH 4 data for the last deglacial and Holocene, we used two end members (3132 yr BP, representing the age of total modern thermokarst lake emissions, and 7828 yr BP, representing the age of modern yedoma-region thermokarst lake emissions) to account for variation in emitted methane radiocarbon ages over time.Such variation could have arisen from (1) increased mineralization of contemporary and/or accumulating younger carbon sources, including Holocene-aged permafrost, thick peat, and lake sediment packages; (2) a progressive decrease in the spatial extent of unmodified Pleistocene-aged yedoma permafrost due to thermokarst and thermo-erosion 56 ; or (3) changes in relative contributions from yedoma and non-yedoma thermokarst lakes caused by shifts in emission strengths or lake distributions.Today, the composite age of all thermokarst lakes (3132 ± 731 yr BP, Table 1) is significantly younger than yedomatype thermokarst lakes alone (7828 ± 3899 yr BP) because 74% of present-day thermokarst lake area is outside the yedoma region, where thermokarst lakes form in comparatively less carbon-rich and/or younger permafrost soils (Table S2).The two forward modeled atmospheric Δ 14 CH 4 curves based on each end-member age (Fig. 5) bound a range of possible atmospheric Δ 14 CH 4 values for comparison to the ice core Δ 14 CH 4 data.
Thermokarst lakes as a positive feedback to climate warming.While we have tentatively ruled out that lakes were the primary source of abrupt AMC rise, our data do suggest climate-driven lake feedback to warming.Episodes of abrupt warming at 14.7 ka and 11.6 ka (Fig. 3a) triggered ice sheet retreat and permafrost thaw that accelerated both thermokarst and glacial lake formation 23 .Methane release resulting from a rapid expansion in lake area following the Younger Dryas (Fig. 2) likely contributed to a more gradual increase in NH temperature, which followed the initial abrupt rise (11.6 ka, Fig. 3).Radiative forcing of 0.03 W m −2 resulting from a 13 Tg CH 4 yr −1 increase in thermokarst lake emissions would have warmed the planet by ~0.024 °C over the deglacial period.This constitutes a small positive feedback to warming that likely helped sustain high, early Holocene temperatures during a diminishing phase of solar insolation that began 11 ka 30 (Fig. 3).
Contrary to recent suggestions, our findings provide specific evidence that the emission of methane from permafrost carbon sources functioned as a small positive feedback to deglacial climate warming.Thermokarst lake emission ages represent a mixture of methane derived from substrates of different ages (late Pleistocene, Holocene, contemporary, Fig. 4).Despite uncertainties in carbon-source ages and mixing ratios, net emission ages of thousands of years old (Table 1), and methane sample ages that predominantly exceed lake age (Fig. 6), identify thermokarst lakes as a conduit for release of previously sequestered carbon, and distinguish this source from other components of the contemporary carbon cycle.Regardless of how long this carbon was withheld from contemporary cycling, its reintroduction through climate-driven thermokarst lake formation and subsequent methane production constitutes positive feedback to warming.Furthermore, even if all of the 14 C-depleted methane emitted from non-yedoma thermokarst and glacial lakes originated from active layer soil carbon, and not from thawing permafrost soils themselves (false scenario), the release to the atmosphere, especially as methane, still constitutes a permafrost carbon feedback.This is because it represents the mineralization of permafrost region soil organic matter 57 facilitated by the creation of anaerobic environments by permafrost thaw and active layer saturation.
Our bottom-up estimates suggest that thermokarst lakes contributed up to 15 ± 6 Tg CH 4 yr −1 to the atmosphere, and constituted a maximum of 9% of the total northern hemisphere source (compared to 17% for all lakes) during the deglacial period.Comparison with ice core data (Fig. 5) shows that at its current analytical precision, Δ 14 C-CH 4 data from ice cores, unfortunately, cannot provide strong constraints on younger (but still aged relative to contemporaneous 14 CO 2 ) methane emissions from permafrost systems.We showed that methane emissions from lake sources (Figs. 2, 3) had a delayed response relative to abrupt northern hemisphere warming (as recorded in Greenland ice cores 9 ) and peaked at around 11-9 ka.Our bottom-up estimates (Figs. 2, 3) and box modeling results (Fig. 5) provide a target time-period for future ice core measurements (e.g., a more accurate AMC and stable isotopes interpolar gradient that is free from excess CH 4 artifacts from dusty Greenland ice 55 ) to test against.Future studies will have greater success disentangling source contributions as these datasets continue to improve.
While our methods aptly constrain both total carbon loss from yedoma permafrost through the formation of new lakes, and the timing of peak methane release, we found that the distribution of this carbon loss (i.e., relative peak magnitude) during deglaciation was strongly influenced by lake size (Supplementary Information Section 1.2).Changes in thermokarst lake size over time, evidenced by the superposition of smaller lake basins over larger ones throughout the yedoma region 56,58 , likely occurred as icerich permafrost land surfaces evolved from contiguous uplands favoring lateral expansion, to complex, fragmented networks that provided greater opportunity for lake modification and drainage 33,56 .A better understanding of this phenomenon, and quantification of systematic changes in lake size over time, are necessary to better assess the importance of yedoma thermokarst lake methane emissions during the early deglacial period.
Future permafrost carbon feedback.In light of these new paleoenvironmental records, radiocarbon methane age syntheses, and modeling, we conclude that thermokarst lakes were a positive feedback to climate warming during the last deglaciation, albeit small compared to other new atmospheric methane sources at that time.The last deglaciation provides an imperfect analog for modern warming, however.Today, permafrost is less extensive and permafrost landscapes are geomorphically more mature than at the end of the last glacial period.While some areas are more carbon-rich now compared to the end of the last ice age 37,59 permafrost soils are also, in many places, more protected from warming due to insulation by accumulated peat and vegetation 60 .Still, current rates of anthropogenic warming are two orders of magnitude faster than during the previous glacial-interglacial transition; and 50-140% of the 6 °C global warming that occurred over 8000 years from 18 ka to 10 ka 7,8 is projected to occur in the Arctic in just 80 years under current emissions projections 2 .Hence, models agree that 21st-century warming will likely overwhelm contemporary boundary conditions and result in widespread permafrost loss, leading to the mobilization and decomposition of a large fraction of the permafrost carbon pool 4 .Models of abrupt thaw estimate average 21st-century thermokarst lake emissions of 70 Tg CH 4 yr −1 21 , equal to more than four times the maximum emissions from this source during the last deglaciation.Thus, we maintain that thermokarst lakes could again function as conduits for permafrost carbon release that will create a positive feedback to warming this century.Such feedback is already underway in some locations 13 .
Ultimately, much of the 21st-century permafrost carbon feedback will depend on how many new lakes the landscape can support.Models consistently predict an increase in precipitation relative to evapotranspiration in the Arctic, especially in summer 61 , suggesting that hydrological conditions will become more favorable for thermokarst-lake formation 33 .However, some regions are experiencing a net lake area loss 62 due to permafrost thaw and the creation of thermo-erosional drainage networks [63][64][65] or through-going taliks 66 .Since emissions from new thermokarst lake areas are much higher than those from old lakes and drained basis, modeling suggests gross drainage will have to exceed gross lake formation and expansion by a factor of four to seven before the positive permafrost carbon feedback will be reversed 20 .Better predictions of thermokarst lake formation and persistence, and associated methane release, will therefore rely upon enhanced modeling of climate, hydrological and geomorphological interactions under future climate change scenarios.
Future formation of glacial lakes will be constrained to land areas exposed by continued glacier and ice sheet melt.However, increases in methane emissions projected to occur as lake bottom waters warm 67 and ice free-seasons lengthen could be larger for glacial lakes than for thermokarst lakes alone due to their much larger spatial extent.

Conclusion
Here we show that our bottom-up, panarctic lake methane emission estimates are consistent with top-down polar ice core 14 CH 4 constraints, especially when revised lake flux ages based on a new synthesis of radiocarbon data are applied.Lake formation accelerated in response to climate warming, releasing old methane from a mixture of stored (permafrost), and contemporary carbon sources that created a positive climate feedback helping to sustain early Holocene temperature increases.Thermokarst lakes emitted disproportionately more old carbon than glacial lakes.However, new emission estimates from glacial lakes, which dominate the panarctic lake area, were more than double Fig. 6 Ages of thermokarst lake methane samples relative to lake age.The black line depicts a 1:1 relationship between lake age and gas age.Basal lake ages (x-axis) are from Walter ref. 13 and ref. 16 ; Lake methane ages (y-axis) are from and ref. 15 .
those from thermokarst lake emissions throughout the majority of the Holocene.

Methods
Lake methane age calculations.To assess the contribution of old carbon to deglacial and Holocene lake emissions, we compiled extensive radiocarbon data 15 associated with various modes of methane emission (diffusive, storage, and ebullition fluxes) for different Arctic lake types and settings (thermokarst lakes formed in Pleistocene-aged yedoma permafrost of North Siberia, Alaska, and Northwest Canada; thermokarst-lakes formed in Holocene-aged permafrost throughout the remainder of the Arctic permafrost region, and panarctic glacial lakes) (Table S1 and Fig. S2).We used these data, along with a previously developed mixing model 6 based on Antarctic ice core Δ 14 C-CH 4 data 6,47 to reevaluate potential methane contributions from old carbon sources associated with permafrost thaw and lake formation during the deglacial period.
Lake basal age extraction.We used a recent compilation of 1207 basal ages and lake origins from past and extant lakes (HiLLBAO), within the domain of glaciation and permafrost extent during the Last Glacial/Permafrost Maximum (https:// doi.org/10.1594/PANGAEA.894737) 16to reconstruct past lake methane emissions.This dataset was gathered primarily from scientific literature descriptions of lake cores, peat cores underlain by lake sediments, and riverbank exposures.Though 95% of records were radiocarbon-dated, robust dates obtained by a variety of methods were accepted.We evaluated the effect of assuming a normal (1σ = 200 years) or positively skewed (Figs.S4, S5) distribution of dating error for all basal ages included in the dataset on modeled methane emissions (scenarios A and B, respectively, in Figs. 2, 3).We chose this approach to avoid heavy weighting of older radiocarbon dates based on their larger analytical error, and to explore the impact of non-normally distributed dating error on the timing of emission estimates.We acknowledge that basal ages in the dataset could have a bias towards being younger than they really are due to (1) the failure of some coring efforts to penetrate stiff glaciolacustrine sediments overlying true basal sediments; (2) the common interpretation of mineral-organic sediment transitions as indicators of lake formation within the literature; and (3) the inability of a single coring location to account for transgressive lake formation/expansion (i.e., lake formation may have occurred earlier than inundation of the coring site).Applying a positive skew to the distribution of basal age errors allows some assessment of the impact of this possible bias.Other dataset attributes, potential sources of dating error, and spatiotemporal patterns of lake formation are discussed in the original presentation of the dataset 23 .All radiocarbon dates of lake sediments were recalibrated for this study using Calib 7.0 68 .
Modeled lake methane emissions presented here are based on distributions of lake basal ages ≤18,000 years old (n = 1112 out of the full dataset of 1207 dates).We excluded 95 basal ages for lakes formed earlier than 18,000 years ago and large ice-dammed, proglacial lakes from our database and subsequent lake area and methane emission estimates due to their high spatiotemporal complexity.
Modern lake area extraction.We defined the study domain (Fig. 1) as the panarctic region, including all areas that (1) are affected by permafrost today 69 , (2) were affected by permafrost during the Last Permafrost Maximum 26 , or (3) were covered by ice sheets during the LGM 70 .This region, falling primarily northward of 35°to 45°, largely overlaps the high latitude region northward of 30°defined in interpolar gradient studies 27,28,54 .Extracted areas of modern lakes present on both formerly glaciated and unglaciated landscapes within our study domain are based on the Global Lakes and Wetland Database (minimum lake size: 0.1 km 2 ) 71 .We then used the current extent of coastlines and political/geographic or physiographic regions (ESRI Basemap) to separate our study domain into subregions suitable for regional scaling of lake methane emission characteristics (Fig. S3).For each subregion, we determined the areal extents of modern lake distribution based on the Global Lakes and Wetland Database.We used the distributions of different lake types within the dataset 23 to proportionally adjust lake areas by type (Table S2).
Glacial lake methane flux estimates.The majority (88%) of non-thermokarst lakes included in our dataset were formed by a variety of mechanisms directly associated with deglaciation 23 .The remainder (12%) were formed by volcanic, tectonic, carbonate weathering, or other processes.For simplification, we refer to all non-thermokarst lakes as "glacial" lakes in this analysis.To estimate fluxes from these lakes, we adopted a regionalized approach, accounting for strong spatial dynamics in patterns of glacial lake formation related to ice sheet retreat 23 .We used modern lake extents (Table S2) and distributions of glacial lake basal ages 16 to linearly reconstruct glacial lake areas in 100-year timesteps from 18 ka to the present for each subregion of our study domain (Fig. S3).This approach assumes the establishment of differently sized lakes with a consistent mean size (i.e., no systematic variation) over time (Supplementary Information Section 1.2).To partially account for glacial lake area loss caused by climate drying and isostatic rebound since deglaciation, we modified rates of glacial lake formation using a lake drainage function.We derived this function from glacial lake termination ages within our dataset obtained from dated transitions in peat core or sediment exposure stratigraphy (n = 40).This allowed reconstructed glacial lake areas used for scaling methane emissions to exceed modern lake extents during the early Holocene.
Following approaches previously adopted for Holocene peatland emissions 29 , we applied a temperature correction to measured, present-day annual methane fluxes (sum of measured annual diffusion, storage, and ebullition fluxes) from glacial lakes (10.1 ± 4.8 g CH 4 m −2 yr −1 , Table S1) in 100-year timesteps to reflect the ecosystem-level influence of reconstructed changes in temperature 72,73 on methane production using the relationship derived by ref. 36 .We assumed these temperature-corrected areal flux rates were representative across the domain, and multiplied the areal flux rate by the reconstructed lake areas in 100-year timesteps to estimate methane emissions from panarctic glacial lakes through time.
Thermokarst lake methane flux estimates.Thermokarst lakes form when icerich permafrost thaws, causing ground subsidence and the creation of a waterimpounding depression.We used an existing mass balance-constrained approach 37 as the basis for estimating methane emissions from thermokarst lakes in the yedoma region; however, we updated the model to reflect the temporal distribution of 60% more thermokarst lake formation dates available through the High-Latitude Lake Basal Ages and Origins (HiLLBAO) database 16 .Briefly, we used a massbalance derived estimate of 95 ± 34 Pg C 35 to constrain cumulative Pleistoceneaged methane produced in thermokarst lakes from 18 ka to the present.Following long-term laboratory incubation observations and earlier modeling 37,74 , we assumed a 1:1 stoichiometric relationship between carbon dioxide and methane produced from the labile fraction of thawed yedoma.This ratio is also consistent with the theoretical chemistry of steady-state methanogenesis with cellulose as the primary substrate 75 .We then distributed the pool of yedoma permafrost carbon lost as methane in 100-year timesteps to capture end-member scenarios of carbon release trajectories: (1) solely based on rates of new thermokarst lake formation in each timestep following observations that the majority of yedoma-carbon derived carbon is emitted within a century of thermokarst-induced inundation along shorelines 13 (scenario B in Figs. 2, 3) or (2) in combination with a numerically simulated flux trajectory that extends methane produced from thawed yedoma deposits in thermokarst lakes for the first 1200 years after lake initiation following ref. 18(Fig. S6, scenario A in Figs. 2, 3).Though we identified changes in yedoma thermokarst lake size as an important source of potential uncertainty in the timing and magnitude of modeled deglacial lake emissions (Supplementary Information Section 1.2), we could not reliably estimate its impact on our modeling results due to a paucity of data.
We estimated Holocene-aged methane emissions that dominate in late-stage yedoma lakes and non-yedoma thermokarst lakes in 100-year timesteps using modern annual fluxes observed in these environments (6.8 ± 4.4 g CH 4 m −2 yr −1 , Table S1) in the yedoma region 38,44,45 .This value was modified in 100-year timesteps to account for changes in temperature as in methods section 1.4.After 3570 ± 1160 years, lakes were assumed to drain.The 3600-year average "lifespan" of thermokarst lakes was determined by stratigraphic analyses of multiple lake initiation and drainage sequences in yedoma exposures (n = 6 date pairs; data from ref. 37 , ref. 76 and Walter Anthony, unpublished data).Post-drainage emissions from alases (drained lake basins), often fall under a wetland classification.To avoid double counting with peatland methane inventories, we excluded post-lake drainage emissions and emissions from vegetated lake margins from our analysis.Our analysis applies only to open-water lake areas.
We used a similar method to calculate emissions from thermokarst lakes that formed on various other permafrost substrates outside of the yedoma region (e.g., in non-glaciated areas of western Siberia, central Asia, northern Alaska, and Europe situated within the last permafrost maximum extent), here termed "non-yedoma" thermokarst lakes, with a modern area of 126,685 km 2 (Table S2).We linearly reconstructed non-yedoma lake area in 100-year timesteps using the distribution of thermokarst lake initiation ages.We also assumed that the average age of existing thermokarst lakes in the dataset represents the typical "lifespan" of these lakes (8700 ± 5200 years).We applied an annual flux rate of 6.8 ± 4.4 g CH 4 m −2 yr −1 (Table S1) mineralized from predominantly younger permafrost and contemporary carbon substrates fueling methane production in non-yedoma environments (Fig. 6).This flux rate was temperature-adjusted and multiplied by reconstructed lake areas in each 100-year timestep to estimate non-yedoma thermokarst lake emissions over time.
Forward atmospheric Δ 14 CH 4 box model.A detailed description of the box model is given in the Supplementary Material of ref. 6 .In brief, atmospheric CH 4 burden is reconstructed from CH 4 mole fraction data in ice cores.The lifetime of atmospheric CH 4 is assumed to be 10.05 years.From 18 to 7.7 ka, CH 4 mole fractions from NEEM and WAIS ice cores 77 are used to calculate the global CH 4 source strength (Q tot ).From 7.7 ka onwards, we used CH 4 mole fraction compilation from Beck et al. 78 The isotope mass balance of 14 C-CH 4 in the atmosphere is given by 14 where 14 C atm is the atmospheric Δ 14 CH 4 , Q tot is the global CH 4 emissions, 14 C TKL is the 14 C activity of CH 4 emissions from thermokarst lakes (Q TKL , which is estimated in this study), and 14 C intcal is the 14 C activity of contemporaneous/ modern CH 4 emissions (Q mod ).All 14 C terms in Eq. 1 are in pMC (percent modern carbon) unit while the flux terms (Q) are in Tg CH 4 yr −1 .We assumed that at the time of emission, 14 C mod is equal to 14 C activity of 14 CO 2 derived from IntCal20 79 .
As discussed above, we ran two scenarios for 14 C TKL , one assuming 14 CH 4 activity of 3132 years and 7828 years old relative to contemporaneous emissions.This assumption is implemented using Eq. 3.
14 C TKL ðtÞ ¼ 14 where t TKL is either 3132 yr or 7828 yr and λ is the 14 C decay constant in yr −1 .
Climate forcing.We used equations published in the 2013 IPCC report 80 to calculate the radiative forcing and global temperature increase resulting from lake methane emissions estimated in this study.We assumed a global mean equilibrium temperature change following a doubling of atmospheric CO 2 of 3 °C.

Fig. 5
Fig.5Comparison of expected atmospheric Δ 14 CH 4 based on bottom-up thermokarst lake emission estimates A and B with published ice core Δ 14 CH 4 data.Ice core data and the IntCal20 curve representing the contemporary radiocarbon age of atmospheric CO 2 are presented on their respective timescales (black line).14CH 4 data from ice core studies are shown as blue6 -and aqua 47 -filled diamonds.Expected atmospheric Δ 14 CH 4 based on bottom-up thermokarst lake emission estimates for scenario A using 14 CH 4 from lakes with a mean age of 3132 yrs (dotted yellow line) and 7828 yrs (dotted pink line).Expected atmospheric Δ 14 CH 4 based on bottom-up thermokarst lake emission estimates for scenario B using 14 CH 4 from lakes with a mean age of 3132 yrs (solid brown line) and 7828 yrs (solid red line).Error bars represent the 95% confidence interval.

Table 1
15ke areas, the relative contribution of emission modes to total annual methane emissions, and methane emission14C ages (mean ± SD in yrs BP) compiled from published literature15.
A summary of emission modes by lake type are provided in TableS1.ND no data.