Decadal-scale hotspot methane ebullition within lakes following abrupt permafrost thaw

We mapped historic lake margins using high-resolution (1 m) aerial photos from 1949, 1967, and 1985 and 2009 lake margins from 2.5 m resolution Satellite Pour l'Observation de la Terre (SPOT) imagery to quantify lake area change and thermokarst shoreline expansion rates (supplementary information 1.2 and 1.3). To relate lake area changes to CH₄ fluxes, we classified lakes in which one or more shoreline had retreated as 'expanded lakes' and lakes that had formed entirely since 1949 as 'new lakes.'

We used a combination of field measurements and remote sensing to estimate CH₄ ebullition in lakes. In the field, we took advantage of seasonal ice cover on lakes to map discrete points of re-occurring ebullition (i.e. seeps) that have been shown to dominate annual CH₄ emissions in yedoma thermokarst lakes (64% to 72% of total lake CH₄ emissions in Siberia and Alaska; Walter et al 2006, Greene et al 2014, Sepulveda-Jauregui et al 2015). In winter, downward growing lake ice temporarily traps upward rising CH₄ bubbles, revealing the locations of ebullition seeps (figure 2(a)). Due to lower ebullition rates, bubbles of A-, B-, and C-type seeps are trapped in and beneath winter lake ice. Methane that does not dissolve out of the under-ice bubbles and get oxidized in the water column is subsequently released during spring ice-melt (Greene et al 2014). The strongest seeps, termed hotspots, release CH₄ to the atmosphere year-round (Walter et al 2006). Hotspots prevent ice-formation at the seep location due to water column convection caused by high bubbling rates (figure 2) (Zimov et al 2001, Greene et al 2014, Daanen 2020). Using a hand-held Garmin global positioning system (GPS) and a differential GPS (dGPS) we conducted ice-bubble surveys to map the locations of 11 696 A-, B-, C- and hotspot ebullition seeps within 92 survey transects on 29 lakes (table 1; supplementary information 1.4). Through analysis of high-resolution, early-winter aerial photos we mapped the locations of 496 open-hole hotspot seeps in 22 lakes (table 2; supplementary information 1.6). In aerial photos, dark, round hotspot holes in ice (usually <1 m diameter) had a strong contrast against white, snow-covered lake ice (figure 2(c), supplementary figure 2). On the ground, we mapped hotspots across whole-lake surfaces in 21 lakes, 14 of which were also aerially photographed. Seep bubble densities were converted to ebullition rates using a transfer function based on sim213 000 bubble-trap flux measurements associated with different seep types (Walter Anthony and Anthony 2013). Distinct ebullition events (n = 275) were sampled for analyses of gas composition, stable isotopes and radiocarbon dating (supplementary information 1.5).

Zoom In Zoom Out Reset image size

We compared spatiotemporal hotspot ebullition patterns and radiocarbon ages to taliks mapped with geophysical surveys (supplementary information 1.8) and lake sediment coring (supplementary information 1.9) at two intensive study lakes, Goldstream Lake (expanding lake) and Big Trail Lake (new lake) in a similar geologic setting.

Statistical methods are provided in supplementary information 1.10.

The number of lakes in our study area identified in high resolution imagery more than doubled from 1949 (132 lakes) to 2009 (278 lakes). Formation of 202 new lakes in yedoma permafrost deposits offset the disappearance of 71 existing lakes (figure 3(a)). Fifteen non-yedoma lakes formed between 1949 and 2009, offsetting the drainage of three existing lakes and doubling the total number of non-yedoma lakes in the study extent. Non-yedoma lakes comprised 6.8% of the 278 lakes in 2009.

Zoom In Zoom Out Reset image size

Total lake area growth accelerated from 1949 to 2009 (figure 3(b)). New lakes comprised 35% of total lake area in 2009, compensating for shrinkage and drainage of some existing lakes and contributing to a net lake area increase of 37.5% from 1949 (71.6 ha) to 2009 (98.5 ha).

Changes in lake number and area across the remote sensing time series (figure 4.) highlight the dynamics associated with recently developed thermokarst lakes in ice-rich permafrost in the study area. Lakes larger than 5 ha remained fairly stable between 1949 and 2009. Exceptional lake area expansion occurred in the 2.5–5 ha class between 1967 and 1985, and a large increase in lake number occurred in the 1.0–2.5 ha size class during the following time period (supplementary information 1.10). The observed lake dynamics highlight the role of recently formed and maturing lakes on the spatio-temporal patterns of CH₄ ebullition identified in the study area.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

The average expansion rate of shorelines in lakes that existed prior to 1949 was 0.27 ± 0.09 m yr⁻¹ (mean ± standard deviation) for yedoma lakes and 0.27 ± 0.08 m yr⁻¹ for non-yedoma lakes (1949–2009 period). Expansion rates of new yedoma lakes, which formed after 1949, were 0.63 ± 0.27 m yr⁻¹.

Ebullition fluxes estimated from ice-bubble surveys of A-, B-, C- and hotspot seeps coupled to bubble-trap fluxes revealed >2 times higher ebullition in new lakes compared to older, expanding lakes in both yedoma and non-yedoma permafrost soil types (table 1; Welch two-sample test, p < 0.001). Ebullition in yedoma lakes was 3.0-fold and 3.7-fold higher than in non-yedoma lakes for new lakes and expanding lakes, respectively (Two-sample test, p < 0.001).

In new yedoma thermokarst lakes, hotspot seeps, which were ⩽0.004% of lake area (based on bubble-cluster diameter of sim5 cm measured with a 1 cm grid at the lake surface; ten seeps measured on two lakes) comprised 5 ± 1% of total lake seep ebullition. In older expanding yedoma thermokarst lakes, hotspot seeps, which were ⩽0.005% of total lake area, were 7 ± 2% of total lake seep ebullition (table 1). Hotspot density was nearly two times higher in the new thermokarst areas compared to stable open water areas that were already lakes in 1949 (table 3). Only one of the five non-yedoma lakes had hotspots (table 1), and hotspot density in the expanded area of this lake was 4–6 times lower than in new yedoma lake areas (table 3).

Methane dominated ebullition bubbles (sim80% by volume), while bubble CO₂ concentrations were typically 1% to 2%. Bubble CH₄ concentration did not significantly vary by lake type, but it was higher in hotspot seep bubbles compared to non-hotspot-seep bubbles (ANOVA; p < 0.001) (table 4). Carbon stable isotope values of CH₄ (δ¹³C_CH4: −77‰ to −54‰, min to max) were consistent with microbial methanogenesis. The radiocarbon age of CH₄ (¹⁴C_CH4) in hotspot seeps was younger for new yedoma lakes (8526 ± 741 ¹⁴C years BP) compared to expanding yedoma lakes (20 079 ± 1227 ¹⁴C years BP) (table 4). ¹⁴C_CH4 in non-hotspot seeps of expanding yedoma lakes were younger (7473 ± 1762 ¹⁴C years BP) than from hotspots in these lakes. Non-hotspot seeps were also younger in new yedoma lakes (4742 ± 803 ¹⁴C years BP) compared to hotspots in new lakes. Non-yedoma lake hotspots had the youngest ¹⁴C_CH4 age (2450 ± 25 ¹⁴C years BP).

Spatial patterns of hotspot seeps derived from aerial photo analysis were highly correlated with those based on ground surveys (Kendall's rank correlation tau = 0.6277; p < 0.0001) (figure 6, supplementary figures 3–5). In expanding yedoma lakes, hotspots were most common near the shorelines that had changed from land to open water since 1949. We observed hotspots across much of the lake surfaces in new yedoma lakes. Very few hotspots were observed in non-yedoma lakes, regardless of lake expansion. Figure 7 shows photographs of expanding and stable interior Alaska lake shorelines.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Of the lake hotspot seeps identified through our remote sensing analysis, 86% (76% in field work, FW) occurred in areas of lakes that had changed from land to water between 1949 and 2009, with the remaining 14% (24% FW) occurring within 29 m (36 m FW) of the 1949 shoreline towards the center of the lake (figures 8(a) and (b)). Excluding Big Trail Lake, a large new lake with sim200 hotspots, we found hotspots were still predominately (78% and 60%) in newly formed lake areas based on aerial photo and field work measurements, respectively.

Zoom In Zoom Out Reset image size

We used geophysical measurements of talik depths and lake sediment coring to improve understanding of lake-surface ebullition patterns and CH₄ ages in relation to abrupt thaw in expanding versus new lake types. Geophysical measurements, including both ground-based electrical resistivity tomography (ERT) and airborne electromagnetics (AEM) revealed differences in talik depths between two intensively studied lakes, Goldstream Lake (expanding lake) and Big Trail Lake (new lake) (figure 9). In Goldstream Lake, where radiocarbon dated lake sediment cores suggest the center of the lake originated approximately 850 years ago (figure 10), talik depth extended through the sim20 m thickness of in-situ thawed yedoma (i.e. taberites), which had late-Pleistocene age-dates (>39 600 to >54 585 years BP) to 40 m as inferred by the region of low electrical resistivity (figures 9(c) and (e)). The talik was not as deep (<30 m) closer to the expanding eastern thermokarst margin (figure 9), where hotspot ebullition was concentrated and where bubble ¹⁴C_CH4 ages were up to 30 800 years. BP (figure 10(e)). In contrast, Big Trail Lake, formed from a wetland between 1949 and 1967, had a shallower talik (<15 m) (figures 9(d) and (f)) and younger ¹⁴C_CH4 in hotspots (5440–13 650 years BP) that extend across the whole lake surface (figure 10(f)).

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

The 37.5% net increase in lake area in our interior Alaska study area from 1949 to 2009 is large, but not unprecedented among other ice-rich, permafrost regions of Quebec (50% increase 1957–2003, Payette et al 2004) and Siberia (18% increase 2001–2009; Boike et al 2016). Mean lake expansion (0.27 m yr⁻¹) rates we observed tended to be lower than those observed on the windier Arctic coastal plain of Alaska (sim0.7 m yr⁻¹; Arp et al 2011, Lewellen 1970) and other pan-Arctic thermokarst lakes of the tundra (Jones et al 2011).

Our multi-temporal analysis demonstrates fluctuations in lake number and area across the study period that coincide with expected effects of climate on lake development. We did not see the formation of many large lakes that would otherwise point to non-climate related drivers. While we cannot rule out unknown lags in the impacts of road construction, agriculture, fire, and mining on lake development, these anthropogenic activities largely occurred in the first half of the 20th century, prior to our earliest record of lake areas (year 1949). If human activities were a primary driver of lake expansion, then we would not expect the large decrease in lake number and area that occurred between 1949 and 1967. Hereafter we discuss lake development in relation to climate. Non-climate related factors (e.g. humans, beavers) and wildfire are further discussed in supplementary information 2.

Long-term multi-decadal warming corresponded with lake area increase (figures 11(b) and (e)), but short-term pulse events may also have contributed. For instance, summer air temperatures warmed since 1949 (0.03 °C yr⁻¹) with each period showing a warmer average temperature, while the trend within each time period was inversely correlated to lake expansion. An increased frequency of imagery could provide an assessment of the role of extreme weather events on lake expansion and abrupt thaw. Literature on abrupt ice-wedge degradation in the continuous permafrost has emphasized the impact of one summer's extreme warmth (Jorgenson et al 2006, Liljedahl et al 2016), which may be equally important for lake abrupt thaw. Thirty-three of the warmest summers in the entire time series have occurred in the two most recent time periods (1967–2009) when lake area increase was most dramatic.

Zoom In Zoom Out Reset image size

Changes in autumn air temperature may also have affected lake area change. While p-values were too high to confirm statistical significance (supplementary table 1), we noticed corresponding trends in lake area growth and autumn (September–October) air temperature during the three study intervals (figure 11(d)). Lake area increased with warming fall air temperatures of 0.02 °C yr⁻¹ (1967–1985) and 0.06 °C yr⁻¹ (1985–2009), although the average fall temperature was relatively unchanged (1.9 ± 1.5 and 1.5 ± 2.2 °C, mean/standard deviation, respectively) between the two time periods (figure 11). Lake area decreased (−6%) when fall temperatures decreased (−0.04 °C yr⁻¹ 1949–1967), although expansion of floating mats in some lakes by 1967 also contributed to the shrinkage of open water lake area. Warmer autumn seasons suggest an increase in the rain-to-snow precipitation ratio (Bintanja and Andry 2017) during a time when evapotranspiration losses are relatively low (Nakai et al 2013). Any increased groundwater storage and subsequent groundwater flow may therefore promote thermal erosion of permafrost, ultimately favoring shoreline erosion of thermokarst lakes. Increased rainfall has stimulated increased active layer depths across a variety of interior Alaskan boreal ecosystems in recent years (Douglas et al 2020).

Lake area changes also followed precipitation patterns. A 6% decrease in regional lake area from 1949 to 1967 was concurrent with a persistent decrease in summer precipitation from 1930 to 1967. However, regional lake area subsequently increased 13% from 1967 to 1985 and then 30% from 1985 to 2009, while summer precipitation increased only during the last period (figure 11(c), supplementary table 1).

We found that lake and climate changes were coupled in the same ways for existing lakes versus new yedoma lakes and non-yedoma lakes from 1967 to 2009; however, increases in new lake area were most pronounced. When new ponds first form, their development may be slow when shallow water freezes to the bed in winter since bedfast ice impedes talik development (Soloviev 1973). However, warming during the summer and fall accelerates development. Once lakes become deep enough to retain floating ice in winter, talik growth continues year-round and at a faster rate annually (Burn and Smith 1990, Arp et al 2015). While, the number of new yedoma lakes increased by only one lake from 1985 to 2009, the area of the 201 previously formed, new yedoma lakes doubled during that period (figure 3(a)). It is likely that this rapid expansion was caused by warmer, wetter summers that accelerated degradation of ground ice surrounding the newly formed lakes.

We attribute the 2.6-fold higher ebullition in new lakes compared to older, expanding lakes (figure 5, table 1) to the fraction of lake surface actively expanding into yedoma permafrost deposits. Radial expansion of taliks into labile yedoma substrates beneath new lakes leads to high emissions across entire lake surfaces. In contrast, the labile permafrost carbon pool has been previously exhausted beneath more of the older expanding lakes, such that high ebullition is now concentrated in areas of expanding shorelines (Walter Anthony and Anthony 2013).

Our interpretation that CH₄ production and ebullition decrease over time as thawed permafrost soil organic carbon substrates are used up is also supported by the empirical relationship observed between lake age and ebullition among age-dated thermokarst lakes in the yedoma regions of north Siberia, the northern Seward Peninsula, Alaska, and interior Alaska (figure 12).

Zoom In Zoom Out Reset image size

Goldstream Lake and Big Trail Lake further support the concept that permafrost soil organic carbon fueling methanogenesis in taliks diminishes over time since thaw (figure 10). Because it was older and already had a deep talik, thaw beneath Goldstream Lake expanded radially at greater depths during the past 60 years than thaw beneath Big Trail Lake. The deeper radial thaw beneath older lakes with deeper taliks provides microbes access to older, more deeply buried syngenetic permafrost soil organic carbon. This results in older ¹⁴C ages in hotspot ebullition seeps from existing lakes compared to new lakes (figure 5(c), table 4), a pattern that we attribute to talik depth (figure 10).

In both lakes, the ¹⁴C age of CH₄ in hotspot seeps decreased with distance from the thermokarst margin (figures 10(e) and (f)). This is likely explained by mixing during bubble ascent through sediments. As bubbles rise, ¹⁴C-depleted CH₄ formed deep within taliks mixes with younger CH₄ formed in surface sediments, which are thicker in the lake center than along margins (Farquharson et al 2016). Conversely, non-hotspot ebullition seeps, formed closer to the sediment surface, get older with distance from the expanding thermokarst margins due to an increasing contribution of ¹⁴C-depleted soil organic carbon that is reworked from expanding lake margins and deposited as laminated lake sediments in lake centers (Walter Anthony et al 2014, Farquharson et al 2016).

Based on higher bubbling rates, higher CH₄:N₂ ratios and older ¹⁴C ages in hotspot ebullition bubbles compared to other ebullition, hotspot-type seeps are thought to originate from the deepest, organic-rich talik horizons beneath thermokarst lakes (Walter et al 2006, 2008, Walter Anthony and Anthony 2013). If this is true, then locations of hotspots across lake ice surfaces should indicate the spatial boundary of organic-rich permafrost thaw beneath lakes. Conversely, the absence of hotspots seeps in yedoma lakes would indicate complete vertical degradation of yedoma permafrost.

We do not attribute the paucity of CH₄ hotspot seeps in the oldest areas of thermokarst lakes to microbial anaerobic oxidation of methane (AOM), which occurs in thermokarst lakes (Martinez-Cruz et al 2017, 2018), but which appears not to significantly impact CH₄ in ebullition hotspots (Winkel et al 2019). Rather, in cases where vertical thaw exceeds the thickness of the organic-rich yedoma soils, it is the lack of substrate for methanogenesis that explains the paucity of hotspot seeps in lake centers. Since we observed ebullition hotspots strongly associated with abrupt-thaw areas of lakes, including actively expanding lake margins and newly formed lakes (figure 6, supplementary figures 3 and 4), but fewer hotspots in older, stable lake areas where yedoma permafrost was presumably degraded long ago and its labile carbon fraction microbially exhausted, we propose that hotspot seep occurrence can be used as an indicator of thawing labile, permafrost carbon within taliks. Based on this conceptual framework, we used the distribution of mapped hotspot seeps relative to the 1949 shoreline (figure 8), combined with our observations of thermokarst expansion rates, as a proxy for the temporal limit of labile organic carbon availability after thaw (supplementary information 1.7). Since sim80% of hotspots occurred in areas where land changed to lake since 1949, we conclude that (a) hotspot CH₄ release from lakes occurs largely within decades (60 years) of permafrost thaw, and (b) hotspot seep emissions cease altogether within 130 years following the transition of land to thermokarst lake.

Our observations of a large hotspot pulse of CH₄ emissions occurring within decades of land transition to thermokarst lakes is consistent with numerical modeling of permafrost carbon mineralization in northern Seward Peninsula yedoma thermokarst-lake taliks, whereby 75% of the labile organic carbon fraction is mineralized within the first century following thaw (Kessler et al 2012). It is also consistent with laboratory incubations of yedoma taberal sediments, which showed enhanced methanogenesis at the base of the downward expanding talik lasting up to 40 years (Heslop et al 2015, 2020).

Additional ¹⁴C-depleted CH₄ emissions observed in other, smaller, non-hotspot ebullition seep types within and beyond the spatial limits of hotspots in lakes (e.g. figures 10(a) and (c)) indicate that permafrost-soil derived CH₄ emissions persist beyond the decadal scale following thaw. It is likely, however, that the source of the ¹⁴C-depleted substrate to these non-hotspot ebullition seeps is not limited to in situ thaw beneath lakes, but is derived also from permafrost thaw and erosion along shorelines and within the watershed. Pleistocene-aged, organic-rich silt is a major component of the laminated, surface-lake sediment facies and is interbedded with terrestrial and aquatic carbon accumulations (Walter Anthony et al 2014, Farquharson et al 2016). It is in these organic-rich, yet ¹⁴C-depleted surface lake sediments that shoreline- and watershed-derived permafrost soil carbon is mineralized to CH₄ (Heslop et al 2015). Thus, ¹⁴C-depleted CH₄ ebullition from non-hotspot seeps in stabilized lake areas implies the contribution of shoreline and watershed permafrost soils deposition in lakes to the persistence of ¹⁴C-depleted CH₄ emissions well beyond the decadal time scale of hotspot emissions following lake formation.

Up to one million square kilometers of new, abrupt-thaw thermokarst lake areas are expected to form in the Arctic this century in response to RCP4.5 and RCP8.5 warming (Schneider von Deimling et al 2015, Walter Anthony et al 2018, Turetsky et al 2020). Our work provides input for yedoma-region modeling. We expect the lower level emissions from older, thermokarst-lake areas to persist as a background emissions floor while the newly formed thermokarst-lake areas and their higher CH₄ emissions will represent the 'new CH₄' signal and add to the emissions from the older lakes over a 60+ year time scale. Non-yedoma lakes have lower ebullition.

Regional differences in climate, hydrologic, and ground-ice controls on talik development rates should impact the duration of hotspot emissions following permafrost thaw beneath lakes. Higher spatial resolution data on ground-ice distributions in future work would particularly benefit our understanding of thermokarst and CH₄ response. We also acknowledge that regional differences in surficial geology imply variability in the duration of the hotspot CH₄ pulse. For instance, in ice-rich, thermokarst-susceptible portions of the vast non-yedoma permafrost region (Olefeldt et al 2016), organic-rich permafrost soils are thinner than yedoma deposits (Hugelius et al 2014, Strauss et al 2016). This implies faster thaw and release of permafrost carbon. Consideration of surface geology is also important for understanding the role of thermokarst lakes in IPCC scenarios and partitioning ¹⁴C-depleted sources of atmospheric CH₄ between natural versus anthropogenic sources (Dyonisius et al 2020, Hmiel et al 2020), since non-yedoma permafrost soil carbon and thermokarst-lake CH₄ ages are younger than those in yedoma soils (Walter Anthony et al 2016, Elder et al 2018, Matveev et al 2018). Differences in surficial geology could also give rise to differences in bubble gas composition, preferential modes of CH₄ release (diffusion/ebullition), and wintertime sources of water column dissolved gases (sediments vs ebullition) that are released in spring (Greene et al 2014, Matveev et al 2018, 2019, Elder et al 2018, 2020).

Landscape drainage capacity (i.e. how long a thermokarst lakes lasts before it drains) is a factor that is not well accounted for in models. The magnitude of future abrupt-thaw thermokarst CH₄ emissions may be moderated if landscape topography and hydrology cannot sustain a large increase in lake abundance (van Huissteden et al 2011). However, other processes affecting lakes will exacerbate the permafrost carbon feedback. For instance, longer ice-free seasons (Surdu et al 2014) will lead to a higher proportion of sediment-produced bubbles being directly released to the atmosphere without impediment by the seasonal lake ice sheet (Greene et al 2014). This reduces dissolution and oxidation of bubble CH₄ in the lake water column. Warmer winter temperatures and thicker snow depths, observed in many regions, result in thinner lake ice and thus contribute to new talik growth in previously shallow bedfast-ice lakes (Arp et al 2016). The same snow—pond ice—permafrost feedback may further enhance new lake formation. Nutrient enhancement of ecosystem productivity under warmer climate scenarios could also increase the supply of contemporary substrates to methanogenesis, increasing future lake CH₄ emissions (Wik et al 2016), particularly by ebullition (Delsontro et al 2016). Climate warming leading to warmer sediment temperatures can increase total annual emissions by enhancing microbial respiration and making mineralization of recalcitrant sediment organic carbon pools more energetically favorable (Heslop et al 2019). Finally, extreme weather events leading to ice-wedge degradation (Jorgenson et al 2006, Liljedahl et al 2016) likely enhance lateral lake expansion and should be taken into account, with higher temporal resolution remote sensing imagery, to investigate sub-decadal rates of lake change and answer the question, 'How abrupt really is abrupt thaw?'