The evolution of methane production rates from young to mature thermokarst lakes

Thermokarst lakes, formed by permafrost thaw in the Arctic, are hotspots for methane (CH₄) and carbon dioxide (CO₂) emissions and are expected to double permafrost carbon emissions by the end of the century. While the implications of ongoing permafrost thaw on CH₄ dynamics in these lakes have been modeled, here we provide empirical data on CH₄ production dynamics as lakes evolve from young recently formed lakes to older lakes that have been present for hundreds of years. Sediment cores (up to 4 m long) were collected from the centers and thermokarst margins of a new thermokarst lake (Big Trail Lake (BTL), < 70 years old) and from an older thermokarst lake (Goldstream Lake (GSL), ∼ 900 years old) from the same interior Alaskan watershed. The highest CH₄ production rates were observed in the uppermost sediments near the sediment–water interface at the thermokarst margins of both lakes, with a steep decrease with sediment depth into the talik. BTL exhibited elevated CH₄ production rates, correlated with higher carbon lability for thermal-induced reactions measured by Rock-Eval analyses, suggesting its potential use as a proxy for organic carbon breakdown by methanogenesis. In contrast, GSL displayed lower CH₄ production rates, likely due to a longer period of organic carbon degradation and reduced carbon lability. The integrated sediment-column CH₄ production rates were similar (around 7 to 10 mol m⁻² yr⁻¹), primarily due to the thinner talik at BTL. Our data support the predictions that the formation and expansion of thermokarst lakes over the next centuries will increase CH₄ production in newly thawed Yedoma permafrost sediments, while CH₄ production will decrease as taliks mature and labile organic carbon is used up.

1 Introduction

Permafrost covers one-quarter of the Northern Hemisphere (Obu, 2021; Zhang et al., 2008). Significant warming of the Arctic and subarctic regions, with temperature increases of 2 to 5 °C relative to pre-industrial levels (Post et al., 2019), exacerbates the thawing of this permafrost. In turn, Arctic soils, which are currently a store of carbon, become a source of carbon to the atmosphere by emitting greenhouse gases (Schuur et al., 2015). It is estimated that until the end of this century thaw, rapidly thawed permafrost will become an important source of greenhouse gases (Turetsky et al., 2020). Among these gases, methane (CH₄) is expected to be a dominant driver of the circumpolar permafrost carbon radiative, responsible for up to 70 % of this effect (Walter Anthony et al., 2018).

Methane release is expected to increase in rapidly thawing permafrost beneath thermokarst lakes, which are CH₄ “hotspots” on the landscape, and they are expected to double permafrost carbon emissions and increase associated radiative forcing effects by 130 % by the end of the century (Elder et al., 2021; Hugelius et al., 2014; Olefeldt et al., 2016). Thermokarst lakes are formed when ground ice melts and water accumulates in subsidence areas (Hopkins, 1949). As thaw continues beneath the lake, previously frozen organic carbon in the thawed permafrost soils, called taberites (Farquharson et al., 2016; Strauss et al., 2013), becomes available for microbial degradation, which produces CO₂ and CH₄ (Freitas et al., 2025; Heslop et al., 2015; Walter Anthony et al., 2020). In particular, the organic-rich permafrost Yedoma soils of Alaska and Siberia are noteworthy reservoirs of old ¹⁴C-depleted labile soil carbon that is quickly degraded into greenhouse gases upon thawing (Dutta et al., 2006; Estop-Aragonés et al., 2020; Knoblauch et al., 2018; Zimov et al., 1997).

Models have simulated the changes in methane production as thermokarst lakes evolve (Kessler et al., 2012). These models show that while older lakes may have lower surface production rates, their taliks can still contribute significantly to total methane flux. Over centuries, the talik deepens through the Yedoma and into the bedrock, while intra-talik organic carbon in the in situ thawed Yedoma sediments (taberites) gradually becomes refractory. At later stages, thawed permafrost organic carbon is no longer available for CH₄ production, and the lake is no longer a significant source of permafrost-derived CH₄ to the atmosphere, unless permafrost soil carbon is re-transported from lake margins or the watershed to surface lake sediments (Walter Anthony et al., 2014).

However, while the implications of ongoing permafrost thaw on CH₄ dynamics in these lakes have been modeled, few empirical data exist on changes in CH₄ production in relation to the stages of thermokarst lake evolution. The limited studies on CH₄ production rates in the talik of thermokarst lakes in the Arctic exhibit substantial variability, with fluctuations of 3 to 4 orders of magnitude observed across different environments and studies. For example, the top sediment layer of a non-Yedoma lake on the North Slope of Alaska shows CH₄ production levels of 1000 nmol cm⁻³ d⁻¹ (de Jong et al., 2018). CH₄ production rates range from 2 to 35 nmol cm⁻³ d⁻¹ in the Doughnut and Vault lakes in discontinuous Yedoma permafrost in Alaska (Martinez-Cruz et al., 2018), while values as high as 350 nmol cm⁻³ d⁻¹ have been reported in Vault Lake, with the lowest in Goldstream Lake (Sepulveda-Jauregui et al., 2015). Our recent study in this region on short sediment cores (up to 1 m depth) constrained CH₄ production rates in the upper sediments, similar to Martinez-Cruz et al. (2018), based on radiocarbon and CH₄ accumulation incubations (Pellerin et al., 2022). Freitas et al. (2025) showed – by using radiocarbon dating, sediment incubations and sediment facies classifications – that CH₄ production can also occur deep (∼ 20 m) beneath Yedoma thermokarst lakes in sand and gravel layers. Altogether, the data raised several questions about the evolution of CH₄ production rates throughout the talik, the role of CH₄ oxidation and the lability of organic carbon in thermokarst lake systems, warranting further investigation.

In this study, Yedoma thermokarst lake sediment cores were retrieved up to 4 m in depth to measure, calculate and compare CH₄ production rates along the talik between two distinct lake systems. Then, the CH₄ production rates, indicating the susceptibility of the organic carbon to microbial degradation, were compared to the total organic carbon (TOC) and its lability based on thermally induced breakdown. The thermal lability was deduced from a Rock-Eval analysis, which involves gradual heating under pyrolysis conditions, followed by combustion of the residual sample (Behar et al., 2001). With the gradual heating, the generation of hydrocarbon, CO and CO₂ are monitored. Following pyrolysis, the residual organic matter is lean in hydrogen and structural changes, making it more aromatic and refractory. The long cores enabled a direct comparison of CH₄ production rates, porewater geochemistry and organic carbon lability along geographically proximate yet geomorphologically distinct lakes – a young lake (< 70 years) and an old lake (∼ 900 years) (Fig. 1). These lakes also improved our understanding of how CH₄ production rates and carbon lability are expected to change over time, providing the first empirical data for comparison to model predictions.

2 Methods

2.1 Region of study

The Goldstream Valley watershed is located about 10 km northeast of Fairbanks, Alaska (Fig. 1). This area has a subarctic continental climate with an average annual temperature of −3.3 °C and annual precipitation of 280 mm (Douglas et al., 2020). The vegetation is primarily composed of boreal lowland species. The emergent lake and nearby thermokarst environments have been described by Elder et al. (2021). Recent studies (Hasson et al., 2022) show that the mean annual ground temperature of permafrost is −0.26 °C at a depth of 7.2 m. The Goldstream Valley features discontinuous ice-rich Yedoma-type permafrost, originating from Late Pleistocene loess deposits that were remobilized during the Holocene. These deposits form thick layers over relic Goldstream formations, altering the soil chemistry and ice volume distribution (Péwé, 1975). Between 1949 and 2009, the number of thermokarst lakes in the valley has doubled and their total area increased by about 40 % (Walter Anthony et al., 2020). These lakes are interconnected by a shifting watershed that feeds into the Tanana River, which is part of the Yukon River basin.

Figure 1 Study sites: the young (∼ 70 years old) Big Trail Lake (BTL) and the old (850–900 years old in its center) Goldstream Lake (GSL) in the Goldstream Valley (© Google Earth), including the locations of the collected cores in both lakes (edge (E), center (C)).

In this study, two markedly different thermokarst lakes were studied. Big Trail Lake (BTL) (64.9189° N, 147.8212° W, 609 m² in 2009) is an actively expanding lake formed from a wetland and possibly a migrating fluvial channel sometime between 1949 and 1967 (Walter Anthony et al., 2020). Extensive geophysical surveys at BTL showed massive ice (e.g., foliated ice wedges) starting roughly at 10 to 15 m below the irregular talik shape (Walter Anthony et al., 2020). Permafrost and valley hydrology were investigated on the valley scale (Emond et al., 2018). BTL is surrounded by valley-bottom creeks and streams; some supply water to the lake by draining the surrounding upland fens and the historical channel. These tributaries sometimes run dry and isolate the lake hydrologically at the surface. The outlet of the lake feeds into Eldorado Creek.

Goldstream Lake (GSL) (64.9156° N, 147.8495° E, 1278 m²), located about 1.5 km from BTL, is actively eroding into relic Yedoma permafrost with likely much less reworked Yedoma, due to its location at the base toe slope of Goldstream Valley and slightly elevated above the Goldstream Creek watershed at 196 m. Geophysical surveys estimated a talik thickness range of 30 to 40 m (Emond et al., 2018; Péwé, 1975). A previous study cored 20 m of Goldstream Lake sediments and found ice-free sand and gravels beneath 16 m of thawed silt (Walter Anthony et al., 2020). Radiocarbon dates suggest that the oldest part of the basin is around 850 to 900 years old; however, rapid expansion of the eastern margin into Yedoma permafrost has occurred since 1949.

2.2 Sampling and profiles

We collected sediment cores using a vibro-corer deployed directly from the lake ice in March 2022. In each lake, one core was collected from the center of the lake and from the edge within the water-saturated zone with an aluminum liner. It was winched out of the sediment, tilted on its side and cut into ∼ 1 m sections (marked in Table S1 in the Supplement). These were quickly capped and transported by a snow machine and vehicle to the University of Alaska Fairbanks for further analysis. In the lab, each section was cut lengthwise, and one half was sampled horizontally at intervals of 15 to 25 cm. First, a subsample of sediment was taken immediately after slicing the core liner for CH₄ concentrations, and the stable carbon isotope composition (δ¹³C${}_{{\mathrm{CH}}_{\mathrm{4}}}$ and δ¹³C${}_{{\mathrm{CO}}_{\mathrm{2}}}$), about 3 mL of sediment, was taken in a cut-off 3 mL syringe and inserted directly into 20 mL vials fully saturated with an anoxic 5 M sodium chloride solution. Another cut-out 3 mL syringe was inserted into the sediment at each depth to extract sediment for incubation. For density measurements, samples were taken into 10 mL vials, weighed and then re-weighed after drying an aliquot of sediment at 60 °C for 4 d. Total organic carbon (TOC) concentration, the carbon stable isotope composition (δ¹³C_TOC) and Rock-Eval analyses were performed on a sediment subsample was freeze-dried. Porewater was extracted by Rhizons (Dickens et al., 2007). The Rhizons were inserted into the sediment, and a vacuum was created with a 10 mL syringe, which accumulated the extracted porewater after filtration through 0.22 µm. The filtrated porewater was then stored at 4 °C in 2 mL amber glass vials without headspace for the measurement of dissolved inorganic carbon (DIC) concentrations and its stable carbon isotope composition (δ¹³C_DIC) within a week from sampling.

2.3 Production rate experiments

A 3 mL cut-out syringe was inserted into the sediment in each core at 20 to 30 cm intervals to retrieve about 2 mL of sediment. The sediment was added to a 20 mL serum bottle, which was sealed with a butyl rubber stopper and crimped with an aluminum cap. The bottles were immediately vigorously shaken and purged with 99.999 % N₂ gas for 15 min to remove oxygen and other gases. Three serum bottles were taken at each depth. After weighing, the samples were stored in the dark at 4 °C. The increase in CH₄ concentrations in the headspace was recorded after 80, 130 and 160 d by GC-FID (see below), allowing the back-calculation of CH₄ production rates.

Methane production rates were measured at each time point (80, 130 and 160 d) using a small (100 µL) aliquot of gas from the headspace of the serum bottles. The headspace CH₄ concentration was converted to total CH₄ in the bottle based on the concentration, the volume of sediment and the volume of headspace in each serum bottle. The CH₄ production rate was then taken as the increase in CH₄ concentration over time. Since each depth had three serum bottles, the average CH₄ production rate was reported for each depth and the uncertainty of the CH₄ production rate was reported as the standard deviation of the mean. After 160 d, δ¹³C${}_{{\mathrm{CH}}_{\mathrm{4}}}$ was also measured. Since any CH₄ dissolved in the porewater at the time of sampling had been removed during the N₂ purge, the δ¹³C${}_{{\mathrm{CH}}_{\mathrm{4}}}$ values represented the newly accumulated CH₄ during the incubation.

Total profile CH₄ production rates, reported in mol m⁻² yr⁻¹, were calculated using the CH₄ production rates obtained from the incubation experiments and depth integration throughout the thawed talik. The inferred talik thickness was based on geophysical surveys and measurements of the taberite depth (Freitas et al., 2025; Walter Anthony et al., 2020).

2.4 Analytical methods

The headspace of the serum bottles of the CH₄ rate incubations was measured for CH₄ concentrations at each time point using the procedure described by Pellerin et al. (2022). In general, 250 µL aliquot of gas from the headspace was inserted into a gas chromatograph (GC, Thermo) equipped with a flame ionization detector (FID) and a packed ShinCarbon ST column (Restek) at 120 °C. The GC-FID was calibrated with a standard curve prepared with a CH₄ concentration of 1 %. CH₄ concentrations in the profiles were measured by injecting 3 mL of 99.999 % N₂, while simultaneously removing 3 mL of the saturated 5 M sodium chloride solution. After equilibrating for 2 weeks upside down, bottles were vigorously shaken and analyzed for CH₄, as described above. This method had a precision of ±2 µM. The δ¹³C${}_{{\mathrm{CH}}_{\mathrm{4}}}$ and δ¹³C${}_{{\mathrm{CO}}_{\mathrm{2}}}$ values were measured together by PreCon and GasBench II, interfaced with a Delta-V gas source isotope ratio mass spectrometer (GS-IRMS, Thermo). Internal standards of Tiso-2 (isometric instruments, 2.5 % CH₄) δ¹³C −38.3 ‰ and in-house (99.8 % CH₄) δ¹³C −48.6 ‰ were used, with an analytical error of 0.2 ‰. For δ¹³C${}_{{\mathrm{CO}}_{\mathrm{2}}}$ determination, a 99.8 % CO₂ lecture bottle (Scott Mini-Mix) with δ¹³C of −25.5 ‰ was used, and the results were reported with an analytical error of 0.3 ‰. All δ¹³C results were reported compared to Vienna Pee Dee Belemnite (VPDB).

The TOC in the freeze-dried sediment was measured after the removal of carbonates by the addition of 1 % H₃PO₄ and drying at 40 °C in silver cups until the sample stopped reacting with the acid. Triplicate samples and internal standards were then packed in tin cups and measured on an elemental analyzer and HS2022 IRMS (Sercon). TOC concentrations were measured using concentration calibration and peak heights, and δ¹³CTOC was reported vs. VPDB with the reference materials IAEA-600 (δ¹³C −27.7 ‰ VPDB), USGS62 (δ¹³C −14.8 ‰ VPDB) and USGS63 (δ¹³C −1.2 ‰ VPDB). The precision was ±0.3 ‰.

Organic carbon lability for thermally induced chemical reactions was characterized by a Rock-Eval analysis (Rock-Eval 6, Vinci Technologies). The technique determines the proportion of pyrolyzable C (PC) and residual carbon (RC). PC is composed of the sum of three pyrolysates: S1, composed mostly of small volatile molecules; S2, larger hydrocarbon molecules thermally cracking like algal cell walls; and S3, derived from oxygen-containing molecules. The residual carbon (RC) is released from the sample during the combustion cycle (Carrie et al., 2012). The indices used as proxies to the organic lability are (1) hydrogen index, which is calculated as HI = S2/TOC × 100 (Behar et al., 2001). Higher HI values indicate greater hydrogen-rich organic compounds, implying lability. (2) The ratio between PC and RC – a higher ratio indicates that the OM is richer in hydrogen and is more aliphatic; hence, it might be more available to microbial respiration. About 20 mg of the prepared samples were placed in the RE6, which was then ramped at a predetermined rate (25 °C min⁻¹) from 200 to 650 °C in the pyrolysis oven. The oven was cooled down, and the sample was transferred to the combustion oven, where it was ramped from 200 °C up to 850 °C at a rate of 25 °C min⁻¹.

The DIC and δ¹³C_DIC were analyzed after acidification of the porewater sample to convert DIC to CO₂. The δ¹³C_DIC was measured with GS-IRMS (Thermo) interfaced to GasBench II. Values were reported relative to VPDB with a precision of ±0.1‰. DIC concentrations were measured by integrating the signal of the sample on the IRMS. The signal was calibrated using peak heights NaHCO₃ solutions, prepared at concentrations of 3 to 10 mM. The precision was ±0.2 mM.

The comparison of δ¹³C${}_{{\mathrm{CH}}_{\mathrm{4}}}$ between the CH₄ in the profiles and the CH₄ produced in the rate incubations was used to quantify the contributions of CH₄ in the shallow sediments relative to the transport of CH₄ migrating or diffusing from deeper in the talik. This is assuming that the CH₄ measured in the profiles represents CH₄ that is produced in situ, plus CH₄ that is produced in greater depths and migrates upward, while that in the incubations (after purging) represents only the in situ CH₄ production. The bubble composition was assumed to be the most negative δ¹³CCH₄ value measured in the profiles because it migrates from deeper in the sediment, as was shown previously (Pellerin et al., 2022).

Equation (1):

$$\begin{array}{}\text{(1)}& \begin{array}{rl}& \mathrm{new}{\mathrm{CH}}_{\mathrm{4}}\mathrm{production}\mathrm{fraction}=\\ & {\displaystyle \frac{{\mathit{\delta }}^{\mathrm{13}}{\mathrm{CH}}_{\mathrm{4}}\left(\mathrm{PW}\right)-{\mathit{\delta }}^{\mathrm{13}}{\mathrm{CH}}_{\mathrm{4}}\left(\mathrm{bubble}\right)}{{\mathit{\delta }}^{\mathrm{13}}{\mathrm{CH}}_{\mathrm{4}}\left(\mathrm{incubation}\right)-{\mathit{\delta }}^{\mathrm{13}}{\mathrm{CH}}_{\mathrm{4}}\left(\mathrm{bubble}\right)}}.\end{array}\end{array}$$

3 Results

3.1 Methane production characteristics from porewater profiles

Microbial respiration was intensive in general, as indicated by the high DIC concentrations in all sites. In the center of BTL, DIC concentrations increase linearly from the surface towards 150 cm depth, with a peak of 37 mM at 200 cm, followed by a subsequent decrease back to 25 mM (Fig. 2a). On the other hand, at the edge of BTL, DIC concentrations ranged from 10 to 15 mM (Fig. 2b). Along the talik in the center of the lake, the δ¹³C_DIC values increased from 3 ‰ to 16 ‰ (Fig. 2a). On the edge of the lake, the δ¹³C_DIC values increased in the upper 50 cm from 5 ‰ to 17 ‰ and then remained similar (Fig. 2b).

Figure 2 Pore water profiles of GSL (circles) and BTL (diamonds) cores. Panels (A) and (B) correspond to the BTL center (full) and the edge (empty), respectively, while panels (C) and (D) correspond to the GSL center and the edge, respectively. Uncertainty on the measurement is smaller than the symbol, unless indicated.

The BTL center exhibited CH₄ concentration peaks of 70 and 300 cm (about 1 mM), while CH₄ concentrations in the remaining core sections ranged from 0.3 to 0.5 mM. The δ¹³C${}_{{\mathrm{CH}}_{\mathrm{4}}}$ values ranged from −55.6 ‰ to −69 ‰ in the deeper section (Fig. 2a). CH₄ concentrations at the BTL edge (Fig. 2b) exhibited values of 0.2 to 0.5 mM at depths of 0 to 150 and 200 to 300 cm, with peak values of up to 1.4 mM observed at depths of 150 and 350 cm. The shallow part exhibited δ¹³C${}_{{\mathrm{CH}}_{\mathrm{4}}}$ values of −55 ‰, which became more negative, with −71.5 ‰ in the deeper section.

At the GSL center, DIC concentrations remained at around 20 mM, with a maximum of 32 mM observed at a depth of 100 to 150 cm (Fig. 2c). The δ¹³C_DIC values decreased in the upper 100 cm, from −2 ‰ to −15 ‰, followed by a slight increase to −12 ‰ (Fig. 2c). At the edge of GSL, DIC concentrations exhibited a linear increase with depth, from 20 to 46 mM (Fig. 2d). The δ¹³C_DIC values (Fig. 2d) exhibited different trends to BTL, with negative values, decreasing from approximately −2 ‰ to −6 ‰ within the first 100 cm. Below this depth, there was a slight increase in δ¹³C_DIC values to −3 ‰.

Methane concentrations in GSL remained relatively steady at the center (1 mM), with a slight decrease at 150 cm to 0.7 mM (Fig. 2c). The CH₄ concentrations edge exhibited a varied range of 0.5 to 1.5 mM (Fig. 2d). The highest CH₄ concentration was observed in the deepest and shallowest parts of the core (close to 1.5 mM). The lowest concentration was found at a depth of 150 cm. The δ¹³C${}_{{\mathrm{CH}}_{\mathrm{4}}}$ values at the center of GSL ranged from −68.4 ‰ in the shallow part to −78 ‰ in the deeper part (Fig. 2c). At the edge, they ranged from −65 ‰ in the shallow part to −73 ‰ in the deepest section (Fig. 2d).

The surficial sediments of the lakes' edges were organic rich (TOC ∼ 10 %), with relatively low δ¹³CTOC (Table S1 in the Supplement) and low CH₄ and DIC concentrations. The TOC content decreased significantly with depth on the edge, while the center of the lakes had low TOC content. In BTL in both the edge and center, the δ¹³CDIC increased significantly with depth, with values typical of methanogenesis. In GSL, on the other hand, there was a significant decrease in δ¹³C_DIC with depth but relatively constant DIC and CH₄ concentrations and isotopes in the rate incubations, which suggest a small role of CH₄-related processes with some signature of methane oxidation in the upper sediments.

3.2 Methane production rates in sediment incubations

Methane production rates were measured by sediment incubation batch experiments. The sediment from BTL and GSL had the highest CH₄ production rates near the sediment–water interface. In BTL, the rates within the upper meter ranged from 2 to 20 nmol cm⁻³ d⁻¹ in the core taken from the center of the lake (Fig. 3a) and 7 to 35 nmol cm⁻³ d⁻¹ in the core taken from the edge of the lake (Fig. 3b). In the deeper sediments, the rates decreased to about 1 nmol cm⁻³ d⁻¹ in both sites.

Figure 3 CH₄ production rates (nmol CH₄ cm⁻³ d⁻¹) from incubation experiments of BTL sediments: center (A) and edge (B) cores. A steep decrease in rates with depth is observed in both the edge and center sites. Depth to permafrost beneath BTL is thought to be 10 to 15 m (Walter Anthony et al., 2020). CH₄ production rates (nmol CH₄ cm⁻³ d⁻¹) in GSL from incubation experiments of the center (C) and edge (D) cores. A noticeable decrease in rates with depth is observed at both the center and edge of the lake. The talik depth in GSL is estimated in the center to be between 15 and 40 m (see below).

At the center of GSL (Fig. 3c), the CH₄ production rates were 4 to 6 nmol cm⁻³ d⁻¹ in the upper 1 m of sediment. Deeper, the rates decreased, ranging from 2 to 4 nmol cm⁻³ d⁻¹. At the edge, higher rates were measured in the upper 50 cm (yet lower than in BTL), followed by a sharp decrease below 50 cm depth (Fig. 3d).

3.3 Total thawed talik CH₄ production

The total profile-integrated CH₄ production rate throughout the thawed talik indicates the potential flux of CH₄ out of the sediment into the lake water column. This accumulated rate depends on CH₄ production rates at individual depths, facies thicknesses (surface sediments, taberites, etc.) and talik thickness. Since we did not have samples below 4 m depth, we extrapolated our data to the depth range of the known taberite thickness, which was 10 to 15 m at BTL and about 16 m at GSL (Walter Anthony et al., 2020). Different fits were evaluated for calculating the total CH₄ production (Tables S3–S5), and a power law regression was chosen as the best fit for extrapolating the total CH₄ production rate (Table S4, Fig. 4).

Figure 4 CH₄ production rates as measured in the different cores (symbols) and the power law extrapolation (curves, equations in Table S5) down to 4 and 15 m (upper-right side). A constant low rate of less than 1 nmol cm⁻³ d⁻¹ below 3 m down to 15 m is also marked in the upper-right side, as well as the potential role of the high production of CH₄ near of the base of the young talik.

The total 12 m column CH₄ production calculated for the BTL edge and center cores was about 8.5 (±12 %) and 7.4 (±37 %) mol m⁻² yr⁻¹, respectively. In GSL, talik thickness is greater (< 40 m); however, the volumes of thawed silt are potentially comparable at both lakes, and the main difference is in time since thaw. Because BTL is a younger lake, we assume the talik sediments have thawed within the last 70 years, whereas talik sediments beneath the center of GSL are thought to have been thawed for 8 to 9 centuries (Walter Anthony et al., 2020). At GSL, our 15 m profile-integrated CH₄ production rates from the center and the edge were 7 to 11 mol m⁻² yr⁻¹ (Table S4). Large uncertainties stem from the extrapolation (Table S5), the exceedingly long tail followed by the real measurements (Fig. 4), the potential variability in the composition and thickness of the surface organic-rich sediments and the actual talik depth across both lakes. We also assume that the rates are extremely low down to the base of the talik, whereas there is a potential for high production rates near the base of the young taliks. These high rates together with the high surface rates in the edges of the lakes and the center of BTL can explain the ebullition there and the lack of ebullition at the center of GSL (Walter Anthony et al., 2020).

3.4 The source of CH₄

The source of CH₄ in the sediment was estimated by comparing the δ¹³CCH₄ values in the profiles to those measured in the rate incubations (after purging and waiting several months for in situ production) (Fig. S1). The δ¹³C${}_{{\mathrm{CH}}_{\mathrm{4}}}$ values in the profiles from BTL were about −55 ‰ to −75 ‰ (edge) and −60 ‰ to −70 ‰ (center), whereas the values from the incubations of sediment from BTL were about 10 ‰ to 20 ‰ all along the profiles. The same observation is made in GSL, where the δ¹³C${}_{{\mathrm{CH}}_{\mathrm{4}}}$ values of the dissolved CH₄ in the profiles were about −60 ‰ to −70 ‰ (edge) and −70 ‰ to −80 ‰ (center), with more positive values for incubations. The significant difference between the incubations and in situ profiles points towards an additional deep source of ¹³C-depleted methane in the natural environment, consistent with talik-sourced ebullition observations (Walter Anthony et al., 2020).

3.5 Organic carbon characteristics in the sediments

The TOC profiles in the center of the BTL sediments exhibited a gradual decline from 2 % to less than 1 % and in GSL from approximately 1 % (Fig. S2). At the edge of the lakes, a higher TOC was observed near the top of the cores, up to 13 % (BTL) and 9 % (GSL), gradually decreasing with depth. The δ¹³CTOC in the center was constant near −27 ‰ in the upper 150 cm (Table S1). An increase to −25‰ was observed around 200 and 350 cm. At the edge, δ¹³CTOC values were around −28 ‰ at the upper part, with an increase observed below 200 cm to −25 ‰.

Rock-Eval data obtained at both the center and edge sites of both lakes exhibit a correlation between organic index values and sediment depth. The HI of BTL was mostly much higher than that of GSL and is indicative of a mix of Type I–III kerogen for the young lake (BTL) and Type III kerogen for the mature lake (GSL). The OI, on the other hand, was extremely high for both lakes, exceeding 150 mgCO₂ g⁻¹ TOC. In both the edge and center cores of BTL, there was an inverse dependency between the HI and OI, suggesting that with the loss of H, the OM became more oxidized. This relationship is missing from the mature lake, in which the OM lost most of its H, presumably due to microbial degradation (Fig. S3). The $\mathrm{PC}/\mathrm{RC}$ ratio in both sites of BTL decreased with depth and stabilized at 250 cm, with the center core generally exhibiting lower ratios compared to the edge. The $\mathrm{PC}/\mathrm{RC}$ ratio at the edge of GSL slightly decreased with depth. When comparing the two cores, the lower values were observed at the GSL center core, meaning it is more refractory than the edge (Fig. 5).

Figure 5 Sediment profiles of the ratio of pyrolyzable carbon to refractory carbon ($\mathrm{PC}/\mathrm{RC}$). The ratio decreases not only with depth but also with the evolution of the lakes, as organic carbon becomes more refractory. The uncertainty of the data is smaller than symbol sizes (circles – BTL, diamonds – GSL, full symbols – center, empty symbols – edge of lakes).

4 Discussion

4.1 Total CH₄ production in the talik

The CH₄ production rates observed in the upper sediments of this study are similar in magnitude to the ones observed in our previous studies (Lotem et al., 2023; Pellerin et al., 2022). Previous studies in which cores depth was limited to about 1 m depth were unable to provide a full understanding of the characteristics of CH₄ dynamics within thermokarst lake taliks. In Vault Lake, another thermokarst lake in central Alaska, a decrease in CH₄ production was observed, with depth down to 6 m, but CH₄ production rates were significantly higher than most reports (Heslop et al., 2015). Freitas et al. (2025) also showed low but significant cumulative anaerobic respiration throughout the taberal sediments, down to the gravel zone at around 16 m in GSL, with scattered values below.

Since we were able to sample up to 4 m in the talik of both BTL and GSL and in both the center and the edge (near the lakeshore) locations, our rate measurements provide an important confirmation of conceptual and numerical model predictions. They show high CH₄ production rates in the surface with a steep decrease with depth in the talik as the taberite organic carbon becomes more refractory over time. Our findings are also consistent with the findings from an 8 m permafrost sediment core from the Lena delta, where the highest CH₄ production rates were observed in the first 125 cm (Wagner et al., 2007).

The study also serves as the first empirical test of models predicting changes in CH₄ production rates with the evolution of lakes. It shows that CH₄ production rates are highest in the top sediments and are low but still significant in the deep talik. We did not measure significant differences between CH₄ production rates at depths deeper than 100 cm at BTL and GSL (Fig. 4). However, other lines of evidence such as the different isotopic enrichments of the DIC pools strongly suggest that there should be differences in CH₄ production rates between BTL and GSL. It may be that the low rates and higher uncertainty of the measurements in the deep talik made it impossible to differentiate between the rates in the two lakes.

The long cores and the steep CH₄ production decrease with depth enable the estimation of the total production rates with higher certainty than the accumulated rates calculated in the 1 m cores and the estimated constant rates along the talik (Pellerin et al., 2022). It should be noted that despite reaching down to 4 m in the talik, we still needed to extrapolate the rates to the entire thickness of thawed taberal sediments, which has a depth of 10 to 16 m (Walter Anthony et al., 2020). It has been suggested previously that CH₄ production rates may increase with the depth of the talik because more recently thawed permafrost might release more labile organic carbon available for degradation (Walter Anthony et al., 2014) due to a rapid turnover time upon thaw (Schädel et al., 2014; Shaver et al., 2006). However, previous studies have not found a significant rise in CH₄ production rates with depth, except near the thaw front at the base of the young talik (Heslop et al., 2015). This potential rise in young talik is marked in Fig. 4 and may explain the lowest accumulated production rate in the center of GSL. The recent study of Freitas et al. (2025) showed similarly low values in the taberites overlying the sand/gravel zone in GSL at about 16 m and then scattered high values (in part because they were normalized to the total organics that was very low). To extrapolate our rates deeper into the talik, we experimented with several approaches such as using power law decay or constant rate, which reflect the lower rates obtained in the incubations. Using the integrated rates with a power law equation, we calculated total talik CH₄ production (Table S4). We decided to present a power law fit since it indicates a declining rate in the deeper talik and reflects the overall trend observed in our incubations (Table S5). However, it is important to highlight the fact that all models that we experimented with yield similar conclusions regarding the total talik CH₄ production. The model used to extrapolate CH₄ production rates deep in the talik does not affect our conclusions.

An important finding is that facies thickness and talik depth play a significant role in determining total talik CH₄ production, instead of just lake age and location in a lake. This is because, despite measurably higher CH₄ production rates in the upper meter of sediment cores, as well as near the edges of the lakes, the low but relatively constant CH₄ production rates observed at depths in all the sites (5 % to 10 %) impact the integrated talik CH₄ production rates (Fig. 4).

The new total thawed talik CH₄ production rates are of the same magnitude as the CH₄ emissions measured by previous studies in BTL (Elder et al., 2021) and other discontinuous permafrost thermokarst lakes (Fig. 6). This similarity supports only a minor role for CH₄ oxidation (aerobically and anaerobically) in the lakes. It also fits the finding in the upper 1 m of several thermokarst lakes, which shows that the anaerobic oxidation of CH₄ rates, as deduced from batch experiments, are 2 orders of magnitude lower than CH₄ production and not a significant sink of CH₄ (Lotem et al., 2023).

Figure 6 Total talik CH₄ production in BTL and GSL (mol m⁻² yr⁻¹), obtained by integrating the CH₄ production rates over the potential range of thawed talik depths in the sediments (from the surface to the marked depth) using a power law extrapolation. The middle of the box corresponds to the mean estimate, and the uncertainty is the length of the boxes.

The implication of this observation is that as lakes mature, total thawed talik CH₄ production rates can remain similar or even increase, although the carbon becomes less available for microbial degradation. This is when the talik deepens fast enough and offsets the overall drop in carbon lability. This is the case observed for both BTL and the eastern thermokarst margin of GSL. The question is whether there are control and a link between the organic carbon, its nature and CH₄ production rates in BTL and GSL, as discussed below.

4.2 Factors influencing CH₄ production rates and potential net emissions

We further investigated the controls of CH₄ production in the BTL and GSL taliks. These lakes are in the same valley system, and it seems likely that they have the same sources and that CH₄ production will be influenced from in situ differences. As the isotopic composition is much more sensitive for quantifying sedimentary processes rather than the concentration changes (e.g., Sivan et al., 2011), δ¹³CDIC and δ¹³CCH₄ can be used to determine the controls on CH₄ production in the lakes. The first clear observation is that methanogenesis is strongly imprinted in the δ¹³CDIC of the young recently thawed talik of BTL and not in older GSL talik, which has been thawed longer (Fig. 2). The δ¹³C_DIC values of BTL are positive (around 10 ‰ to 20 ‰) as methanogenesis produces ¹³C-depleted CH₄ and, by mass balance, ¹³C-enriched DIC. The δ¹³C_DIC values in the sediments of GSL are about −5 ‰ to −10 ‰ at depth (Fig. 2).

The δ¹³C${}_{{\mathrm{CH}}_{\mathrm{4}}}$ values of both lakes of −78 ‰ to −53 ‰ (Fig. 2) support this picture. They are typical of freshwater sediments, which are usually dominated by microbial acetoclastic methanogenesis (e.g., Whiticar et al., 1986). Methanogenesis pathways and their importance in this system were quantified by radiocarbon probing (Pellerin et al., 2022) and by the δ¹³CCH₄ vs. δDCH₄ (Liu et al., 2025). The relatively heavy δ¹³C${}_{{\mathrm{CH}}_{\mathrm{4}}}$ values suggest low CH₄ oxidation (e.g., Whiticar, 1999; Sivan et al., 2011), which results in ¹³C enrichment in CH₄. The control of methanogenesis on BTL profiles is evident from looking at the δ¹³C${}_{{\mathrm{CO}}_{\mathrm{2}}}$ values and plotting them against δ¹³C${}_{{\mathrm{CH}}_{\mathrm{4}}}$, similar to Krause and Treude (2021; Fig. S7). The data support acetoclastic methanogenesis and CH₄ oxidation in all sites.

The main difference between the lakes is likely related to the fundamental difference in age of the two taliks, which influence the lability of organic carbon and the resulting biogeochemical cycles. The lability difference of the organics was further investigated by assessing the properties of the organic carbon in the talik of BTL and GSL.

Methane production rates are highly correlated to TOC and indices of thermal-induced lability in the upper sediments (Fig. 7a), but methane production rates and TOC are not correlated in the deep talik of BTL (Figs. S4–S5). The correlation is between two independent methods that estimate the lability of organic carbon to microbial degradation – the lability based on thermal degradation and the CH₄ production rate. The most labile fraction of pyrolyzed carbon is correlated with the CH₄ production rate in both the center and edge sites (Fig. 7b–d). The results demonstrate that the highest CH₄ production rates observed in our study are likely due to both the lability and quantity of organic matter.

Figure 7 Methane production rates (including its uncertainty) as measured in the top upper 100 cm of both the center and edge cores compared with the TOC content (A), the hydrogen index (B) (R² = 0.67) and the most labile organic carbon compounds (Figs. S1, S2, (C, D), as measured using Rock-Eval). A positive relationship (R² = 0.67 to 0.98) is evident in the upper 100 cm for BTL (circles) but not for GSL (diamonds).

The pyrolyzed carbon to residual carbon ratio ($\mathrm{PC}/\mathrm{RC}$) is a simple measure of how labile the organic carbon is. As the $\mathrm{PC}/\mathrm{RC}$ ratio is higher, the sample has more hydrogen and has a more pyrolyzable fraction (Carrie et al., 2012). Throughout the four cores taken for this study, the $\mathrm{RC}/\mathrm{PC}$ ratio was strongly dependent on depth in the talik (Fig. 5), which is consistent with older less labile organic matter deeper in the sediment because of the depositional history. However, the most interesting observation was the decreasing trend of lability to thermal-induced reactions between sites. The BTL edge had the highest $\mathrm{PC}/\mathrm{RC}$ values, which correspond to organic carbon with the highest lability, followed by the BTL center and the GSL edge, with the lowest value observed at the GSL center (Fig. 5). The lake edges exhibit the highest $\mathrm{PC}/\mathrm{RC}$ ratios compared to the center of lakes, suggesting that the center of the lakes contain more refractory organic carbon, which may be a result of the edges of the lakes being “younger” (the age since thaw). Alternatively, it could be that the edges have additional input of organic carbon from vegetation or runoff from land, as supported by the labile fractions found in the edges (Fig. 7c). Both the edge and center of GSL, on average, had lower $\mathrm{PC}/\mathrm{RC}$ ratios than both the center and edge of BTL. This is interpreted as GSL containing a greater fraction of less labile organic carbon, consistent with the longer time since the permafrost thawed and formed the talik and lake at this site. Carbon thus becomes less labile for microbial degradation as the lake evolves. In permafrost environments, organic carbon lability is not always correlated directly with what is termed the “age” of the organic carbon, such as in most marine sediments where lability, depth and age are often directly correlated. In the talik of thermokarst lakes, time since thaw seems to have some control on the lability of the organic carbon and exerts control on the CH₄ production rate. Altogether, our study emphasizes the potential use of the lability of thermal-induced reactions as a proxy for organic lability for methanogenesis.

It therefore seems that at the onset of thawing, during the formation of a young thermokarst lake, a high amount of highly labile organic carbon accumulates and begins to degrade. As thawing progresses, the lake margins expand and deepen and the refractory organics remain, leading to lower CH₄ production rates. A proposed conceptual mechanism is presented in Fig. 8. This mechanism can be tested in future studies, which include expanding the study to additional lakes, incorporating seasonal sampling and integrating microbial community analyses to better understand CH₄ cycling processes.

Figure 8 Proposed conceptual mechanism of OM lability and CH₄ emission as thermokarst lakes mature. Laboratory incubations support the conceptual and numerical models of decreasing carbon lability leading to lower CH₄ production as Yedoma thermokarst lakes mature. The gray sediment denotes Yedoma permafrost soil. The red scale color denotes organic carbon lability with talik expansion into in situ thawed Yedoma (taberites). (A) Young active lake: labile carbon results in high CH₄ production rates, with CH₄ production in taberites proportional to the talik thickness. (B) Expanding lake: over decades, new organic carbon from surrounding soils and plants is added to surface lake sediments and to freshly thawed taberites at the base of the talik, while the availability of labile organic carbon attenuates in formerly thawed taberites. (C) Mature lake: over centuries, the talik deepens (into bedrock) beyond the depth of taberites, while taberite organic carbon becomes refractory. CH₄ emissions may still be significant, depending on exogenous (non-taberite) inputs of organic carbon. At later stages (as suggested by Walter Anthony et al., 2014), the organic carbon is no longer available for methanogenesis, and the lake thus is no longer a significant source for permafrost-derived CH₄ to the atmosphere.

4.3 Methane origin in the talik

The comparison of δ¹³C${}_{{\mathrm{CH}}_{\mathrm{4}}}$ values between the profile measurements and those obtained from the incubation experiment shows less-enriched δ¹³C${}_{{\mathrm{CH}}_{\mathrm{4}}}$ values in the sediment cores across all samples (Fig. S6). These values in the sediment profiles show that the pore water CH₄ is not entirely produced in situ, suggesting that the new CH₄ produced in the lake sediments contains more isotopically enriched δ¹³CCH₄ values compared to the CH₄ diffused from deeper layers, which could be produced by a different CH₄ production pathway (Pellerin et al., 2022; Liu et al., 2025). CH₄ generated at greater depths within the talik, which subsequently ascends through diffusion or bubble transport, exhibits a more negative δ¹³C${}_{{\mathrm{CH}}_{\mathrm{4}}}$, probably due to different conditions from those prevailing in the top meter of the talik, possibly under lower metabolic rates and using different metabolic pathways (Maltby et al., 2016; Berberish et al., 2020; Pellerin et al., 2022).

By considering the depleted δ¹³C${}_{{\mathrm{CH}}_{\mathrm{4}}}$ values in the sediment cores and the corresponding incubation values at each depth, we can quantitatively assess the fraction of modern CH₄ present in the sediment profile. A comprehensive analysis of the young production fraction in all four cores shows that as depth increases, the contribution of newly produced CH₄ decreases and tends to approach zero across all cores (Fig. S6). Additionally, the new CH₄ fraction observed in Goldstream Lake (both the edge and center) is smaller compared to that of BTL. This finding suggests a greater influence of newly produced CH₄ in younger lakes compared to those that have formed and thawed longer. A shift in the methanogenesis pathway in the sediments vs. the incubations can be an alternative explanation, however less likely, given the observations here and in previous studies (Pellerin et al., 2022).

5 Conclusions

This study presents the first empirical data quantifying CH₄ production and the organic carbon degradation of thermokarst lakes from young to mature. This was achieved by quantifying the evolution of organic carbon degradation and CH₄ production rates throughout the evolution of lakes from a young dynamic lake to a mature one. The findings underscore the vertical variations in CH₄ production rates, the influence of permafrost thawing on microbial activity and the divergent patterns observed among lakes of different development stages.

Our high-resolution profiles and long-term incubations show the highest CH₄ production rates on the edges of the young BTL, then the center of BTL, then the edge of the mature GSL, with the lowest at the center of GSL. The higher rates coincided with higher TOC levels, more labile carbon content for thermal-induced chemical reactions, simpler carbon compounds and a young CH₄ source. These factors provided probably favorable conditions for microbial populations to decompose carbon, resulting in elevated CH₄ production. The high correlation between the two different methods to estimate organic carbon lability for microbial degradation CH₄ production rates and the lability to thermal-induced chemical reactions emphasizes the potential use of the Rock-Eval analyses to estimate the susceptibility of organic carbon for microbial degradation.

Our proposed conceptual mechanism, as depicted in Fig. 8, considers these two parameters (lake age and thawed talik thickness) when discussing accumulated CH₄ production rates. It can therefore be expected that the expansion of thermokarst lakes in the Arctic will continue to influence CH₄ production as the younger lakes expand into deeper permafrost layers.