Seasonal patterns in greenhouse gas emissions from thermokarst lakes in Central Yakutia (Eastern Siberia)

In the ice‐rich permafrost area of Central Yakutia (Eastern Siberia, Russia), climate warming and other natural and anthropogenic disturbances have caused permafrost degradation and soil subsidence, resulting in the formation of numerous thermokarst (thaw) lakes. These lakes are hotspots of greenhouse gas emissions, but with substantial spatial and temporal heterogeneity across the Arctic. We measured dissolved CO2 and CH4 concentrations in thermokarst lakes of Central Yakutia and their seasonal patterns over a yearly cycle. Lakes formed over the Holocene (alas lakes) are compared to lakes that developed over the last decades. The results show striking differences in dissolved greenhouse gases (up to two orders of magnitude) between lake types and seasons. Shallow lakes located in hydrologically closed alas depressions acted as CO2 sinks and strong sources of diffusive CH4 during some seasons (ebullition was not assessed). Recent thermokarst lakes were moderate to extremely high sources of diffusive CO2 and CH4, with considerable accumulation of greenhouse gas under the ice cover (winter) or in the deepest water layers (summer), highlighting the need to include spring and autumn as critical periods for integrated assessments. The water column was stratified in winter (all lake types) and especially in summer (recent thermokarst lakes), generating anoxia in bottom waters and favoring CH4 production and storage, particularly in the most organic‐rich lakes. The diffusive fluxes measured from thermokarst lakes of this typical taiga alas landscape of Central Yakutia are among the highest presented across Arctic and subarctic regions.

Permafrost occupies more than 20 million square kilometers and represents 24% of land cover within the northern hemisphere (Brown et al. 1998). It is especially abundant in Siberia, Alaska, and Canada, and its spatial extent, thickness, and ground ice content can vary widely across landscapes Strauss et al. 2017). Studies on permafrost landscape dynamics during the Holocene and in recent decades have shown that areas dominated by ice-rich permafrost are very sensitive to changes in temperature and other local disturbances Ulrich et al. 2019). Rising air and ground temperatures result in permafrost thaw, which can have widespread implications for local and regional hydrology (Biskaborn et al. 2019). Permafrost thaw can release substantial amounts of organic and mineral matter to aquatic ecosystems, causing profound changes in their biogeochemistry and their role in the global carbon cycle (Vonk et al. 2015;Tank et al. 2020 and references therein). The current rate and magnitude of temperature rise in the Arctic is disproportionately high compared to global averages, with mean annual air temperature predicted to increase by as much as 5.4 C within the next century in the absence of significant and directed global effort to reduce greenhouse gas emissions (IPCC 2019). This will likely herald a period of dynamic changes within permafrost landscapes.
A prevalent pathway of permafrost degradation in areas of ice-rich permafrost is the initiation of thermokarst processes, which eventually result in the formation of numerous lakes in regions where topography is flat .
Thermokarst processes begin when disturbances such as warming or forest removal cause deepening of the active layer (the layer which thaws each summer), which results in either thaw subsidence or thermal erosion depending on the relief of the terrain and the ground ice content. In areas that are dominated by low relief terrain, such as Central Yakutia, thermokarst subsidence often causes ponding and shallow depressions (French 2017). The coalescence and expansion (both laterally and vertically) of these ponds eventually results in the development of larger and deeper lakes, with a portion of unfrozen water remaining under the ice cover in winter . Once formed, these lakes profoundly change the local ground thermal regime, sometimes increasing surrounding sediment temperatures by as much as 10 C above the mean annual air temperature (Brouchkov et al. 2004). Thermokarst processes, especially the presence of a lake or pond, increase the rate of permafrost thaw significantly compared to what would be predicted from increases in air temperature alone (Brouchkov et al. 2004;Schuur et al. 2015). Continued lateral lake expansion can cause slope slumping (so-called "retrogressive thaw slumps") or further ground surface subsidence and erosion (Séjourné et al. 2015;Bouchard et al. 2017). Expansion and deepening continue until the lake depth surpasses the depth of ice-rich permafrost or thaw propagation beneath the lake is impeded by lake sediment insulation properties, halting further deepening. After this phase, lake size and depth are controlled mostly by the ratio between precipitation and evaporation (Soloviev 1973;Grosse et al. 2013;French 2017). Although some of the mechanisms involved are not well understood, drainage, evaporation and infilling eventually result in lake disappearance (Soloviev 1973;Desyatkin et al. 2009).
It is estimated that soils in northern permafrost landscapes contain twice as much carbon as currently exists in the atmosphere (Hugelius et al. 2014). Specifically, thermokarst lakes are considered as biogeochemical hotspots that play a pivotal role in the processing of permafrost organic matter and therefore potentially on the global climate (Walter Anthony et al. 2016). Anoxic sediments at the bottom of these lakes are sites of CH 4 production that can be released through ebullition and diffusion. The CH 4 released though ebullition in shallow lakes has little opportunity to be oxidized within the water column, so it typically represents a direct flux to the atmosphere (Bastviken et al. 2004;Bouchard et al. 2015). However, the CH 4 released through diffusive flux, or through ebullition in deeper stratified lakes, remains in the water column long enough so that a significant fraction can be oxidized to CO 2 by methane-oxidizing bacteria (Vachon et al. 2019). In winter, the CH 4 emitted through both processes can be oxidized under the ice cover (Matveev et al. 2019). Significant release of CO 2 from thermokarst lakes in northern latitudes has also been documented (e.g., Abnizova et al. 2012). There are three main pathways by which CO 2 can be released from thermokarst lakes, besides when it originates from CH 4 oxidation: (1) respiration within the oxic water column and sediments (MacIntyre et al. 2018), (2) lateral transfer of dissolved inorganic carbon from surrounding soils by groundwater flow and carbonate weathering (Zolkos et al. 2018), and (3) photo-oxidation of dissolved organic carbon (Ward et al. 2017). Lateral transport of CH 4 from the active layer was also shown to represent a significant allochthonous source in Toolik Lake, Alaska (Paytan et al. 2015). Overall, the seasonal release of greenhouse gases produced from these multiple sources is largely controlled by the ice cover and the mixing regime of lakes (Sepulveda-Jauregui et al. 2015;Matveev et al. 2019).
Little is known about the seasonality of greenhouse gas emissions from thermokarst lakes developed within permafrost. This is particularly true for winter and spring seasons that present logistical constrains (Matveev et al. 2019;Préskienis et al. 2020). Here we report on dissolved greenhouse gas concentrations, and derived diffusive fluxes, in thermokarst lakes of different development stages in Central Yakutia, Eastern Siberia, a region experiencing rapid degradation of ice-rich permafrost. The limnological properties and greenhouse gas concentrations were quantified four times throughout a full annual cycle (2018-2019). We tested the two-fold hypothesis that (1) marked local landscape heterogeneities (permafrost conditions, lake morphology and water balance) strongly affect greenhouse gas concentrations and fluxes in thermokarst lakes of Central Yakutia, and (2) these spatial patterns are further influenced by strong seasonal changes in lake water temperature, mixing regime, and dissolved oxygen.

Study area
Central Yakutia is characterized by extreme subarctic continental climate with long, cold and dry winters (January mean temperature around −40 C) and warm summers (July mean temperature around +20 C), resulting in notably strong seasonal variability. The winter season (defined by the presence of an ice cover) usually lasts from late September until early May. The low annual precipitation (190-230 mm) is generally restricted to summer. Average snow depth for winter months (January to April) ranges from 24 cm in January to a maximum of 30 cm in March, and then decreasing to 10 cm at the end of April (1980-2020 recorded values from Yakutsk weather station). Yearly evaporation rates exceed total precipitation in this region (Fedorov et al. 2014b). Central Yakutia is no exception to the observed trend of high latitude regions warming disproportionately quicker than lower latitudes. Between 1996 and 2016, the mean annual air temperature of Central Yakutia has increased by 0.5-0.6 C per decade (Gorokhov and Fedorov 2018).
The study site (62.554 N; 130.982 E) is located near the rural village of Syrdakh, which lies on the lowland plain between the Lena River to the west and the Aldan River to the east, approximately 130 km north-east of Yakutsk (Fig. 1). The region is predominantly covered by late Pleistocene sediments, including silty clays and sandy silts of fluvial, lacustrine or eolian origin (Ivanov 1984). The activity of the Lena and Aldan rivers, as well as their smaller tributaries, has resulted in the formation of numerous fluvial terraces during the Pleistocene (Soloviev 1959). Lakes later developed within late Pleistoceneaged fluvial terraces: the Tyungyulyu terrace, which covers the western section of the study area, 50-200 m above sea level (asl), dated 14-22 kyr BP, and the higher Abalakh terrace in the eastern sector of the study area, 200-280 m asl, dated 45-56 kyr BP (Soloviev 1959(Soloviev , 1973Ivanov 1984). This region is characterized as a middle taiga landscape regime (Fedorov et al. 2014a) and is dominated by larch, pine and birch forests (Ulrich et al. 2017a). Grasslands are abundant in unforested areas, such as land previously cleared for farming, ranching or in the remnant depressions of old thaw lakes known as "alases" (see below). They consist of halophytic steppe-like and bog plant communities (Ulrich et al. 2017b).
Permafrost in this region is continuous (Fig. 1a), thick (> 500 m deep), and the upper 30-50 m (Pleistocene-age fluvial and eolian sediments called "Yedoma") can be extremely rich in ground ice (50-90% by volume) (Ivanov 1984). The amount of organic carbon stored in Yedoma varies widely. For example, deep cores from Northern Siberia and Alaska yielded organic carbon pool estimates of approximately 10 + 7/−6 kg m −3 (Strauss et al. 2013). A 22 m-deep core in Central Yakutia (Yukechi) on the Abalakh terrace yielded a much lower value of organic carbon content of~5 kg m −3 (Windirsch et al. 2020), while another Central Yakutian study (Spasskaya Pad/Neleger site) of a shallow core (2 m) showed a notably higher organic carbon content of 19 kg m −3 for the top two meters of larch forest-covered Yedoma deposits (Siewert et al. 2015).
Our study site is underlain by Yedoma silty loams, which are common to the Lena-Aldan interfluve, with abundant ground ice in the form of 1.5-3 m-wide ice wedges. Active layer depth in the region generally ranges between~1 m below forested areas to > 2 m in exposed grassland areas (Desyatkin et al. 2009). Zones of unfrozen ground (or taliks) exist underneath major rivers and lakes deeper than the ice-cover thickness. Nearly half of the landscape has been affected by thermokarst since the early Holocene, resulting in the formation of thousands of partly drained alas depressions (Soloviev 1959;Brouchkov et al. 2004). However, recent thermokarst activity related to natural landscape evolution, increasing air temperatures and/or human-induced landscape modifications (agriculture, clear-cutting, and infrastructure) is also occurring in the region, as shown by the presence of numerous small and young, fast-developing lakes and retrogressive thaw slumps along lake shores (Fedorov et al. 2014b;Séjourné et al. 2015). The landscape is thus highly dynamic; yet the competing driving factors (climate vs. local geomorphology or vegetation development) and the timing of such changes (gradual vs. rapid/threshold) are complex to characterize and quantify at the regional scale.

Lake types
To understand the spatial and temporal heterogeneities in greenhouse gas concentrations from thaw lakes in the area, it was first necessary to classify the lakes at the study site (Fig. 2). Differentiation was done based on field observations, past radiocarbon dating of lake sediments, geochemical signatures of lake waters, and a multiple-stage development model (Soloviev 1973;Desyatkin et al. 2009). Relevant statistics and properties for all lake types are given in Table 1. Based on the compiled information, the lake types are: • Unconnected alas lakes: residual water bodies located within endoreic or hydrologically closed basins (Desyatkin et al. 2009). Most of these lakes likely formed during the transition between the Pleistocene and Holocene, approximately 8-10 thousand years ago (Ulrich et al. 2017b). However, the Holocene Thermal Maximum (~6700-5000 cal. yr BP) was also a time of prolific formation and continued development of new and existing alas lakes (Biskaborn et al. 2012;Ulrich et al. 2017b). These lakes can be up to a few meters deep but are typically very shallow (1 m deep or less) and are thus generally frozen to the bottom in winter. Despite a lack of liquid water in the winter, these frozen lakes have heat fluxes that can be one order of magnitude higher than those of the surrounding permafrost . The ancient lake depressions surrounding the small residual lakes of this type can be up to several kilometers wide and several meters deep. These alas lakes have already undergone much of the thermokarst processes and very little ground ice typically remains beneath the residual lake. Therefore, the thaw potential and resulting input of stored carbon to these lakes is low compared to recently formed thermokarst lakes (Ulrich et al. 2019). Yakutian people have used these depressions for agricultural purposes for centuries (Crate et al. 2017). • Connected alas lakes: lakes connected hydrologically to the watershed by streams or rivers. These lakes are consistently larger (several hundreds of meters across) and deeper (up tõ 10 m). Most of them were probably formed during the late Holocene, approximately 5-3.5 thousand years ago, although detailed chronology about their inception is still incomplete (Soloviev 1973;Ulrich et al. 2017b). Local people currently use some of these lakes for fishing (e.g., Syrdakh Lake; Fig. 1c).
• Recent thermokarst lakes: thaw lakes formed over the last several decades mostly from human activities (e.g., forest fire and forest removal for agriculture, pipelines, or road construction) and rising temperature (Fedorov et al. 2014b). These lakes are generally small (meters to tens of meters across) and relatively shallow (generally one to two meters deep), and are still expanding downwards and laterally due to active layer deepening and thermokarst processes. Their waters can be notably rich in dissolved organic carbon (DOC), with concentrations typically ranging from~100 to 300 mg L −1 (Hughes-Allen et al. 2020).  Table 1; Supp. Fig. S1). Physicochemical measurements were conducted during each of the four field campaigns using a YSI Pro DSS multiprobe sensor, including profiles in temperature (accuracy of AE 0.2 C), dissolved oxygen (accuracy of AE 1% of reading or 1% saturation, whichever is greater), specific conductivity (accuracy of AE 1% of reading) and pH (accuracy of AE 0.2 pH units). Vertical profiles of the water column were obtained from 17 lakes during the winter through a hole drilled in the ice cover using a hand-held auger, and from 3 lakes during the summer from an inflatable boat (Table 1). All other data were collected from the lake surface as discrete point measurements, mostly from the lakeshores.

Dissolved greenhouse gas measurements
In spring, summer, and autumn, the water was mostly sampled from the lakeshore at the surface. Summer samples were also collected from the pelagic section of the lake at the surface and bottom using a horizontal van Dorn sampler (integrating 15 cm depth) for recent thermokarst Lake 1 only. In winter, samples were collected from the center of each lake through the ice cover using a vertical van Dorn sampler (integrating~40 cm depth), including a few supplementary samples taken from bottom waters. Water was quickly transferred into a 2-liter gas exchange LDPE bottle (if not directly sampled from this bottle), and a headspace of 30 mL of atmospheric air was created (bottle with an inlet at the bottom, connected to a syringe). The bottle was shaken vigorously for 3 min, and 20 mL of the gaseous headspace was then transferred to a 12 mL pre-vacuumed gastight Exetainer glass vial previously flushed with nitrogen (Hesslein et al. 1991). Samples were later analyzed by gas chromatography (Thermo Trace 1310, TRI-Plus Head-Space auto-sampler, HSQ 80/100 1.6 mm × 30.5 m in series with MS 5A 1.6 mm × 18.3 m, thermal conductivity and flame ionization detectors). The dissolved greenhouse gas concentration (C sur ) was calculated using Henry's Law, and departure from saturation was obtained subtracting the gas concentration in the water at equilibrium with the atmosphere (C eq ) using global annual average values of 410 ppm for CO 2 and 2 ppm for CH 4 (IPCC 2019).
Diffusive fluxes of dissolved greenhouse gas were estimated as in Bouchard et al. (2015) based on a gas transfer coefficient (k 600 ) standardized to a Schmidt number (Sc) of 600 (Wanninkhof 1992), and calculated with the wind-based model of Vachon and Prairie (2013): k 600 = 2:51 + 1:48 U 10 + 0:39 U 10 log 10 LA where U 10 is the wind speed at 10 m above the ground from the closest meteorological station (at Yakutsk) and LA is the lake surface area (in km 2 ). The wind speed at Yakutsk airport was very well correlated to the weather station at Borogontsy, which is less than 15 km from most lakes in the study site (monthly r ! 0.9 for years 1960-1987; unpublished data from A. Fedorov). The flux (in mmol m −2 d −1 ) was then calculated by applying the equation: where k is the transfer coefficient for a given gas calculated as: Due to generally high pH (> 9) and low wind turbulence during the ice-free seasons, chemical enhancement could further increase the high levels of CO 2 uptake estimated for some lakes (Wanninkhof and Knox 1996). The exchange rate of gases from the atmosphere through the boundary layer of a liquid surface is regulated by molecular diffusion and accelerated by increasing turbulence (wind). However, the exchange rate of CO 2 molecules from the atmosphere through the boundary layer of a lake surface is also regulated by water temperature, pH and ionic strength. When favorable conditions are met (high pH and temperature, low winds), diffusion across the boundary layer can be enhanced by the hydration of CO 2 to HCO 3 − . The CO 2 uptake is therefore a function of both diffusion and chemical reaction. The chemical enhancement factor for CO 2 (α) is defined as: where k enh is the enhanced gas transfer velocity, and k is the gas transfer velocity for a nonreactive gas with the same diffusion coefficient as CO 2 as calculated above (Wanninkhof and Knox 1996).

Statistical analysis
The dissolved greenhouse gas data were not found to be normally distributed and did not meet the assumption of homoscedasticity based on the Shapiro-Wilk normality test. Typical data transformations to improve normality, such as a Box-Cox and logarithmic transformations, were attempted with unsuccessful results. The small sample sizes likely impede the ability to achieve normality, and we did not expect the data to be inherently normal. Therefore, we used Kruskal-Wallis one-way ANOVA to test the relationships between dissolved CO 2 and CH 4 concentrations vs. lake types and seasons. Significance between groups was tested for dissolved CO 2 and CH 4 concentrations as follows: the values for dissolved greenhouse gas concentrations for each lake type were compared across seasons (e.g., dissolved CO 2 concentrations for unconnected alas lakes between fall, winter, spring, and summer) and between lake types for a specific season (e.g., dissolved CO 2 concentrations for unconnected alas lakes, connected alas lakes, and recent thermokarst lakes during the fall season). The null hypothesis that the medians of the compared groups are equal was rejected if p < 0.05. All of the statistical analysis was completed using the Python programming language (Python Software Foundation, https://www.python.org/).

Seasonal conditions
Sampling began in September 2018, at the end of the icefree season when the depth of the active layer is at its maximum. During this period, nighttime air temperatures drop near 0 C, initiating lake water mixing (fall overturn).
Although limnological profiles were not done during the fall sampling season, water surface temperatures (ranging from 4.7 to 12.5 C depending on the lake; Hughes-Allen et al. 2020) indicate that some lakes may have been sampled before the water column had completely mixed to bottom. During the winter, lakes were covered by 0.7-1.2 m of ice and 30-50 cm of snow (measured in March-April 2019). Shallow lakes, mostly unconnected alas lakes, were frozen completely to the bottom during the winter, whereas conditions were met for inverse stratification in deeper lakes (defined as the water bodies keeping liquid water during winter, i.e., of a maximal depth > 1.2 m). Sampling occurred late in winter (< 3 April), but well before the thaw season, which usually begins in May in this region. Many lakes (larger and/or deeper) were still covered in ice at the beginning of May when the site was visited for the spring sampling, but the ice cover melted progressively and initiated the spring overturn until surface water reached 4 C. The surface temperature of most lakes was already above 4 C at the time of sampling in spring (ranging from 2.8 C to 17.5 C; Hughes-Allen et al. 2020), with the smaller unconnected alas lakes showing consistently higher surface temperatures (mean 10.9 C) than the connected alas lakes (5.4 C) or recent thermokarst lakes (6.8 C). Therefore, spring sampling likely occurred after the initiation of the thermal stratification. The warming of surface waters during the summer initiated the thermal and chemical stratification of the water column that may have lasted throughout the summer for deeper lakes.

Physicochemical characterization of lake water Broad trends and seasonal averages
The hydrochemical characteristics of the sampled lakes differed notably between lake types and seasons (surface values presented in Fig. 3). Regardless of the season, unconnected alas lakes generally had the highest electrical conductivity. This inter-lake difference was particularly marked during the winter, when the mean conductivity value for unconnected alas lakes was 5711 μS cm −1 . The conductivity for unconnected alas lake #3 reached 8000 μS cm −1 , four times higher than the average conductivity for connected alas lakes. Connected alas lakes showed the lowest values (winter mean = 563 μS cm −1 ), with recent thermokarst lakes being intermediate (winter mean = 4228 μS cm −1 ) ( Fig. 3a-d). As expected, conductivity reached much lower values (averages < 1000 μS cm −1 ) for all lakes in the springtime, during and shortly after ice-cover breakup. Conductivity increased slightly during the summer, reaching mean values similar to the previous fall (~1000-2000 μS cm −1 ) and showing comparable values between recent thermokarst and unconnected alas lakes.
The percent saturation of dissolved oxygen for all three lake types was lowest in the wintertime, with mean values at lake surface between 0% and 20% (Fig. 3e-h). Unconnected alas lakes maintained the highest levels of dissolved oxygen during fall, spring, and summer compared to the other two lake types (with more pronounced differences in spring), reaching values above 200% in summer (Fig. 3h). Recent thermokarst lakes maintained the lowest concentrations of dissolved oxygen, rarely reaching saturation levels even in the springtime. Connected alas lakes followed a similar trend as the other two lake types, but with intermediate values.
The pH values exhibited less seasonal variability compared to the other hydrochemical properties analyzed (Fig. 3i-l). Unconnected alas lakes consistently had the highest pH levels relative to the other two lakes types, maintaining an average pH between 9.4 and 10 for all four seasons. Connected alas lakes and recent thermokarst lakes exhibited similar seasonal trends, with average pH values in fall, spring and summer near 8.8, and with lower values during the winter season (average 8.7).

Seasonal profiles
Although temperature and dissolved oxygen profiles are not available for all lakes for each season, analysis of the profiles for selected representative lakes (unconnected alas lake #3, connected alas lake #2, recent thermokarst lake #1) (see lake locations in Supp. Fig. S1) illustrates some broad seasonal trends. Profiles of temperature and dissolved oxygen show that the deeper lakes (i.e., connected alas and some recent thermokarst lakes) were stratified during winter and summer (Figs. 4, 5). Inverse stratification was observed in winter, with bottom water ranging from near zero (at 2.1 m in unconnected alas #3) to 4 C (at 8.2 m in connected alas #2). During this period, dissolved oxygen saturation levels ranged between~10% and 40% in surface waters (0-50 cm), while it decreased to values < 5% in bottom waters (Fig. 4a,c,e). Conductivity and pH profiles showed smaller differences between surface and bottom waters during winter (Fig. 4b,d,f). These profiles nevertheless illustrate the broad range of conductivity (496-8066 μS cm −1 ) and pH (7.4-9.4) that can be found between individual lakes during the winter on such a lake-rich landscape.
In August, available temperature and oxygen profiles for recent thermokarst lakes (lakes #1, #8, and #18) show strong stratification (Fig. 5). At the time of sampling (0945-1200), the difference between surface and bottom temperature was 7.1 C (lake #1), 5.8 C (lake #8), and 7.3 C (lake #18). Dissolved oxygen saturation levels ranged between~28% and 54% in surface waters and between~4% and 11% in bottom waters (Fig. 5). All three recent thermokarst lakes showed similar profiles in pH and conductivity. Profiles showed decreasing pH values with depth (ranging from approximately 9 at the surface to 7.5 in bottom waters). Conductivity profiles showed an increase with depth from approximately 2500 μS cm −1 at the surface to 6000 μS cm −1 in bottom waters for all three sampled recent thermokarst lakes.

Dissolved greenhouse gas concentrations
Surface concentrations of dissolved greenhouse gas, expressed as departure from saturation, varied strongly between lake types and seasons (Fig. 6). Broadly, unconnected alas lakes and recent thermokarst lakes had stronger seasonal differences in CO 2 concentrations, while connected alas lakes showed less seasonal difference in CO 2 concentration. Winter CO 2 concentrations for all lake types were generally significantly different from all other seasons. Differences between unconnected alas lakes and recent thermokarst lakes were generally significant for most seasons (Table 2; Table 3). Overall, recent thermokarst lakes consistently showed the highest CO 2 concentrations between lakes types, yet values varied considerably among the 15 lakes sampled and between seasons, with departure from saturation averages ranging from slight CO 2 sink during spring sampling (− 10 μM) to very strong CO 2 source during winter sampling (+ 1440 μM).
During the fall sampling period, unconnected alas lakes were under-saturated in CO 2 , and thus acting as CO 2 sinks relative to the atmosphere, whereas recent thermokarst lakes were generally super-saturated and therefore sources of CO 2 (Fig. 6a). The CO 2 concentrations in connected alas lakes were close to equilibrium with the atmosphere and thus neither a source nor a sink at time of sampling. A remarkable increase in CO 2 concentrations occurred under the ice cover for all lake types. For example, recent thermokarst lake #34 had 20 times more CO 2 (1295 μM) in the winter than in the previous fall (62 μM) and connected alas lake #11 had 10 times more CO 2 in the winter compared to the previous fall (190 μM in winter vs. 27 μM in fall) (Fig. 6c). Bottom water concentration of dissolved CO 2 was only sampled in the connected alas lake #2 during winter, showing a concentration similar than at the surface (155 μM in bottom water vs. 168 μM at the surface).
During spring sampling, nearly all lakes (28 out of 33) had CO 2 concentrations below saturation, with relatively low variability (Fig. 6e). During late summer however (August), saturation levels at the lake surface returned to similar values as observed during the previous fall (early fall, September); unconnected alas lakes reverted to slight CO 2 sinks, connected alas lakes returned to equilibrium, and recent thermokarst lakes became the largest CO 2 sources again, with one particular lake (#35) showing a notably high value of nearly 500 μM (Fig. 6g) (Supp. Fig. S1). During summer sampling, the bottom waters of recent thermokarst lake #1 (Supp. Fig. S1) showed extremely high saturation levels (> 1500 μM), nearly two orders of magnitude higher than at the surface, representing a substantial storage of CO 2 . Differences between surface and bottom water were not assessed in other lakes.
General trends in the concentrations of CH 4 showed that recent thermokarst lakes typically had stronger seasonal differences in CH 4 concentrations compared to the other two lake types. The summer sampling season showed the highest variability in CH 4 concentration between lake types than the other three seasons ( Table 2; Table 3).
During fall sampling, unconnected alas lakes were the strongest sources of CH 4 , while the other lake types were small or negligible sources (Fig. 6b). During winter, as for CO 2 , dissolved CH 4 showed a large increase in concentration for recent thermokarst and connected alas lakes (by one order of magnitude) compared to the previous fall (Fig. 6d). Bottom water concentration of dissolved CH 4 was only sampled for connected alas lake #2. The bottom water concentration was one order of magnitude higher than surface concentration during the winter sampling period (133 μM in bottom water vs. 15 μM at the surface).
During spring sampling, all lakes (excluding a small number of outliers) were moderately saturated in CH 4 compared to the values observed during winter sampling (nonoutliers < 60 μM), with similar saturation levels between lake types (Fig. 6f). During summer sampling, observed CH 4 saturation levels at the surface were lower (~10 μM), with unconnected alas lakes as the largest CH 4 sources compared to the other lake types, which is a similar trend to the previous fall (Fig. 6h). As observed for CO 2 during summer, recent thermokarst lake #35 showed notably higher surface CH 4 concentrations (saturation levels > 700 μM, identified by an arrow in Fig. 6h) compared to all other sampled lakes. Bottom waters of recent thermokarst lake #1 (Supp. Fig. S1) showed extremely high saturation levels during summer sampling (> 40 μM), nearly 130 times higher than at the surface (~0.4 μM), representing a substantial storage of CH 4 . Differences between surface and bottom water were not assessed in other lakes.

Diffusive greenhouse gas fluxes
Considering only the ice-free season (from spring to fall; Table 4), the diffusive CO 2 fluxes ranged from a minimum of − 13.6 mmol m −2 d −1 (or~− 600 mg C m −2 d −1 ) recorded in spring for unconnected alas lake #3 to a maximum of 355 mmol m −2 d −1 (or~1200 mg C m −2 d −1 ) recorded in summer for recent thermokarst lake #18. The ice-free season median CO 2 fluxes were relatively low for connected alas lakes and negative for unconnected alas lakes (1.5 and − 4.8 mmol m −2 d −1 , respectively), and slightly positive for recent thermokarst lakes (1.8 mmol m −2 d −1 ). The diffusive CH 4 fluxes were highest from recent thermokarst lakes in summer, reaching~560 mmol m −2 d −1 (or 6720 mg C m −2 d −1 ). It is noteworthy that the average ice-free season CH 4 flux for all studied lakes was substantially higher than the median value (14 mmol m −2 d −1 compared to 3.1 mmol m −2 d 1 ), reflecting the skewing effect of exceptionally high fluxes in summer from a few recent thermokarst lakes. Unconnected alas lakes were regularly large CH 4 sources during the ice-free season (median flux of~7.0 mmol m −2 d −1 , compared to 1.0 mmol m −2 d −1 for recent thermokarst lakes). It should be noted that these results are estimations based on wind data from the only available meteorological station in the region (Yakutsk Airport), located~120 km from our study site.
It is possible that the chemical enhancement effect contributed to increased CO 2 uptake rates by some lakes during the ice-free seasons. Unconnected alas lakes consistently had the highest α values (mean value of 9.9 for unconnected alas lakes, reaching up to 27 for lake #20 in August). These large values occurred when pH was particularly high, lake water was warm and wind speeds were low. Alpha values calculated for the other two lakes presenting a CO 2 uptake during August were lower (connected alas lake #11 α = 1.6; recent thermokarst lake #37 α = 5.3). During September, there were only two connected alas lakes and two recent thermokarst lakes that were acting as CO 2 sinks (mean α = 1.4), while the unconnected alas lake #10 presented the highest α value (3.8). The spring sampling season showed low variability in the calculated α values between lakes types, although unconnected alas lakes still had the highest mean value (2.8, compared to mean values~1.0 for connected alas lakes and recent thermokarst lakes). It is important to note that during the spring season, most lakes (26 out of 31) were acting as CO 2 sinks and the magnitude of CO 2 uptake was highest compared to the other seasons (Fig. 6). However, the highest chemical enhancement of CO 2 flux was calculated for the summer season, when CO 2 uptake is less prevalent and smaller in magnitude (Fig. 6).

Discussion
The study area represents the typical taiga alas landscape found in Central Yakutia, a region encompassing approximately 17,000 km 2 in the Lena-Aldan interfluve (Ulrich et al. 2017a). Nearly 10 % of this landscape is comprised of lakes at varying stages of thermokarst development (Ulrich et al. 2017a). Lake morphology is directly linked to the evolution stage within the thermokarst sequence (Soloviev 1973), with rates of change substantially lower for older lakes (Desyatkin et al. 2009). The morphology of unconnected alas lakes has remained relatively stable during the last decades, with little lateral expansion or deepening. In contrast, recent thermokarst lakes have been notably dynamic, experiencing subsidence and lateral expansion consistent with the rapid landscape evolution observed since the mid-20th century (Fedorov et al. 2014a;Ulrich et al. 2017a). Connected alas lakes, while they did not notably expand or deepen in the recent past, have shown some signs of thermokarst activity such as thaw slumps (Séjourné et al. 2015). This spatial heterogeneity resulting from the above chronological sequence can contribute to notable differences in lake water physicochemical characteristics, such as timing of spring ice-cover breakup, water temperature and lake water mixing regime, electrical conductivity, and dissolved oxygen concentration.
Depending on the lake type, we also found substantial seasonal variations in the amount of dissolved greenhouse gas in lake water, spanning two orders of magnitude throughout the year. Since most greenhouse gas studies are limited to summer quantification (e.g., Abnizova et al. 2012;Bouchard et al. 2015), such seasonal assessments are precious (Langer et al. 2015;Matveev et al. 2019). The wide seasonal ranges observed here underscore the need to consider such heterogeneities when attempting to estimate global greenhouse gas emissions from thermokarst lakes and permafrost landscapes. Developmental stage as a driving factor on lake greenhouse gas concentrations and fluxes We measured notable differences in physicochemical properties (Fig. 3) and dissolved greenhouse gas concentrations (Fig. 6) between thermokarst lakes of different development stages. The clear (statistically significant) difference between unconnected alas lakes and the other two types of lakes underlines the unique role played by these shallow lakes, acting as CO 2 sinks most of the year and substantial CH 4 sources, especially in fall. Several studies have concluded that unconnected alas lakes are generally not sites of current thermokarst activity (e.g., Brouchkov et al. 2004;Ulrich et al. 2017a), which is consistent with our observations. Little ground ice remains below alas lake depressions, which significantly reduces the potential for input of stored carbon from surrounding permafrost (Ulrich et al. 2019). They are stable both laterally and vertically, and any variations in lake volume can be attributed to changes in the hydrological regime. As they are hydrologically isolated, the interannual water balance in these lakes is entirely dependent upon the ratio between annual precipitation and evaporation (Brouchkov et al. 2004;Ulrich et al. 2017a). Based on these characteristics, evaporative enrichment during summer and substantial solute exclusion under winter ice cover are likely to play a role in the high concentration of nutrients and minerals observed within these lakes, such as found in Western Siberia (Manasypov et al. 2015). The notably high organic carbon concentrations (generally DOC > 100 mg L −1 except in spring) are characteristic of these lakes (Hughes-Allen et al. 2020).
In the open-water seasons, unconnected alas lakes are hotspots of biological activity, sometimes resulting in the entire surface of the lake being covered by floating aquatic plants (e.g., Lemna sp.) (Malyschez and Peschkova 2001). These lakes are particularly shallow (< 2 m), generally warmer than the other two lake types, and rich in nutrients. These characteristics favor high rates of photosynthetic activity, causing unconnected alas lakes to be efficient CO 2 sinks (or negligible CO 2 sources) during the open-water season. Negative fluxes are  g and h). Scales for CO 2 concentrations appear on the left, and for CH 4 concentrations on the right. Number of lakes sampled for each lake type for each season is presented in Table 1. The box indicates the interquartile range. The horizontal line within the box is the median and the small x is the mean. The whiskers indicate minimum and maximum. most prominent during the spring and fall seasons, probably linked to the more turbulent conditions over these periods leading to an efficient CO 2 transfer from the atmosphere to the lake water. It should be noted that the CO 2 uptake rates for these lakes might be underestimated. Chemical enhancement calculations indicate that CO 2 uptake could be much higher during the summer (on average by a factor of 10) when conditions are met for this effect to be significant (high pH, high temperature, and low winds). These lakes are presenting the highest pH values (generally above 9 and up to 10.5).
Unconnected alas lakes also act as CH 4 sources throughout the year, particularly in fall with flux reaching~16 mmol m −2 d −1 (Table 4). These fluxes are generally much higher than what has been reported elsewhere across the Arctic (e.g., Abnizova et al. 2012;Bouchard et al. 2015;Matveev et al. 2016). High production of CH 4 throughout the openwater season is likely linked to the high levels of primary productivity fueling microbes with recent, labile organic carbon to the system. Recently produced, autochthonous organic carbon may be more readily decomposed under anoxic conditions than organic carbon from a terrestrial source (Grasset et al. 2018). Consequently, we hypothesize that CH 4 production is positively linked to the presence of autochthonous organic carbon in these lakes.
Unconnected alas lakes are shallow water bodies less strongly stratified during winter compared to the other lake types (deeper and more humic). Although limnological profiles were not collected during summer for unconnected alas lakes, it is likely that they are not as stratified as other lake types in summer as well due to their very shallow depth and susceptibility to wind-induced mixing. These conditions likely contribute to an earlier release of the greenhouse gas stored in the water column or sediments in spring and during summer. The mean surface temperature of unconnected alas lakes in May was 11 C compared to 5 C (connected alas lakes) and 7 C (recent thermokarst lakes), indicating that summer stratification started earlier for these lakes (Hughes-Allen et al. 2020) despite they are likely polymictic lakes. It is possible that they were presenting positive CO 2 flux and higher CH 4 flux earlier in spring, but since most unconnected alas lakes freeze to bottom in winter, storage flux would be mainly through greenhouse gas stored in the ice itself, or in the interstitial water of the sediment. In summer, periods of wind-induced mixing would likely generate regular venting of the greenhouse gas produced in summer by these shallow lakes, although the very high CH 4 concentrations observed in fall indicate either some release of CH 4 stored at the bottom or very high production rates during this period. Other studies on thermokarst lakes have documented that the timing and magnitude of greenhouse gas storage flux can vary between lakes even when they are subject to the same local climate conditions (Matveev et al. 2019). Table 2. Results of the Kruskal-Wallis significance tests which compared dissolved greenhouse gas concentrations for each lake type between seasons. The results for CO 2 and CH 4 are grouped by lake type. The left-most and right-most columns indicate which seasons are being compared (e.g., winter CO 2 was larger than fall CO 2 for unconnected alas lakes). Significance (p < 0.05) is indicated by > or < signs, while = signs indicate that differences are not significant.

Unconnected alas lake
Connected alas lake Recent thermokarst lake Table 3. Results of the Kruskal-Wallis significance tests which compared dissolved greenhouse gas concentrations of the three lake types during one season. The leftmost and rightmost columns indicate the lake types being compared. Significance (p < 0.05) is indicated by either > or < signs, depending on which lake type has the highest dissolved greenhouse gas concentration for that particular season, while = signs indicate that differences are not significant. Between the three lake types, recent thermokarst lakes generally showed the highest variability in dissolved greenhouse gas concentrations among lakes of this type (Fig. 6). This potentially results from the substantial diversity in morphology and thermokarst developmental stage for this lake type, with some lakes being much deeper and larger than others (depths from 1.2 to 4.6 m, surface areas from 0.1 to 2.0 ha) ( Table 1). Morphology has been suggested to affect the mixing regime, and thus the seasonal patterns in greenhouse gas production vs. consumption, storage and emissions (Matveev et al. 2019). For example, the spring overturn of shallower lakes may represent a more complete outgassing of all greenhouse gas stored under the ice in winter, while deeper lakes may only experience outgassing of the greenhouse gas stored at the lake surface, with deeper storage only venting out in the fall 2019. Morphology as thus been identified as a major characteristic controlling seasonal emission patterns, and deserving attention in future studies.
High concentrations of both CO 2 and CH 4 were measured in recent thermokarst lakes, especially under the ice cover (in winter) but also during spring and summer in some cases (Fig. 6). Notably, thermokarst lake #35 (Supp. Fig. S1), currently developing within a forested ground and away from any apparent human influence, had the highest surface concentrations observed in summer (> 450 μM of CO 2 and 1000 μM of CH 4 ) (Fig. 6g,h). Although a limnological profile is not available for this lake, strong summer stratification is expected such as is observed in other recent thermokarst lakes (Fig. 5). Finding such high concentrations of CO 2 and CH 4 at the surface of a strongly stratified lake would indicate that greenhouse gas production was elevated in the pelagic portion of the lake, or that lateral transfer of littoral benthic production was efficient. It is also possible that this specific sampling day was following a partial mixing event bringing greenhouse gas enriched waters at the surface. The wind speed during the 2 d prior to greenhouse gas sampling reached sustained speeds up to 10 m s −1 , which is much higher than typical wind speeds in this area, usually near 1 m s −1 (Yakutsk Airport Metrological Station). These conditions could have contributed to partial mixing of the lake waters and the observed high concentrations of greenhouse gas in the surface waters. Overall, recent thermokarst lakes were particularly high CO 2 emitters in summer and fall (Fig. 6), and lower CH 4 emitters than unconnected alas lakes (particularly in fall). In fact, CO 2 /CH 4 molar ratios stayed relatively high over the summer and fall in recent thermokarst lakes (> 40) (Table 5), potentially indicating efficient methanotrophy in these systems (Crevecoeur et al. 2017). It is also possible that such elevated ratios were caused by a lateral transfer of CO 2 produced in soils under aerobic conditions (Campeau et al. 2017) or from the chemical weathering of carbonate (Zolkos et al. 2018), sources that would be particularly significant under active thermokarst erosion. The greenhouse gas production in recent thermokarst lakes is likely fueled by the thawing of ice-rich Pleistocene permafrost. Recent thermokarst lakes experience high rates of lateral and vertical expansion (Fedorov et al. 2014a;Ulrich et al. 2017a), and once initiated, this expansion occurs relatively continuously. However, the rate of expansion is strongly dependent on vegetation cover type, the presence of disturbed landscape surface (i.e., clearing of land for agriculture or infrastructure), and local permafrost conditions (Lara et al. 2015;Ulrich et al. 2017a). The expansion of these lakes into surrounding permafrost results in the mobilization of organic and mineral matter which was sequestered for thousands of years by freezing temperatures (Hugelius et al. 2014;Strauss et al. 2017). For example, DOC values ranging from 8 to 838 mg L −1 were measured in 2018-2019 (Hughes-Allen et al. 2020). Previous work in the region suggests that old carbon is released within these lakes in the form of century-to millennia-old dissolved inorganic carbon, which is resulting from variable mixtures of late Pleistocene (20 kyr BP) and more recent or modern carbon (Soloviev 1959(Soloviev , 1973Ivanov 1984). Radiocarbon dating and stable isotopic signatures of gases, lake sediments and surrounding soils should provide more insight about the sources and pathways of the greenhouse gas emitted from recent thermokarst lakes of this region (Bouchard et al. 2015).

Seasonal variations in greenhouse gas concentrations
Seasonal differences in dissolved greenhouse gas concentrations were particularly striking at the studied site. The most notable seasonal change was the substantial increase (up to two orders of magnitude) in both CO 2 and CH 4 during the winter season (Fig. 6). Increasing concentrations under the ice cover have been noted in other Arctic, subarctic and boreal regions (Langer et al. 2015;Matveev et al. 2019). In productive systems, the isolation of the water mass under an ice cover quickly generates anoxic conditions favoring methanogenesis and suppressing methanotrophy. Our studied lakes showed median CO 2 /CH 4 molar ratios that were much lower in winter (2.2 for all lakes; ranging from 0.2 for unconnected alas lakes to 10 for connected alas lakes) compared to summer (33 for all lakes; ranging from 4.4 for unconnected alas lakes to 56 for recent thermokarst lakes) (Table 5). This suggests a strong production of CH 4 and the predominance of hydrogenotrophy (consuming CO 2 for methanogenesis) in the anoxic conditions of winter, such as proposed in other thermokarst lakes of subarctic peatlands (Matveev et al. 2019).
The greenhouse gases accumulated in the water column in late winter will be partly or completely released in spring, as soon as the ice cover degrades, a period when sampling is particularly challenging logistically. The spring overturn period has been shown to be particularly important for greenhouse gas release from lakes to the atmosphere (Langer et al. 2015). Phelps et al. (1998) found that greenhouse gases released during spring overturn accounted for as much as half of the total yearly emissions from Alaskan lakes. Our sampling occurred in May, likely after the major venting period, since the ice cover had already disappeared in many lakes (especially in unconnected alas lakes, which were all ice-free in early May) and surface waters had already started to warm (26 lakes were presenting temperatures above 5 C, reaching up to 17.5 C). This likely explains why lake water did not present high saturations levels during this sampling period despite the high concentrations observed 1 month earlier under the ice cover. In fact, CO 2 concentrations were mainly below saturation in May, indicating that primary producers had already started to be active.
In summer, several recent thermokarst lakes were strongly stratified (examples given in Fig. 5), isolating bottom waters where oxygen depletion occurs and greenhouse gases accumulate. During August sampling, sufficient time had passed after initiation of summer stratification to observe a significant accumulation of greenhouse gases in bottom waters; for example, more than 1500 μM of CO 2 and 50 μM of CH 4 were measured at the bottom of recent thermokarst lake #1. This summer storage of greenhouse gas will likely be released to the atmosphere during the fall overturn period, in addition to greenhouse gas produced during the fall in sediments and water column, which can be significant considering the time lag before water and sediment temperatures cool down after air temperature has dropped (Serikova et al. 2019). This could explain the super-saturation in CO 2 (recent thermokarst lakes) and CH 4 (unconnected alas lakes) measured in fall (Fig. 6). Fall sampling may have been too early for the highly stratified thermokarst lakes, and thus the cumulated CH 4 at their bottom would be venting later. If this is the case, surface concentration (and emissions) of CO 2 would also be higher at the surface later in fall for these lakes. By contrast, connected alas lakes are relatively well mixed during the summer, as they are large water bodies influenced by wind-induced turbulence and water flows from regional streams and rivers (Fig. 2). These characteristics likely lead to the continuous release of CO 2 and CH 4 produced within the lake (and intermediate gas concentrations, Fig. 6), rather than an isolated release of the storage flux during the fall overturn period. It is important to note that the intervals between sampling periods were not symmetrical throughout the year, and particularly short between summer (mid-late August) and fall (early September), explaining their similarities in greenhouse gas concentrations. These differing patterns underline the importance in placing automated sensors in lakes to better identify critical periods of greenhouse gas emissions (Matveev et al. 2019).

Diffusive greenhouse gas fluxes: Comparison across highlatitude regions
Diffusive greenhouse gas fluxes from lakes studied across Arctic and subarctic regions are generally limited to the icefree season, with only a few studies providing estimations from wintertime or over a year cycle (Langer et al. 2015;Matveev et al. 2019). Our estimations for diffusive greenhouse gas fluxes are based on four discrete periods over a full year cycle, including spring and fall, which have been previously identified as critical periods for greenhouse gas evasion (Jammet et al. 2015;Matveev et al. 2019). The diffusive CO 2 fluxes estimated for this study site of Central Yakutia (range from − 14 to 350 mmol m −2 d −1 ; Table 4) are clearly outside the range reported for lakes and ponds across the Arctic. The most comparable fluxes are measured from ponds on Bylot Island (Bouchard et al. 2015;Préskienis et al. 2020) and on the Mackenzie Delta (Tank et al. 2008). The estimated diffusive CH 4 fluxes are also one order of magnitude higher (up to 560 mmol m −2 d −1 ) than published values across the Arctic (maximum values below 20 mmol m −2 d −1 for different types of lakes, including thermokarst lakes; Wik et al. 2016) and almost as high as low seep ebullition fluxes from erosive margins of thermokarst lakes in Alaska and Siberia (Walter Anthony et al. 2010). An important point to consider here is that our estimations only include diffusive fluxes, while many studies have shown that ebullition can represent a significant mode of CH 4 evasion (Bastviken et al. 2011;DelSontro et al. 2016), particularly from Arctic lakes (Walter et al. 2007). Some studies report that ebullition fluxes can be equivalent to diffusive fluxes, if not largely dominant in certain conditions (Walter Anthony et al. 2010;Wik et al. 2016), while others found that diffusion dominates in shallow thermokarst lakes (Matveev et al. 2016;Préskienis et al. 2020). The diffusive fluxes presented here already being in the upper range of published values for thermokarst lakes across the circumpolar North, they can only increase when ebullition is included, underscoring the need to consider such landscapes in global assessments.

Conclusions
Our study site in Central Yakutia is within an area of thick, ice rich permafrost, with the potential for substantial carbon mobilization in the face of continued climate change. There is considerable landscape variability and the thermokarst lakes in this region vary in age from early Holocene to the last decades. We observed temporal and spatial heterogeneity (up to two orders of magnitude) in the concentrations and fluxes of dissolved greenhouse gases from the three development stages of thermokarst lakes found in the region. There were also large differences in the relative importance of CO 2 and CH 4 emissions, which is fundamental for assessing their global warming potential. These differences are related in part to the mixing regime that is controlled by lake morphology and leading to variable seasonal patterns. They are also linked to variable contributions from primary producers and thermokarstic erosion shaping light and oxygen availability. We show that diffusive fluxes of both CO 2 and CH 4 from thermokarst and alas lakes of Central Yakutia are among the highest presented across Arctic and subarctic regions. With their extreme bottom greenhouse gas concentrations and large but variable emissions, recent thermokarst lakes need to be closely considered as they likely involve the mineralization of ancient organic carbon that can contribute to the amplification of the greenhouse effect. On the other hand, unconnected alas lakes act as active CO 2 sinks where primary production fuels methanogeny, likely dampening total emissions from the region depending on the overall greenhouse gas budget. The high levels of temporal and spatial heterogeneity between seasons and lake types illustrate the complexities of such permafrost landscapes and their responses to changes in temperature, precipitation, and anthropogenic influences. The fate of the large amounts of organic carbon stored in permafrost landscapes of Central Yakutia needs to be assessed to identify in which cases lakes contribute to its transfer to the atmosphere, especially considering the continued influence of climate change on these sensitive environments.