The lagoons investigated here are the semi-closed Uomullyakh Lagoon (UoL) and the nearly closed Polar Fox Lagoon (PFL) (Table 1), located on the southern coast of the Bykovsky Peninsula (71°51′ N, 129°19′ E) (Strauss et al., 2018, Figure 1D). The peninsula is located in Northeast Siberia (Figure 1B), at the mouth of the Bykovskaya Channel, one of the main outflow channels of the Lena River Delta (Figure 1A). The hinterland of the Bykovsky Peninsula is the Kharaulakh Range with elevations up to 500 m above sea level (asl), which is part of the northern foothills of the Verkhoyansk Mountains. The peninsula consists of the remains of a former late Pleistocene sediment accumulation plain with elevations of up to 45 m asl (Schirrmeister et al., 2002; Schirrmeister et al., 2018). It lies within the continuous permafrost zone and is characterized by thick, ice-rich Yedoma Ice Complex deposits (Grosse et al., 2007). Postglacial sea level rise during the Lateglacial interstadial and the early Holocene inundated the shallow Laptev Sea shelf and transformed terrestrial permafrost into subsea permafrost (Romanovskii et al., 1998; Lantuit et al., 2011). Cliffs and thermokarst basins are the typical backshore coastal landforms. Sixteen large thermokarst depressions are found on the peninsula (Grosse et al., 2005; Fuchs et al., 2018). In total, more than 50% of the Bykovsky Peninsula is strongly influenced by thermokarst and thermo-erosion landforms (Grosse et al., 2005; Lantuit et al., 2011).

According to the Köppen–Geiger world climate classification (Peel et al., 2007), the peninsula is located in the polar tundra climate zone (ET). Despite the coastal location, the climate is strongly continental with long, cold winters (September–May) and short, cool summers. In January, the mean air temperature is −30.0°C and rises to +8.5°C in July (1987–2016). The mean annual precipitation does not exceed 350 mm (Tiksi Hydrometeorological Observatory, WMO 218240). The open-water season in the southern Laptev Sea starts on average in mid-June and extends to mid-October (Lantuit et al., 2011).

During drilling campaigns within the Russian-German science cooperation Lena Expedition in 2017, deep sediment cores from two thermokarst lagoons (UoL and PFL) and one thermokarst lake (Goltsovoye Lake (GoL)) were drilled on the Bykovsky Peninsula (Figure 1D). Furthermore, sediment characteristics of a 6.2 m deep core from the Ivashkina Lagoon (IvL) drilled in 1999 were published by Schirrmeister et al. (2018). The results showed several stages of landscape development from the initial Yedoma Ice Complex through thaw-induced lake development and eventually lagoon formation. In this study, we report on results from the approximately 30 m deep cores from Polar Fox and Uomullyakh lagoons.

The Uomullyakh Lagoon (71°43′51″N, 129°16′21″E, Figure 1D) is a 2.1 km² large, shallow lagoon (water depth 0.8 m at the borehole location) with underlying former Yedoma and Alas permafrost (Strauss et al., 2018). The lagoon is well connected with the Tiksi Bay via a 90 m wide opening in the center of the narrow and flat sand spit, which can be flooded by storm surges. In April 2017, test hole drilling (17 locations) with a Kovacs 5 cm diameter ice auger revealed bedfast ice at all locations. However, bathymetry surveys in summer 2017 (Strauss et al., 2018) suggest that small isolated water pockets might have existed in April 2017. We recorded 40 cm of compacted snow on top of 80 cm ice, which was overlying the frozen lagoon bed at the drilling location (71°43′51″N, 129°16′29″E). In the west, the lagoon is connected to two drained thermokarst lakes basins by a former drainage channel.

The Polar Fox Lagoon (71°44′35″N, 129°20′16″E; Figure 1D) is a smaller sized (0.6 km²), brackish to salty water body (water depth 3.30 m in the lagoon center) formed after partial lake drainage. This nearly closed lagoon is located in a partially drained thermokarst basin (Alas in Figure 1D) and is connected to the sea only during the open water season via an 800 m long and roughly 50 m wide outlet channel. The water level in the lagoon is at about sea level and water flow occurs in both directions depending on wind, wave, and tide conditions. Ice thickness measurements from April 2017 and bathymetry from July 2017 suggest that the lagoon was 75% covered with bedfast ice (Angelopoulos et al., 2020a). In the center of the lagoon, at the borehole location (71°44′58″N, 129°20′30″E), a 175 cm salty water layer was present underneath the 155 cm thick seasonal ice layer in April 2017 (Strauss et al., 2018).

Both lagoons are strongly impacted by warm freshwater discharge from the Lena River into the Buor Khaya and Tiksi bays (Juhls et al., 2019). This discharge causes large seasonal and interannual variations in temperature and salinity (Supplementary Table S1). For example, Polar Fox Lagoon had a salinity of 13 psu (practical salinity unit) and 44 psu below the ice cover in April 2017 and April 2019, respectively. In July 2017, the salinity dropped to less than 1 psu.

Fieldwork on the Bykovsky Peninsula was carried out in April 2017 (Strauss et al., 2018). Over eight days, a 33.5 m deep borehole was drilled from the ice surface on UoL (core number PG2410) and a 31 m deep borehole was drilled from the ice surface of neighboring PFL (PG2411) (Figure 1) using a URB2-4T drilling rig mounted on a tracked vehicle. The sediment core retrieved at UoL was 32.3 m long and at PFL 27.7 m long. Immediately after recovery, the cryolithology of the sediment cores was visually described and the cores were packed and transported frozen to AWI Potsdam. GeoPrecision temperature chains were installed in both boreholes, each with an accuracy of ± 0.1°C at 0°C and a resolution of 0.01°C. The thermistors had been calibrated in a MilliQ water-ice bath to measure their offset from 0°C prior to deployment. The chains remained in the boreholes for 5 days at Polar Fox Lagoon and 11 days at Uomullyakh Lagoon. To account for drilling heat effects, the observed temperature drift for each thermistor was analyzed following Lachenbruch et al. (1982) to estimate the undisturbed temperature. Both sediment cores were mostly unfrozen but had some partially frozen layers. The heat generated during the drilling of the cores caused the cores to partially thaw, so that the on-site haptic observation of the frozen state included some uncertainty. The PFL raw temperature data is available from Angelopoulos et al. (2020a) and the UoL temperature data and processing are described in the supporting information (Supplementary Figure A4).

In the cold lab (−9°C), the sediment cores were opened, cleaned and split, and the cryolithology was described in more detail. Furthermore, the core halves were photographed and subsamples were taken approximately every 50 cm for laboratory analyses. For the Polar Fox core, a higher sampling resolution of about every 10 cm was used for pore water extraction.

Ground temperature was recorded for 11 days in the Uomullyakh borehole and for 5 days in the Polar Fox borehole in April 2017. The freezing point depression ΔT_f of the sediment was calculated from the molality of dissolved ions estimated from the measured pore water electrical conductivity according to Eq. 2 (Atkins, 2018) and then compared with the measured borehole temperatures. Δ T f = R ( T m , f w ) 2 L f c (2)

This equation describes the lowering of the freezing point of free water with the salt content c (mol/kg). R is the universal gas constant 8.314 J/(K mol), T_m,fw is the freezing point of free water (273.15 K) and L_f is the latent heat of fusion of water (334,000 J/kg). For the conversion of the measured EC to molality (mol/kg) the MATLAB implementation of TEOS-10 (McDougall and Barker 2011) was used. It has been demonstrated that TEOS-10 also performs well for hypersaline solutions (Pawlowicz, 2012). Further details can be found in the supplementary material.

Along the coast of the Bykovsky Peninsula six lagoons exist, five of them originated from thermokarst basins (Jenrich, 2020). All of them are located in the southern part of the peninsula (Figure 1). We used a simple geomorphological approach to further differentiate them based on their connectivity with the ocean into open, semi-closed and nearly closed thermokarst lagoons (Table 1).

The extent of each of the five lagoons was determined using the Global Surface Water dataset which is based on Landsat-5, -7, and -8 satellite images from 1984 to 2018 at 30 m resolution (Pekel et al., 2016). Water bodies were defined by a water occurrence threshold of >75% over this time period. The raster dataset was vectorized and smaller geometric errors, which occurred during vectorization, were solved with the Fix Geometry function in QGIS3.6. The five Bykovsky Peninsula lagoon polygons were selected manually and these water bodies were split from the ocean by using the function “split by line”. The resulting lagoon polygons were re-projected in UTM Zone 52N (EPSG:32652) to calculate the polygon area in km².

The organic carbon pools of Uomullyakh and Polar Fox lagoons were estimated after Jongejans et al. (2018) using Equation 3, based on a deposit thickness of 30 m, the lagoon size, bulk density, and TOC content. We assume that there is no form of massive ice in the thermokarst lagoon sediments, so that the TOC pool can be calculated as: T O C p o o l ( M t ) =   t h i c k n e s s   ⋅ c o v e r a g e   ⋅ B D   ⋅   T O C 100 10 6 (3) with thickness and coverage in m, BD in kg/m³ and TOC in wt%. For further calculations, we assumed a TOC of 0.05 wt% for the samples with TOC below the detection limit (n = 21). To compensate for data gaps caused by core loss or larger sampling distances, the BD and TOC values were first multiplied, and the depth interval given as the difference of two adjoining depths was calculated. For core loss sections, this interval was accordingly higher. The product of BD and TOC was weighted and replicated for each depth interval. To interpolate the carbon content of point measurements to greater areas, the data was resampled using the bootstrapping method after Jongejans and Strauss (2020). Assuming that the deposition of organic material in the study area did not show significant differences, the calculated C-budgets of UoL and PFL were averaged and upscaled on the three further lagoons on the peninsula. Therefore, the size of each lagoon was multiplied with the mean of the TOC pool of UoL and PFL, and divided by the mean size of UoL and PFL.

We used the R environment to perform various statistical analyses on 10 parameters (R Core Team (2016); R version 3.6.3). In order to compare the differently scaled parameters, the data was normalized with the function "normalize". To subdivide the cores into stratigraphic zones, the stratigraphically constrained, incremental sum of squares cluster analysis (CONISS) was applied. Thereby, only stratigraphically adjacent clusters are considered for merging. The CONISS clustering was performed in R with the 'chclust' function from the package 'rioja'. The used dissimilarity index was 'bray'. In order to compare the properties at the time of sediment deposition, only sedimentological and biogeochemical parameters were included in the calculation. In total 10 parameters: mean depth, grain size (mean, volumetric proportion of clay, silt and sand), mass specific magnetic susceptibility, TOC, TN, TOC/TN and δ¹³C. Hydrochemical parameters were not used as pore water in unfrozen soil is flowing and mixing constantly and so does not represent the hydrochemical composition at the time of sediment deposition. The results were visualized in a dendrogram. The cutting point was selected based on visual validation and dating results.

A principal component analysis (PCA) was performed on the same parameters to reduce this high-dimensional data set to fewer dimensions and to detect patterns in the distribution of sedimentological and biogeochemical parameters over depth. The results were visualized in a biplot. The arrows show the loading on each variable, with the length approximating the variance and the angle between the arrows their correlation. Points close together correspond to depth that have similar scores on the PCA components. PCA was performed using the ‘prcomp' function from the Package ‘stats'.

At UoL, from core bottom up to 1480 cm depth (UoL-UIII), the grain-size of the sediments alternated several times from silty fine sand to coarse-grained sand. Rounded pebbles up to 2 cm in diameter were found up to 1030 cm depth. In UoL-UII, layers of silty fine sand became thicker and coarse-grained layers thinner. From 820 cm to the top (UoL-UI), the sediment consisted of black to dark grey clayish to silty fine sand with parallel sediment structures. UoL-UIb was partially frozen. From 2100 cm downwards the sediment was rich in structureless ground-ice. Between the frozen parts the core was thawed. Several organic-bearing layers were encountered along the core. At 1530 cm, a layer of large wood pieces (up to 5 cm in size) was observed.

At PFL, the sediment gradually became finer upwards in the core. From core end up to 860 cm (PFL-UIII & II), the sediment was sandy with fine, medium and coarse sand alternating frequently. Rounded pebbles up to 3 cm in diameter were found up to 1120 cm depth. Above 860 cm in PFL-Unit I, the sediment was dark grey and increasingly silty. From 630 cm onwards, shell and plant remains were encountered. Between 2245 cm and 2100 cm two 5–10 cm thick macro organic containing layers were observed. At 1810 cm, a 5 cm layer of large wood pieces (up to 8 cm in size) similar to UoL was found. The upper 1500 cm of sediment contained only sparse organic remains. The core was mostly thawed and had only small frozen sections (1815–1855 cm, 2040–2050 cm and 2970–2985 cm, and ground ice structures between 810 and 1160 cm).

Radiocarbon dates for the Uomullyakh Lagoon core ranged from >53.2 cal ka BP to 2.6 cal ka BP (Figure 2, Supplementary Table S3) and showed several age inversions from 1036 cm depth downwards. The upper two samples until 478 cm depth were of early Holocene age. Both the oldest sample (>53 cal ka BP) at a depth of 2049 cm, and the (even deeper) youngest sample (2.6 cal ka BP) at 3033 cm, were from plant macrofossils. Bulk organic samples also showed inversions although not so extreme, and varied between 18 and 34 cal ka BP. The luminescence ages of quartz (OSL) and feldspar (pIRIR) (Supplementary Table S4) were similar in range to the radiocarbon dates at 1107 cm [37 ± 2 ka (OSL), 29.9 ± 1.6 ka (pIRIR) and 43.0 ± 0.6 ka cal BP (radiocarbon)] and 2108 cm depth [63 ± 3 ka (OSL), 59 ± 4 ka (pIRIR) and >53.2 ka cal BP (radiocarbon)]. Since quartz dose recovery was satisfactory for these samples from UoL, the OSL ages are considered more reliable for these samples than the feldspar and radiocarbon ages. The lowermost luminescence age at 2978 cm was non-finite for both quartz (>100 ka) and feldspar (>300 ka). Since feldspar provided the more limiting constraint, the minimum feldspar age is shown in Figure 2.

Radiocarbon dates for Polar Fox Lagoon ranged from > 53.3 cal ka BP to 33.5 cal ka BP (Figure 3, Supplementary Table S3). For the uppermost sample, the radiocarbon dating failed due to insufficient organic carbon content, so that available dates started at 480 cm sediment depth (corresponding to depth of 813 cm measured from ice surface). As seen at UoL, the age dating revealed no linear trend over depth. The youngest age was found in the greatest depth (2761 cm), while the oldest sample was found at 1130 cm sediment depth. The luminescence ages included only non-finite ages for quartz and feldspar of >100 ka and >300 ka at 2771 cm and of >100 ka and >320 ka at 3002 cm depth (Supplementary Figure S1, Table S4). Therefore, the feldspar pIRIR ages were used for interpretation of the lowest core ages, as with UoL.

The radiocarbon ages of the lowest unit UIII of both cores were not used in the discussion because of low carbon content and high risk of contamination. The results of the luminescence dating, aiming at a depositional age of mineral grains, were considered more reliable. The two samples with quartz OSL ages from the Uomullyakh Lagoon core had ages well within the range of the method and they were in stratigraphic order.

We acknowledge the fact that shore erosion and sediment mixing in such permafrost thaw-affected aquatic environments may easily result in radiocarbon age inversions and thus challenge chronological interpretation.

Detailed dose rate results and luminescence ages are summarized in Supplementary Table S4. Information on radionuclide concentrations are shown in Supplementary Table S5.

Results of the hydrochemical analyses of pore water are shown in Figure 4 and also in Supplementary Table S8. The pore water pH of both lagoons was neutral to slightly alkaline with a median of 7.5 for UoL and 8.0 for PFL. The EC was much higher in the pore water of the UoL-core with a maximum of 108 mS/cm at 1955 cm depth (median: 83.1 mS/cm) and showed higher variations over depth. Relatively low EC values at 2267 and 2745 cm depth were correlated with heavy stable water isotopes (Supplementary Figure S2A). The EC of the Polar Fox core pore water was highest in the upper sediment layer (max: 41 mS/cm) and showed freshwater to brackish conditions from 800 cm depth downwards (median: 3.3 mS/cm). The DOC distribution differed between the two lagoons, while the median value was similar (UoL: 58 mg/L, PFL: 50 mg/L). In UoL, the DOC concentration varied most in the uppermost unit and was highest at 360 cm depth (282 mg/L), while in PFL the greatest variation was found in the lower 5 m with the highest DOC concentration measured at 2835 cm depth (221 mg/L; Figure 4 and Supplementary Table S8).

The stable water isotope composition was generally heavier in the pore water of the semi-closed Uomullyakh Lagoon (mean δ¹⁸O/δD: −20.9‰/−163.4‰) (Supplementary Figure S2A) than in the nearly closed Polar Fox Lagoon (mean δ¹⁸O/δD: −17.7‰/−133.7‰) (Supplementary Figure S3A). In the UoL core, the pore water of the upper sample had the lightest isotope composition with δ¹⁸O of −19.8‰, which is comparable to the upper 5 m of PFL core pore water where its isotope composition was heaviest. There was no significant difference in isotopic composition between frozen and unfrozen samples (Supplementary Figure S2B, S3B). For UoL, we found a depletion of δ¹⁸O with increased salinity. This was not found in the mainly freshwater to brackish dominated PFL core.

High freezing point depressions were observed for hypersaline samples (Figure 4) (e.g. for Uomullyakh at a depth of 1955 cm). For UoL, the borehole temperature was below the calculated freezing points in the lower core (1880 cm sediment depth and deeper) and thus the sediment could have been frozen at the time of drilling. In comparison to Uomullyakh, the Polar Fox borehole temperatures were significantly warmer (up to 5.2°C warmer) and significantly exceeded the calculated freezing point in the upper 480 cm of sediment above the partially frozen layer. Within Unit II and III of PFL, the calculated freezing point was very close to the observed temperature, suggesting that the system was close to equilibrium. This confirms the observations made in the field on the mostly thawed state of the core.

The surface areas for the Uomullyakh, Polar Fox, Ivashkina, Stepanenko, and Safroneeva lagoons were determined with remote sensing to be 2.10, 0.64, 3.87, 2.88, and 1.89 km², respectively. The total thermokarst lagoon area is 11.4 km². Table 2 presents C-budgets and C-density calculated for the upper 3 m and 30 m of the sediment deposited in the Uomullyakh and Polar Fox lagoons. Results were extrapolated for the deposits below the other three thermokarst lagoons on the southern Bykovsky Peninsula using mean values for C-budget and C-density of Uomullyakh and Polar Fox lagoons. Below UoL, a carbon stock of 0.15 ± 0.01 Mt with a carbon density of 24.37 ± 2.09 kg/m³ for the upper 3 m and 1.14 ± 0.2 Mt with a carbon density of 18.03 ± 3.12 kg/m³ for 30 m thick deposits was calculated. The upper 3 m of PFL stored 0.05 ± 0.001 Mt organic carbon with a density of 26.50 ± 0.57 kg/m³. The carbon stock amounts to 0.24 ± 0.04 Mt with a carbon density of 12.55 ± 0.04 kg/m³ for 30 m thick deposits. The total carbon pool of the thermokarst lagoons on the southern Bykovsky Peninsula is 5.72 ± 0.63 Mt with an average carbon density of 15.29 ± 1.94 kg/m³ for the upper 30 m (for the upper 3 m: 0.85 ± 0.09 Mt, 25.44 ± 0.75 kg/m³).

The first two PCA axes explained 88% of the variation in the data for the Uomullyakh Lagoon core and 82% for PF Lagoon core (Figure 5). Though the variance is small, two groups were distinguished for both lagoon cores. Shallow depth correlated with fine deposits and high organic matter content and quality, while in the deeper core, the variables sand content and δ¹³C were decisive. The main variance in the data was explained by differences in grain size and organic matter contents and composition along the core.

At both lagoons, we identified four phases of landscape development based on sedimentological and biogeochemical parameters. These are: 1. Middle/Late Pleistocene fluvial deposition, 2. Late Weichelian (MIS 2) Yedoma deposition, 3. Holocene thermokarst lake formation and 4. Holocene lagoon formation.

The grain size distribution of sediments beneath Uomullyakh Lagoon and Polar Fox Lagoon indicates a changing depositional environment over time. In both lagoons, the sediments represent accumulation over the last glacial and interglacial cycles (Table 3). It is very likely that only the top few meters were accumulated since the previous interglacial. Based on stratigraphically constrained cluster analysis (CONISS) and visual observation, we distinguished four main stages of the late Quaternary history of the study area. The horizontal axis of the dendrogram (heights) represents the distance and therefore dissimilarity between clusters. The CONISS analysis of the 10 sedimentological and biogeochemical parameters at UoL showed a greater dissimilarity between clusters as for PFL. That hints to a more diverse depositional history for UoL. The cut of the dendrogram at height 7 for UoL (Figure 2) and height 4 for PFL (Figure 3), divided the cores into three main units, which is supported by the results of the PCA (Figure 5).

Radiocarbon ages out of stratigraphic order due to reworking of sediment as part of thawing and re-freezing processes are common for sedimentary deposits in thermokarst landscapes (Wetterich et al., 2009; Biskaborn et al., 2013; Schirrmeister et al., 2017, 2018; Jongejans et al., 2018). In the case of UoL, radiocarbon age reversals occur in Unit II (11–17 m) and Unit III (below 17 m). Possible reasons for unexpected high radiocarbon concentrations include processes associated with the melt of large ice wedges and subsequent mixing, which transport younger material into older layers. Also, possible contamination by younger material during the drilling process and sample preparation have a larger effect when sample TOC contents are very low and available material is limited and can lead to dating bias (Oswald et al., 2005; Vyse et al., 2020). Unit II is characterized by thawed Yedoma Ice Complex deposits (see chapter 5.1.2). Existing intact Yedoma deposits, in north-east Siberian Arctic lowlands consist of massive syngenetic ice wedges several tens meters in height (Schirrmeister et al., 2002, 2011, 2013; Wetterich et al., 2014; Strauss et al., 2017; Möller et al., 2019; Fuchs et al., 2020). On the Bykovsky Peninsula only 10 km northeast of the lagoons, the Mamontovy Khayata Yedoma exposure is an example with syngenetic ice wedges up to 40 m high and 5–6 m wide (Schirrmeister et al., 2002). Likely, the same deposits were present at our drilling sites before thermokarst degradation. To estimate total thaw subsidence for the sites we measured the elevation difference of lagoon bottoms and the surrounding Yedoma upland based on height information of the ArcticDEM (Porter et al., 2018). The results suggest that Holocene permafrost degradation subsided the Yedoma surface by up to 25 m for UoL and 27 m for PFL, potentially causing mass flux of sediment and its organic carbon. At many locations across the landscape, large ice wedges thawed and the surface subsided, leaving ice wedge casts and water bodies, both of which infilled with local sediment. Where young near-surface organic matter slumped, and mixed into lower, older layers, age inversions as seen at UoL would have been created.

Though sediment mixing in such thaw-affected aquatic environments challenges chronological interpretation, comparison with deposits at other sites on the peninsula are consistent with our interpretation of the landscape development. In general, the radiocarbon age ranges of both cores correspond to previous dating from Yedoma sequences of Mamontovy Khayata (58.4–12 ka BP), about 10 km northeast (Schirrmeister et al., 2002), of Cape Razdelny (core 103/81: 40.8–21.63 ka BP) about 10 km southwest (Slagoda, 1993; Grosse et al., 2016) and of the Goltsovoye Lake (47.5–21.39 cal ka BP) less than 500 m north (Jongejans et al., 2020). Furthermore, the luminescence dating at 1107 cm and 2106 cm depth is consistent with the corresponding radiocarbon ages, showing an overall chronological succession. This sequence comprises Pleistocene ages in deep fluvial deposits, followed by late Weichelian ages in overlying Yedoma deposits. Therefore, the age range for the upper sections (UII, UI) appear reasonable. We correlate stratigraphic units between cores to assign the late Pleistocene ages of UoL-UII to the upper part of PFL-UII.

For deeper layers, especially Unit III of both cores, the results of the luminescence dating were used. The dose recovery ratios are satisfactory, and the comparison with the more-difficult-to-bleach feldspar IR₅₀ and post-IRIR signals indicates that it is very likely that the quartz was well bleached at deposition (see Supplementary Material Chapter 1.3 for a more detailed technical discussion of these issues).

Unfortunately, the quartz luminescence signals from the bottom Uomullyakh Lagoon sample, and both samples from the Polar Fox Lagoon, are saturated, and so all give minimum quartz age estimates of >100 ka. Interestingly, the feldspar pIRIR signals are also saturated, and since the dose range of the feldspar signal is considerably larger than that for the quartz OSL signal, the implied minimum age is > 300 ka. At least for the two lagoons located in a similar morphological setting and with similar sedimentology (see below) it is reasonable to assume that the bleaching history of the sediment entering the two water bodies is also similar. Luminescence dating from the adjoining Goltsovoye Lake supports this assumption (Supplementary Table S4). The quartz and feldspar signals of the bottom sample (fluvial deposits at 3513 cm depth) are also saturated and in the same age range as for UoL and PFL (quartz age: >90 ka, feldspar age: >280 ka). The implication of the broad agreement between the pIRIR and OSL signals in the top two Uomullyakh samples is that both were well bleached before deposition. If we extrapolate this observation to the bottom-most Uomullyakh sample and the two Polar Fox samples, we may assume that the pIRIR minimum ages should also be reliable, and since they are the most constraining, they are used in further discussion.

Our data suggests that deep (∼20 m) hypersaline taliks (UoL) preserve organic carbon better than taliks that contain mostly freshwater below a thin saline layer (PFL). Drawing on microbial subsea permafrost research (Mitzscherling et al., 2017), hypersaline cryotic taliks can stress microbes and slow the decay of organic matter.

Once a connection from the lake to the sea is established, salts from seawater can diffuse into the sediment (Angelopoulos et al., 2020a) or transport into the sediment more rapidly via convection (Shakhova et al., 2017). Therefore, the frozen/thawed state of the sediment is largely a reflection of recent marine sediment deposition and salt transport. In the Uomullyakh Lagoon core, the EC of pore water generally increased with depth, from 30 mS/cm just below the seasonally frozen layer (240 cm below ice) to 108 mS/cm at 1525 cm depth below ice (mean 77 mS/cm for the entire core). The salinities were significantly higher than the Tiksi Bay water (3.8 mS/cm) (Supplementary Table S1) for the same time of the year. At the Polar Fox Lagoon core, the maximum EC (41 mS/cm) was measured in the uppermost sample just below the sediment surface. The relatively high EC in both lagoons suggests that brine rejection from lagoon ice formation and the seasonal isolation of the water bodies from Tiksi Bay are key processes affecting the pore water salinity.

Uomullyakh is a semi-closed lagoon, especially during the open water season when it may experience flooding and enhanced marine sediment deposition from storm surges (e.g. Nazarova et al., 2017) (Figure 6A). Polar Fox Lagoon, on the other hand, is a nearly closed lagoon (Figure 6B). It is isolated from the sea once bedfast ice forms at the shallow inlet (∼0.5 m) of the connecting channel (Spangenberg et al., 2021).

At Uomullyakh Lagoon, all 17 ice-auger measurements in April 2017 revealed bedfast ice, suggesting that there was nearly no surface water exchange with Tiksi Bay. Despite a seasonally closed system for the surface water, lateral groundwater exchange in the sediment can still occur if the sub-lagoon and subsea permafrost taliks are connected. This talik connection is possible, because the sediment just offshore of the spit is also characteristic of the lagoon talik when the spit was located further south. Furthermore, permafrost aggradation beneath the spit is controlled by limited sub-aerial exposure time as the spit retreats further north with coastal erosion. The large difference in isotopic composition of the pore water between the lagoons confirms that the Uomullyakh Lagoon sediment had a larger source water contribution from Tiksi Bay (Figure 7). Tiksi Bay as a heavy endmember likely caused the enrichment in heavier water isotopes of the Uomullyakh Lagoon pore water. Pedrazas et al. (2020) used electrical resistivity surveys to map continuous unfrozen sediment from an onshore cryopeg to a sub-lagoon talik down to a depth of 20 m in northern Alaska, and suggested it could be a conduit for groundwater and nutrient exchange. The isotopic composition of the pore water, along with the near vertical profile of salinity with depth (5–20 m and 25–33 m; Supplementary Figure S3) suggests that convective processes influenced the salt distribution (Harrison and Osterkamp, 1982). In addition, there was a likely high salt content injection into the sediment following bedfast ice formation. The seasonally frozen layer can still have been permeable and susceptible to brine flow (Osterkamp et al., 1989). Although temperatures were as low as −4.4°C (at 1540 cm), the freezing point of the hypersaline sediment was calculated to have been as low as −5.1°C (at 1520 cm), resulting in a cryotic talik. At a depth of 2040 cm, the observed temperature (−4.7°C) was below the freezing point (−3.82°C), indicating the approximate depth to ice-bearing permafrost. It is unclear how the ice-bearing permafrost table evolved, but the decreasing temperature trend with depth suggests the lagoon sediment column might be warming. Therefore, it is plausible that hypersaline sediment partially froze once sufficient sediment deposition in the lagoon created a mostly bedfast ice regime with some small isolated water pockets. Now, the ice-bearing permafrost table may be degrading under warmer, but still cryotic bedfast ice conditions. This is similar to terrestrial cryopeg expansion under warming subaerial conditions (Streletskaya, 1998).

At Polar Fox Lagoon, the contemporary distribution of ice-bearing permafrost and taliks is somewhat different. Polar Fox Lagoon is presumably younger than Uomullyakh Lagoon and still contains a floating ice area with a maximum water depth of 340 cm. In the floating ice area, coupled heat and salt diffusion models suggest that the thin partially frozen permafrost layer (830 to 480 cm below the sediment surface) developed sub-aquatically after the lake to lagoon transition. The partially frozen layer depth and thickness are transient, as the top-down chemical degradation of newly formed frozen sediment lags behind the talik refreezing (Angelopoulos et al., 2020a). In bedfast ice zones (approximately 75% of Polar Fox Lagoon), electrical resistivity surveys showed that the thickness of the frozen layer and that of the overlying saline layer increased and decreased toward the shoreline, respectively.

The transient nature of frozen sediment and taliks beneath may affect CO₂ and CH₄ fluxes in the sediment column, because partially frozen sediment with an ice saturation exceeding 60% has proven to be an effective barrier to gas flow (Chuvilin et al., 2016). In the case of a shallow lagoon with bedfast ice, the frozen sediment can be saline and have a lower freezing point compared to frozen sediment without salt. Therefore, lagoons may precondition the frozen sediment for rapid thawing (compared to Yedoma permafrost) once it is exposed to warmer coastal waters in the subsea permafrost environment. This, in turn, could contribute to the development of offshore gas migration pathways for greenhouse gas released from thawing subsea permafrost (Shakhova et al., 2019). Any released gas might not reach the atmosphere, as it has to bypass numerous sinks, including anaerobic oxidation in the sediment column (Overduin et al., 2015), dissolution in the pore water and surface water, low-permeability layers, and lagoon ice. Angelopoulos et al. (2020b) provide a review of subsea permafrost within the context of global climate change, and Ruppel and Kessler (2017) provide a review specific to climate change interactions with gas hydrates.

Permafrost thaw beneath thermokarst lakes is more rapid than gradual thaw of near-surface permafrost and could double permafrost carbon emissions this century (Walter Anthony et al., 2018; Turetsky et al., 2020). Field-based observations on carbon turnover and greenhouse gas production in thermokarst lake sediments are scarce (Heslop et al., 2015; Heslop et al., 2019a; Heslop et al., 2019b) and completely missing for thermokarst lagoons. In order to estimate potential future greenhouse gas release from emerging thermokarst lagoons, we assessed the size of the lagoon C-pool by extrapolating our TOC and BD of the two analyzed lagoons to the other three thermokarst lagoons on the Bykovsky Peninsula.

The calculations of the carbon budgets show that the carbon density differs between the neighboring lagoons, especially for greater depths (upper 3 m: 22.5 kg/m³ (Uo) vs. 26.5 kg/m³ (PFL); upper 30 m: 17.4 kg/m³ (UoL) vs. 10.2 kg/m³ (PFL), which demonstrates the small scale spatial variability of the carbon deposits. The input and output of organic matter are strongly dependent on environmental conditions and lake drainage events. Schirrmeister et al. (2018) estimated the carbon density for the upper 6 m of the thermokarst lagoons on the Bykovsky Peninsula at 25 ± 5 kg/m³, which is in the same range as the carbon density we calculated for the upper 3 m (24.5 ± 1.4 kg/m³), but 45% higher than the mean carbon density of the 30 m-long sediment cores (13.8 ± 3 kg/m³). We explain this with the higher TOC values we measured in the lake influenced sediments (UoL-UIa, PFL-UI), which were the stratigraphic class dominating the Schirrmeister et al. (2018) dataset.

A similar carbon storage was estimated for thermokarst lagoons on northern Richards Island, Mackenzie Delta, Canada. Based on an average carbon content of 2.75–3% for marine and lacustrine deposits, Solomon et al. (2000) estimated a carbon deposition of ca. 25 × 10⁴ t in the investigated embayment area (10.3 km²) since thermokarst lake development. This corresponds to a carbon density of 25 kg/m³ for the upper 10 m, which is similar to the carbon density of the lacustrine sediments of the Bykovsky lagoons. Comparing the carbon stored in undisturbed, previously unthawed ice-rich Yedoma deposits of Mamotovy Khayata (Schirrmeister et al., 2011) with the TOC of Yedoma thawed under thermokarst lake and lagoon conditions (Table 4), we would expect to find a gradient of C loss in thawed Yedoma underneath thermokarst landforms that is linked to thaw histories. Using a simple mass balance approach, the difference in mean TOC content between the Mamotovy Khayata Yedoma and the thawed Yedoma beneath thermokarst sites (Table 4) can be expressed as percentage carbon loss. This approach suggests a C loss of 85% for thawed Yedoma under Goltsovoye Lake, 83% C loss for PFL, 48% C loss for UoL, and 44% C loss for IvL due to subaquatic permafrost thaw. These estimations imply that less organic carbon is lost in the more open UoL and IvL compared to the younger, nearly closed PFL and the freshwater GoL. If our assumptions are valid, this would signify that a freshwater dominated lake and the floating ice areas of nearly closed lagoon systems like PFL result in strong organic carbon loss upon thaw while seawater dominated systems seem to better preserve organic carbon. In addition, higher TOC/TN and lower δ¹³C values for thawed Yedoma underneath UoL and IvL indicate lower organic matter degradation both more open lagoons, further supporting this hypothesis. Since approximately 75% of PFL was characterized by bedfast ice in April 2017 (Angelopoulos et al., 2020a), additional C is likely preserved in areas with thick frozen sediment. At UoL, the hypersalinity of the cryotic talik can stress microbial communities, leading to reduced population, diversity, and decay of organic matter (Mitzscherling et al., 2017). Although bacterial communities can adapt to increasingly saline conditions, past refreezing events in the upper 20 m of the sediment core (refer to “Impact of Marine Inundation” Section 5.2), can also preserve C at UoL.

The calculated average C-density of the core section referred to Yedoma deposits of the five thermokarst lagoons on Bykovsky Peninsula is 24.72 ± 4 kg/m³ and therefore in the upper range of the C-density estimated for the 0.41 × 10⁶ km² large Yedoma region (19 + 13/-11 kg/m³) (Strauss et al., 2013) and in the lower range of terrestrial thermokarst deposits of the Yedoma region (33 + 25/-19 kg/m³) (Schirrmeister et al., 2011; Strauss et al., 2013).

Using the lagoon C-density mean, we calculated 5.7 Mt C in the first 30 m of lagoon sediments covering an area of 11.4 km². This is a substantial inventory of formerly frozen C, which is now largely unfrozen and thus available for microbial degradation processes. δ¹³C and TOC/TN values indicate that the degradation of the primarily terrestrial plant material is advanced. Highest organic matter degradation was indicated for the fluvial deposits possibly because the alluvial material was already strongly degraded at the time of deposition.

The comparison of different sites shows that the biogeochemical parameters differ between sediment types rather than between the sites (Table 4). For example, the lacustrine sediments at all four sites (PFL, UoL, GoL and IvL) have very similar TOC and TN contents as well as similar TOC/TN ratios and δ¹³C signals, but they differ from the over- or underlying sediment facies (Schirrmeister et al., 2018; Jongejans et al., 2020). The lacustrine sediments in the deep sediment cores of UoL, PFL and GoL are richer in TOC and TN than the underlying former terrestrial or overlying marine sediments. High carbon accumulation in thermokarst lake sediments can be caused by shore erosion and deposition of terrestrial organic matter, high aquatic productivity, and unique preservation conditions (Anthony et al., 2014). In contrast, the biogeochemical parameters of lacustrine and underlying Yedoma deposits of the Ivashkina Lagoon differ only slightly (Schirrmeister et al., 2018). This is not surprising, as the upper part of the Yedoma deposits (IvL core reaches only to 6.15 m depth) tend to have higher carbon and nitrogen contents and lower degradation rate than the deeper deposits. An exception are the deep fluvial deposits under the GoL, where allochthonous organic material leads to very high TOC and TN contents (Jongejans et al., 2020). The TOC content of thawed talik sediments (former Yedoma deposits) is highest for the UoL, with similar values found in Taberites (2.7 ± 1.4 wt%) and late Weichelian (MIS 2) Ice Complex deposits (2.2 ± 0.9 wt%) of the Siberian permafrost region (Schirrmeister et al., 2011). The organic matter mineralization is highest for former Yedoma deposits below GoL (TOC/TN = 4) and lowest for Yedoma deposits at Mamontovy Khayata (TOC/TN = 20). The high mean TOC/TN ratio of fluvial deposits below GoL is based on five samples with very high TOC/TN ratios. The median TOC/TN ratio of 0 shows that the carbon mineralization is high for most samples of the fluvial deposits below GoL. In general, the carbon mineralization is mostly high (TOC/TN < 12) or in the lower moderate range (TOC/TN 12–25) for thermokarst affected sites. This is a common phenomenon in the warming Arctic, as rising permafrost temperatures result in increasing thermokarst development which favors the microbial decomposition of thawing organic matter (Strauss et al., 2017).

Greenhouse gas production in thermokarst lagoons connected to the marine environment is expected to differ substantially from that in thermokarst lakes: the supply of electron acceptors, specifically of sulfate, will promote sulfate reduction and anaerobic methane oxidation thus favoring the release of CO₂ relative to that of methane (Segers and Kengen, 1998). In line with this, substantial CO₂ release was recently reported for a laboratory study that incubated permafrost with seawater (Tanski et al., 2019). In addition, anaerobic methane oxidation was observed in thawing seawater-affected permafrost (Winkel et al., 2018) and thermokarst lake sediments with temperatures close to 0°C (Winkel et al., 2019). Given both our first estimate of thermokarst lagoon carbon budgets, and the expectation that the rate of formation of such lagoons will increase due to coastal erosion, thermokarst lagoons are considered important biogeochemical processors of former permafrost carbon and a critical component of the future Arctic greenhouse gas system.