Thermokarst Lagoons: A Core-Based Assessment of Depositional Characteristics and an Estimate of Carbon Pools on the Bykovsky Peninsula

Permafrost region subsurface organic carbon (OC) pools are a major component of the terrestrial carbon cycle and vulnerable to a warming climate. Thermokarst lagoons are an important transition stage with complex depositional histories during which permafrost and lacustrine carbon pools are transformed along eroding Arctic coasts. The effects of temperature and salinity changes during thermokarst lake to lagoon transitions on thaw history and lagoon deposits are understudied. We analyzed two 30-m-long sediment cores from two thermokarst lagoons on the Bykovsky Peninsula, Northeast Siberia, using sedimentological, geochronological, hydrochemical, and biogeochemical techniques. Using remote sensing we distinguished between a semi-closed and a nearly closed lagoon. We (1) characterized the depositional history, (2) studied the impact of marine inundation on ice-bearing permafrost and taliks, and (3) quantified the OC pools for different stages of thermokarst lagoons. Fluvial and former Yedoma deposits were found at depth between 30 and 8.5 m, while lake and lagoon deposits formed the upper layers. The electrical conductivity of the pore water indicated hypersaline conditions for the semi-closed lagoon (max: 108 mS/cm), while fresh to brackish conditions were observed beneath a 5 m-thick surface saline layer at the nearly closed lagoon. The deposits had a mean OC content of 15 ± 2 kg/m3, with higher values in the semi-closed lagoon. Based on the cores we estimated a total OC pool of 5.7 Mt-C for the first 30 m of sediment below five mapped lagoons on the Bykovsky Peninsula. Our results suggest that paleo river branches shaped the middle Pleistocene landscape followed by late Pleistocene Yedoma permafrost accumulation and early Holocene lake development. Afterward, lake drainage, marine flooding, and bedfast ice formation caused the saline enrichment of pore water, which led to cryotic talik development. We find that the OC-pool of Arctic lagoons may comprise a substantial inventory of partially thawed and partially refrozen OC, which is available for microbial degradation processes at the Arctic terrestrial-marine interface. Climate change in the Arctic leading to sea level rise, permafrost thaw, coastal erosion, and sea ice loss may increase the rate of thermokarst lagoon formation and thus increase the importance of lagoons as biogeochemical processors of former permafrost OC.

Permafrost region subsurface organic carbon (OC) pools are a major component of the terrestrial carbon cycle and vulnerable to a warming climate. Thermokarst lagoons are an important transition stage with complex depositional histories during which permafrost and lacustrine carbon pools are transformed along eroding Arctic coasts. The effects of temperature and salinity changes during thermokarst lake to lagoon transitions on thaw history and lagoon deposits are understudied. We analyzed two 30-m-long sediment cores from two thermokarst lagoons on the Bykovsky Peninsula, Northeast Siberia, using sedimentological, geochronological, hydrochemical, and biogeochemical techniques. Using remote sensing we distinguished between a semi-closed and a nearly closed lagoon. We (1) characterized the depositional history, (2) studied the impact of marine inundation on ice-bearing permafrost and taliks, and (3) quantified the OC pools for different stages of thermokarst lagoons. Fluvial and former Yedoma deposits were found at depth between 30 and 8.5 m, while lake and lagoon deposits formed the upper layers. The electrical conductivity of the pore water indicated hypersaline conditions for the semiclosed lagoon (max: 108 mS/cm), while fresh to brackish conditions were observed beneath a 5 m-thick surface saline layer at the nearly closed lagoon. The deposits had a mean OC content of 15 ± 2 kg/m 3 , with higher values in the semi-closed lagoon. Based on the cores we estimated a total OC pool of 5.7 Mt-C for the first 30 m of sediment below five mapped lagoons on the Bykovsky Peninsula. Our results suggest that paleo river branches shaped the middle Pleistocene landscape followed by late Pleistocene Yedoma permafrost accumulation and early Holocene lake development. Afterward, lake drainage, marine flooding, and bedfast ice formation caused the saline enrichment of pore water, which led to cryotic talik development. We find that the OC-pool of Arctic lagoons may

INTRODUCTION
Over the last 2 decades, the Arctic has been warming more than twice as fast as the global average (Johannessen et al., 2004;Berner et al., 2005;Notz and Stroeve 2016). This amplified warming has led to rapid surface warming in Siberia with modeled temperature rise of up to 4°C over the last 5 decades (Romanovsky et al., 2010;Biskaborn et al., 2019;Lenssen et al., 2019;GISTEMP Team, 2020). As a consequence, permafrost is thawing in many regions, mobilizing large amounts of sediments including organic carbon and allowing microbial decomposition of previously frozen organic matter (Strauss et al., 2013;Strauss et al. 2021;Schuur et al., 2015). Arctic coastal systems are especially affected by rapid permafrost thaw and mobilization of organic matter by erosion and marine inundation (Fritz et al., 2017). Lagoons, a frontier environment between land and ocean, are transitional and dynamic coastal landforms combining characteristics of both terrestrial and marine systems (Kjerfve, 1994;Tagliapietra et al., 2009;Harris et al., 2017). Globally, coastal lagoons are defined as shallow coastal water bodies separated from the sea by a barrier, connected at least intermittently to the ocean by one or more restricted inlets, and usually have a shore-parallel orientation (Kjerfve, 1994). As lagoons are complex and diverse systems, there are many different approaches to classify them, for example based on geomorphological and geological origin, their morphological properties such as size and water depth (Kosyan, 2016), the degree of isolation from the sea (Kjerfve, 1994), their physicochemical state (for example, salinity, ionic composition, temperature, turbidity, nutrients) (Tagliapietra et al., 2009), or the degree of influence by living organisms such as corals and humans (Brovko 1990).
Permafrost thaw results in landscape change as melting excess ground ice causes surface subsidence, eventually forming thermokarst lakes (Osterkamp et al., 2009;Jones et al., 2011;Lenz et al., 2016). Along eroding sections of Arctic coasts, lowered thermokarst terrain is affected by flooding of nearshore thermokarst lakes and basins with seawater, transforming them into thermokarst lagoons (Ruz et al., 1992;Romanovskii et al., 2000). Thermokarst lagoons are also formed by natural thermokarst lake dynamics in the course of lake drainage at coasts (Arp et al., 2010). However, rates of formation are expected to intensify as a result of ongoing climate change in the Arctic. Increasing coastal erosion Jones et al., 2018), sea-level rise (Proshutinsky et al., 2001;Watson et al., 2015;Nerem et al., 2018), and higher sea water temperatures (Bindoff et al., 2007) are likely to result in an acceleration of erosion of permafrost coasts and in particular lagoon formation. Incubation experiments performed by Tanski et al. (2019) indicate that along eroding permafrost coastlines, large amounts of carbon dioxide can be produced. Therefore, thermokarst lagoon formation is an important process affecting permafrost carbon pools along rapidly changing permafrost coasts. In addition, marine inundation affects the temperature and salinity of taliks formerly beneath thermokarst lakes during their transition to the marine environment (Angelopoulos et al., 2020a). The formation of hypersaline, unfrozen ground beneath thermokarst lagoons may preserve existing thermokarst lake taliks, forming migration pathways for methane (Shakhova et al., 2019;Angelopoulos et al., 2020a). However, lagoon formation can also cause taliks to partially refreeze, provided that heat loss is faster than salt diffusion (Angelopoulos et al., 2020a). Ascending gas can get trapped below low-permeability sediment layers (Ruppel and Kessler, 2017) and be released when conditions warm up. These pathways, which may contain low-permeability ice-saturated zones, become part of the offshore environment with sufficient coastal erosion and so become potential sources of gas release in subsea permafrost areas (Frederick and Buffett, 2014;Shakhova et al., 2017).
Although thermokarst lagoons are widespread along the pan-Arctic coast (Jenrich, 2020), there is still no qualitative or quantitative assessment of their distribution or of the size of the total carbon reservoir they represent. Until now, Arctic thermokarst lagoon research has concentrated on the Mackenzie Delta region in Northwest Canada (Hill, 1990;Ruz et al., 1992;Héquette et al., 1995;Campeau et al., 2000;Solomon et al., 2000) and the Bykovsky Peninsula in Northeast Siberia. On the Bykovsky Peninsula, several lagoons have evolved from thermokarst lake basins in ice-rich Yedoma Ice Complex permafrost, and these are now actively developing in different stages (Romanovskii et al., 2000;Romanovskii et al., 2004;Ulyantsev et al., 2017;Schirrmeister et al., 2018). Still, little is known of the sedimentary history of thermokarst lagoons, the impact of marine inundation, and the amount and quality of carbon stored in greater depths.
To address these research gaps, we applied a multidisciplinary approach to characterize the sedimentological, geochronological, hydrochemical, and biogeochemical properties of two approximately 30 m deep (below lagoon ice surface level in spring) sediment cores from two typical thermokarst lagoons on the Bykovsky Peninsula ( Figure 1). We paired this with remote sensing observations to identify, map, and characterize thermokarst lagoons in a regional context along the coast of the Bykovsky Peninsula. The specific aims of this paper are to characterize (1) the sediments and (2) pore water below the lagoons and (3) to quantify the regional carbon inventory based on spatially explicit lagoon areas in the region.

Study Area
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) , 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., 2018). It lies within the continuous permafrost zone and is characterized by thick, icerich 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 . The mean annual precipitation does not exceed 350 mm (Tiksi Hydrometeorological Observatory, WMO 218240). The openwater 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 2 large, shallow lagoon (water depth 0.8 m at the borehole location) with underlying former Yedoma and Alas permafrost . 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  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 2 ), 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 .
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.

Field Work
Fieldwork on the Bykovsky Peninsula was carried out in April 2017 . 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 TABLE 1 | Classification of lagoons into "open", "semi-closed" and "nearly closed" lagoon systems which occur along a gradient of coastal erosion. Connectivity  with the sea   Impact  of coastal erosion   Examples   Open  Full and broad exchange with the sea  Fully affected and far ahead in the  transition to an open marine system Stepanenko and Safroneeva lagoons (Figure 1), many in the northern Lena Delta, and on the Cape Halkett coast of Beaufort Sea Semiclosed Barrier islands or sand spits partially block exchange with the sea either temporally or spatially Moderately to substantially affected Uomullyakh Lagoon (Figure 1), many on Tuktoyaktuk Peninsula, eastern Alaskan Beaufort Sea (Harris et al., 2017) Nearly closed Narrow outlet/inlet channel connecting lagoon and sea, with the possibility of having lake character temporarily Not yet affected Polar Fox Lagoon (Figure 1), and other lagoons connected through small channels Frontiers in Earth Science | www.frontiersin.org July 2021 | Volume 9 | Article 637899 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).

Laboratory Analyses
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.

Hydrochemical Analyses
Pore water was extracted from thawed samples using Rhizons ™ with a membrane pore size of 0.15 µm and analyzed for pH, electrical conductivity (EC, mS/cm), DOC (mg/L) and stable water isotopes (δD, δ 18 O (‰ vs. VSMOW)). The pH and EC values were measured with a WTW Multilab 540 (accuracy: ± 0.01/± 1 mV). DOC samples were acidified with 50 μL of 30% HCl supra-pure and stored at + 4°C until measured with a Shimadzu Total Organic Carbon Analyzer (TOC-VCPH) (accuracy: ± 1.5%) following Fritz et al. (2015). For analyzing the ratios of δD and δ 18 O, water and ice samples were processed following Meyer et al. (2000) and measured with a Finnigan MAT Delta-S mass spectrometer.

Sedimentology Analyses
Grain size measurements were conducted to gain information about the transportation medium and the level of sediment accumulation at the study areas. Mass specific magnetic susceptibility, describing the ability of a material to be magnetized, gives an indication of the amount of ferri-and ferromagnetic minerals in a sample. Differences in mineral composition between samples usually indicate sediment sources and changes in stratification (Liu et al., 2010). For both measurements the samples were pre-weighed and freeze-dried. The absolute ice or water content was determined based on the difference in weight before thawing and after freeze drying related to the wet weight (Phillips et al., 2015). The bulk density (BD in kg/m 3 ) was calculated from the ice content according to Strauss et al. (2012). The grain size distributions of 74 samples from both cores were measured with a Malvern Mastersizer 3000 laser particle size analyser with a measuring range of 0.01-1000 µm. Additional details on BD calculation and grain size measurement can be found in the supplementary material. Mass-specific magnetic susceptibility (MS in 10 -8 m 3 /kg) was measured with a Magnetic Susceptibility Meter (Model MS2, Sensor Type MS 2B, Bartington Instruments) on the freeze-dried, but not homogenized samples at a frequency of 0.465 kHz (low frequency).

Bulk Biogeochemistry
To gain information about changes in the nutrient cycle during the transformation of terrestrial permafrost into sub-aquatic permafrost, the biogeochemical parameters total nitrogen (TN), total carbon (TC) and total organic carbon content (TOC), and their stable isotopes (δ 13 C and δ 15 N) were measured and the TOC/TN ratio determined. The parameters TOC/TN and δ 13 C are widely used as a broad indicator for the degree of decomposition and source of organic matter (OM) (e.g. Mary et al., 1992;Meyers 1994;Meyers, 1997;Gundelwein et al., 2007;Andersson et al., 2012;Strauss et al., 2015). Generally, a lower TOC/TN ratio and higher δ 13 C indicate more decomposed organic material (Schirrmeister et al., 2011;Strauss et al., 2015). Since microorganisms prefer easily accessible carbon compounds, the organic carbon quality decreases with increasing organic matter decomposition (Schowalter, 2016).
Using δ 13 C and TOC/TN, it is possible to distinguish between land plants and marine, brackish and freshwater algae (decreasing δ 13 C is associated with decreasingly saline environment) (Meyers et al., 1994) (Supplementary Table S2).
To determine TC and TN content, the homogenized, freeze dried and ground samples were measured with a carbon-nitrogen-sulfur analyzer (Vario EL III, Elementar). TOC was measured with a TOC analyser (Vario Max C, Elementar). We measured two replicates of each sample and accepted only <5% deviation for the replicates. The TOC/TN ratio is expressed as atomic TOC/TN value to ensure comparability between studies. The atomic TOC/TN value was calculated following Meyers et al. (1994) by multiplying the weight ratio with 1.167, which is the ratio of the atomic weights of nitrogen (14.007 amu, atomic mass unit) and carbon (12.001 amu). The rate of mineralization was considered high at TOC/TN < 12, moderate between 12 and 25, and low at values >25 (Walhert et al., 2004).
The volumetric TOC content (carbon density, kg TOC/m 3 ) was calculated following Strauss et al. (2013) according to Eq. 1: where V ref is the reference volume of 1 m 3 , BD in kg/m 3 and TOC in wt%. The stable carbon and nitrogen isotopic composition was measured for all samples with TOC or TN values above the detection limit (0.1 wt%) with a ThermoFisher Scientific Delta-V-Advantage gas mass spectrometer equipped with an organic elemental analyser Carlo-Erba NC2500 (accuracy: < 0.2‰). Measured δ 13 C and δ 15 N was compared to the standardized Vienna Pee Dee Belemnite (VPDB) and expressed in per mille (‰ vs. VPDB). TOC and TN values below the detection limit were neglected in further calculations.

Geochronology
For radiocarbon dating, we chose samples about every 3 m (n 19 for both cores). When present, macroplant remains were handpicked under a stereomicroscope, otherwise bulk sediment was selected and dated using Accelerator Mass Spectrometry in the AWI MICADAS (MIni radioCArbon DAting System) Laboratory in Bremerhaven. Radiocarbon dates were calibrated using Calib 7.04 software after (Reimer et al., 2013) and the IntCal13 calibration curve (Stuiver et al., 2020) to calculate calibrated years before present (cal. years BP).
For optically stimulated luminescence (OSL) dating, three core segments were selected in the field from the Uomullyakh core, and two from the Polar Fox core. The core segments were frozen in their original state, wrapped in black foil, and delivered frozen to Aarhus University's Nordic Laboratory for Luminescence Dating (Roskilde, Denmark). The core segments were then thawed in a controlled-light environment and the outer ∼5 mm removed to minimize the effects of smearing, and to ensure that the material used for OSL measurements had not been exposed to daylight during or after retrieval. Luminescence measurements were made using a Risø TL/OSL reader, model TLDA 20. A standard SAR protocol based on Murray and Wintle (2000;2003) was used for quartz dose estimation, and a comparison with feldspar data was used to investigate the degree of bleaching at deposition (e.g. Möller and Murray, 2015). Radionuclide concentrations were measured using high resolution gamma spectrometry (Murray et al., 1987;Murray et al., 2018), calibrated using Certified Reference Materials produced by Natural Resources Canada (NRCAN). For detailed information see supplementary material.

Borehole Temperature Measurements and Calculation of Freezing Conditions
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.
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.

Lagoon Mapping
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, semiclosed 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 2 .

Upscaling Carbon Pools
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: thickness · coverage · BD · TOC 100 10 6 (3) with thickness and coverage in m, BD in kg/m 3 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.

Statistics
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 δ 13 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'.

RESULTS
The cores are described according to the depositional record from bottom to top. All depth specifications refer to the top of the ice cover (contains a compacted snow layer on top of lagoon ice) unless they are specifically stated as sediment depth, in which case we refer to the depth from the top of sediment. The units were separated according to a vertical cut in the dissimilarity dendrogram at height 7 for the Uomullyakh sediment core ( Figure 2) and at height 4 for the Polar Fox sediment core ( Figure 3). Therefore, the Uomullyakh core (PG2410, core recovery of 85%) was divided into three units, whereby Unit I was divided into two subunits: UoL-UIII ranges from 3350 to 1480 cm, UoL-UII from 1480 to 820 cm, UoL-UIa from 820 to 320 cm and UoL-UIb from 320 to 120 cm depth. The Polar Fox FIGURE 2 | Summary of sedimentological, geochronological and biogeochemical parameters for the Uomullyakh Lagoon core: bulk density (upper axis), mass specific magnetic susceptibility (lower axis), grain size (upper axis: volumetric proportions, lower axis: mean grain size), total organic carbon content (TOC), total nitrogen content (TN), organic carbon and nitrogen ratio (TOC/TN), stable carbon and nitrogen isotope composition (δ 13 C, δ 15 N), radiocarbon dating ( 14 C) and luminescence dating (pIRIR, OSL). The top 120 cm is the overlying ice cover. The core was divided into three stratigraphic units based on a CONISS analysis. Color bar at the y-axis indicates frozen (blue), unfrozen (red) and partially frozen (shaded) core segments. Note that "frozen" includes visibly frozen core sections and re-constructed frozen conditions. Arrows indicate non-finite ages.

Core Description
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).
FIGURE 3 | Summary of sedimentological, geochronological and biogeochemical parameters for Polar Fox Lagoon core: bulk density (upper axis), mass specific magnetic susceptibility (lower axis), grain size (upper axis: volumetric proportions, lower axis: mean grain size), total organic carbon content (TOC), total nitrogen content (TN), organic carbon and nitrogen ratio (TOC/TN), stable carbon and nitrogen isotope composition (δ 13 C, δ 15 N), radiocarbon dating ( 14 C) and luminescence dating (pIRIR). The top 330 cm is the overlaying ice and water layer. The core is divided into three stratigraphic units based on a CONISS analysis. Color bar at the y-axis indicates frozen (blue), unfrozen (red) and partially frozen (shaded) core segments. Arrows indicate non-finite ages.  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 nonfinite 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.

Uomullyakh Lagoon
The sediments of UoL-UIII ( Figure 2) were very poorly sorted. Up to 2000 cm depth, sand dominated with a share of up to 96 wt%, mainly fine to medium sand (measured mean grain size of the < 1 mm fraction: 161 µm) and visually observed rounded pebbles up to 2 cm in size. The upper part of UoL-UIII was also characterized by an alternation of coarser deposits (mean grain size: 193 µm) with low MS of 20 × 10 -8 m 3 /kg and fine deposits (mean grain size: 37 µm) with high MS values of 59 × 10 -8 m 3 /kg. At the lower end of UoL-UII, the TN content decreased below the detection limit of 0.1 wt%. In UoL-UII, coarse sand layers with pebbles and silty fine sand layers alternated. The sand fraction decreased (mean grain size: 24.3 µm). At 1110 cm, TOC (6.7 wt%) and TN (0.4 wt%) values were at a maximum and were associated with low δ 13 C (−27.5‰ VPDB) and low δ 15 N (1.3‰) values (Supplementary Table S3). In Unit Ia of UoL, poorly sorted clay and silt deposits were predominant with a share of 16 and 73% respectively (mean grain size: 12.5 µm) (Figure 2). The TOC content ranged between 3.2 and 4.5 wt% and the TN content between 0.2 and 0.3 wt%. UoL-UIb, containing the upper 200 cm, was distinguished from the lower unit UoL-UIa by coarser grain size and lower TOC and TN content (1.9 wt% and 0.15 wt%, respectively). In samples with a mean grain size below 63 µm (mainly Units I and II), the δ 13 C values were mostly below −26‰ VPDB, whereas they were higher in coarse-grained samples. The bulk density decreased from the core bottom to the top (1600-100 kg/m 3 , mean 1200 kg/m 3 ). For further details of the biogeochemical results see Supplementary  Table S6. The MS ranged from 5 to 154 × 10 -8 m 3 /kg, with variations mainly in areas of abrupt changes in grain size. Fine deposits corresponded with high TOC content over the core length.

Polar Fox Lagoon
The lowermost Unit III (PFL-UIII) stretching up to 2080 cm depth was characterized by a large proportion of coarse material (sand and gravel, measured mean grain size of the <1 mm fraction: 230 µm) and a very low organic carbon content of 0.1 wt%. The sediment of this unit was mostly moderately well sorted.
In the overlying PFL-UII the sorting was poor. Silty fine sand and medium sand layers alternated (mean grain size: 72.3 µm) along with MS changes. Changes in MS correspond to rapid increases in grain size (at 630 and 530 cm depth). MS ranged between 6 and 66 × 10 -8 m 3 /kg, with particularly strong variation in the small frozen core sections. Overall, the Polar Fox Lagoon core had coarser material composition of the < 1 mm fraction than the Uomullyakh Lagoon core. The TOC content was low (0.7 wt%) and the TN content decreased to below the detection limit at 1200 cm depth.
The uppermost Unit I (PFL-UI) comprised approximately the uppermost 5 m of sediment (corresponding to a depth of 860 cm) and was characterized by clay, silt and fine sand (mean grain size: 15.1 µm). Samples of the upper and the lower half of this unit show the same pattern of decreasing grain size, associated with increasing TOC and TN contents, as well as an increase in TOC/ TN ratio and a depletion in δ 13 C. In contrast to UoL-UI, TOC, maximal at 4 wt%, was highest in the uppermost sample (Supplementary Table S7). TN ranged between 0.14 and 0.27 wt%. Changes in MS corresponded to rapid increases Frontiers in Earth Science | www.frontiersin.org July 2021 | Volume 9 | Article 637899 in grain size (at 630 and 530 cm). MS ranged between 6 and 66 × 10 -8 m 3 /kg, with particularly strong variation in the small frozen core sections.

Hydrochemistry and Reconstructed Freezing Conditions
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  Figure S3A). In the UoL core, the pore water of the upper sample had the lightest isotope composition with δ 18 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 δ 18 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 FIGURE 4 | Hydrochemical parameters of the pore water and borehole temperature of the Uomullyakh Lagoon (left) and the Polar Fox Lagoon cores (right): pH, electrical conductivity (EC), dissolved organic carbon concentration (DOC), measured borehole temperature (point) and freezing point (line) estimated from pore water salinity. Color bar at the y-axis indicates frozen (blue), unfrozen (red) and partially frozen (shaded) core segments. Note that at UoL "frozen" includes visibly frozen core sections and re-constructed frozen conditions. Frontiers in Earth Science | www.frontiersin.org July 2021 | Volume 9 | Article 637899 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.

Carbon Pool Calculation
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 2 , respectively. The total thermokarst lagoon area is 11.4 km 2 .

PCA
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  Frontiers in Earth Science | www.frontiersin.org July 2021 | Volume 9 | Article 637899 with fine deposits and high organic matter content and quality, while in the deeper core, the variables sand content and δ 13 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.

Depositional History and Sediment Facies
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., 2017Schirrmeister et al., , 2018Jongejans 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., , 2011Wetterich 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 . 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 , 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 . 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-tobleach feldspar IR 50 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.

Middle/Late Pleistocene Fluvial Deposition -UoL-UIII (3350-1480 cm) and PFL-UIII + II (3070-1460 cm)
The first phase of landscape history at the study area is mainly characterized by fluvial deposits. Under the Polar Fox Lagoon, they range in age from >360 to >320 ka (pIRIR) and include the units PFL-UIII and the lower part of PFL-UII (in total 1590 cm). Below Uomullakh Lagoon, these deposits (Uol-UIII) have an age between >300 ka (pIRIR) to 63 ± 2 ka (OSL). This long time period is represented by only about 9 m of sediment, indicating an overall low accumulation rate under a dynamic fluvial regime with frequently alternating phases of accumulation and erosion in Unit III.
The presence of coarse sediments, rounded gravel and pebbles up to 2 cm, as well as driftwood indicates fluvial transportation and deposition (Supplementary Figure S4). Several studies also found fluvial sediments in similar depth underneath the neighboring Ivashkina Lagoon Schirrmeister et al., 2018) and underlying Yedoma deposits on the Bykovsky Peninsula (Slagoda 1993;Slagoda, 2004;Siegert et al., 2002;Grosse et al., 2007;. Kunitsky (1989) and Wetterich et al. (2008) assumed that the paleo-Lena River crossed the present Bykovsky Peninsula in an early Weichselian period (MIS 4 and maybe older). Our luminescence ages suggest that river transport existed over a long period from >360 ka (pIRIR) to 63 ± 2 ka (OSL).
While in UoL-UIII, discontinuous grain sizes, very poor sorting and higher organic matter content reflect frequently varying water runoff in a shallow river branch or near-shore area, PFL-UIII is distinguished by very little organic material, more continuous and coarser grain sizes, and a higher degree of sorting. Such properties reflect stable fluvial current conditions (Wetterich et al., 2008). This could suggest that UoL core was located in the near-shore area and PFL core closer to the center of a paleo river. However, the large distance between both cores Frontiers in Earth Science | www.frontiersin.org July 2021 | Volume 9 | Article 637899 (2.5 km) and the deposition of fine grained sediment underneath the Goltsovoye Lake, located in the middle of both lagoons, makes this scenario improbable. It is more likely that two neighboring river branches, part of a braided river system or delta run through the location of today's lagoons. This assumption is encouraged by finds of fluvial deposits across the peninsula. Alternation of fine and coarse grained layers at UoL-UIII indicate ongoing changes in river morphology. In PFL-UII until 1200 cm depth, organic containing fine sand layers were interbedded with coarser sand and gravel layers similar to the upper part of UoL-UIII also indicating varying water runoff or change in river course. A layer of large driftwood pieces of similar shape and size found in UoL-UIII (at 1530 cm), PFL-UII (at 1810 cm) as well as GoL (at 2045 cm) (Unit II in  are presumably deposited at the same flood event. In comparison to PFL and GoL, UoL core showed coarser grained fluvial deposits after the flood event which go along with increased erosion of overlying sediment explaining the higher elevation of the driftwood layer at UoL. Based on biomarker analyses carried out for the Goltsovoye Lake core,  concluded that the area at time of flooding was a wetland, dominated by low-centered polygons formed during the Kargin Interstadial (MIS 3).

Late Weichelian (MIS 2) Yedoma Deposition -UoL-UII (1480-820 cm) and PFL-UII (1200-860 cm)
In the second phase of the landscape development an increase in silt deposition and organic matter content, as well as a depletion of δ 13 C point to changing environmental conditions away from high energy fluvial processes toward lower energy alluvial and other transport processes. Luminescence dating for UoL at 1107 cm to 37 ± 2 ka (OSL) revealed a late Pleistocene age (beginning of late Weichelian stadial). Radiocarbon dates from sediments of similar depth for both UoL (43 ± 0.56 cal ka BP at 1063 cm) and PFL (36.1 ± 0.39 cal ka BP at 1128 cm) are in a similar time range. Paleo-ecological data from the nearby Mamontovy Khayata Yedoma cliff (Andreev et al., 2002;Schirrmeister et al., 2002) and biomarker analyses from GoL  revealed that the climate was drier and colder, transforming the former wetland into a steppe-like tundra, which probably caused decreasing river run off and changing river morphology.
When comparing the data with those of the Yedoma Ice Complex of Mamontovy Khayata, they appear to resemble mostly the fine-grained layer of the units Uo-UII and PFL-UII. Thus, it can be deduced that ice-rich Yedoma Ice Complex deposits started accumulating when either the previously existing river fell dry or changed its course.
Holocene Thermokarst Lake Formation -UoL-UIa (820-320 cm) and PFL-UI (860-330 cm) The third phase involves the formation of a thermokarst lake, which is a prerequisite for a thermokarst lagoon. High proportions of clay and silt as well as high TOC and TN content and shell remains (PFL at 400-440 cm) indicate that Unit I of both cores represent the lake stage. Based on the results of the CONISS analysis, we discovered ∼5.5 m thick lake sediments for both lagoons, which are likely mixed with Yedoma deposits from the lake shore.
In northeast Siberia, thermokarst lake formation began with rapid warming and wetting between 14 and 11.7 ka (Kaplina and Lozhkin, 1979;Kaplina, 2009;Anthony et al., 2014;Subetto et al., 2017). Grosse et al. (2007) estimated that major lakes and basins on the Bykovsky Peninsula were formed in the early Holocene, probably starting between 12.5 and 9.4 cal ka BP when a sedimentation gap on top of local Yedoma deposits suggests widespread reorganization of the landscape. This is in agreement with a radiocarbon date of 10.5 ka cal BP at 480 cm for UoL-UIa which we interpret as lacustrine deposits.
The gravimetric ice content of Yedoma deposited at Mamontovy Khayata is between 60 and 180% . Assuming that the Yedoma deposits were similarly ice-rich underneath the lagoons, thawing of these ice-rich sediments explains the large surface subsidence discussed before and suggests a substantial deposit compaction during lake stage.
The size of the drained thermokarst lake basins, shown in the landform classification map of Fuchs et al. (2018) (Supplementary Figure S5), indicates that Uomullyakh and Polar Fox basins were once significantly larger than the current water bodies and that drainage events occurred in the past. Lake drainages are common for thermokarst landscapes and are often rapid events with formation of drainage channels cutting deep into ice-rich permafrost .
In UoL-UIa, a rapid increase in grain size at 630 and 530 cm depth corresponds to erratic changes in mass specific magnetic susceptibility indicating an external input of sediment at least twice. Upstream, south-west and north-west of UoL, two drained lake basins (Supplementary Figure S5) are located. Drainage channels connecting the basins with UoL, still visible at satellite images, suggest that UoL has been flooded while these lakes drained. As lake drainage may cause sediment erosion and downstream deposition, it is possibly the reason for the external sediment and organic matter input found in the UoL core.
In PFL-UI, two repetitive patterns of decreasing grain size associated with increasing TOC and TN content and organic matter quality are visible and might indicate drainage events. Therefore, Polar Fox Lagoon may have drained twice during the lake phase. The first drainage could have been northwards toward the adjoining lake basin (Supplementary Figure S5) and the second drainage southwards to the sea leading to lagoon formation.
A wider TOC/TN ratio and lighter δ 13 C composition of the sediment beneath UoL indicates that the organic matter was less mineralized than beneath PFL. That is likely because the UoL-UIb horizon is frozen in winter due to bedfast ice formation and subsequent lagoon bottom freezing, while PFL-UI is mostly thawed throughout the year, allowing more time for microbial degradation of organic matter. Seasonal freezing under floating ice conditions at PFL is also possible but would likely only result in partial sediment freezing and a thinner seasonally frozen layer compared to UoL (Angelopoulos et al., 2020a).
Holocene Lagoon Formation -UoL-UIb (320-120 cm) Indications for the fourth phase of landscape development, the formation of a thermokarst lagoon, were found at the boundary between UoL-UIb and UoL-UIa (between 360 and 280 cm). The increase in mass specific magnetic susceptibility and grain size, the decrease in TOC and TN content, as well as a slight decrease in TOC/TN ratio indicate a system change. In our interpretation, this shift is caused by marine inundation and marks the beginning of the third stage, the lagoon formation. Postglacial rising sea levels caused by melting of continental ice sheets led to the inundation of former terrestrial permafrost on the Laptev Sea shelf until about 5 ka (Bauch et al., 2001) when sea levels stabilized and coastal erosion of ice-rich permafrost started to dominate the coastline (Romanovskii et al., 2000;Romanovskii et al., 2004). Likely, the Uomullyakh thermokarst lake was breached by the transgressing sea and coastal erosion in the early to mid-Holocene. A similar breaching and submergence of thermokarst lakes and drained lake basins has been observed in thermokarst-affected coastal lowlands on the Alaska North Slope (Arp et al., 2010), in the Mackenzie Delta region (e.g. Ruz et al., 1992), and on Banks Island (Grasby et al., 2013).
Satellite images (e.g. Figure 1) indicate that the Uomullyakh Lagoon is characterized by stronger turbidity due to lack of topographical wind shadow (compared to the Yedoma upland surrounded PFL) and its connection to the sea, which enables waves to enter. High connectivity to the sea likely leads to an increased input of marine sediment (Grotheer et al., 2020). This is in agreement with sediment analyses indicating marine deposition in Uomullyakh Lagoon for the upper 2 m mixed with lacustrine deposits.
At Polar Fox Lagoon, this distinction could not be made. This is very likely caused by the different morphology of both lagoons. The Polar Fox Lagoon was formed as a result of thermokarst lake drainage, which led to the formation of a shallow and winding channel connecting lake and sea (Angelopoulos et al., 2020a). The connection to the Tiksi Bay is cut off during winter when ice formation reaches the channel bed, restricting water and sediment exchange (Spangenberg et al., 2021).

Impact of Marine Inundation
Our data suggests that deep (∼20 m) hypersaline taliks (UoL) preserve organic carbon better than taliks that contain mostly Frontiers in Earth Science | www.frontiersin.org July 2021 | Volume 9 | Article 637899 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 icebearing 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 Frontiers in Earth Science | www.frontiersin.org July 2021 | Volume 9 | Article 637899 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 2 and CH 4 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 , dissolution in the pore water and surface water, lowpermeability 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.

Carbon Pool Upscaling
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 3 (Uo) vs. 26.5 kg/m 3 (PFL); upper 30 m: 17.4 kg/m 3 (UoL) vs. 10.2 kg/m 3 (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 3 , which is in the same range as the carbon density we calculated for the upper 3 m (24.5 ± 1.4 kg/m 3 ), but 45% higher than the mean carbon density of the 30 m-long sediment cores (13.8 ± 3 kg/m 3 ). 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  (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 (  (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 3 and therefore in the upper range of the C-density estimated for the 0.41 × 10 6 km 2 large Yedoma region (19 + 13/-11 kg/m 3 ) (Strauss et al., 2013) and in the lower range of terrestrial thermokarst deposits of the Yedoma region (33 + 25/-19 kg/m 3 ) (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 2 . This is a substantial inventory of formerly frozen C, which is now largely unfrozen and thus available for microbial degradation processes. δ 13 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 δ 13 C signals, but they differ from the over-or underlying sediment facies (Schirrmeister et al., 2018;. 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 . 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 .
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 2 relative to that of methane (Segers and Kengen, 1998). In line with this, substantial CO 2 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.

CONCLUSION
This study shows that thermokarst lagoons are highly dynamic landforms at the boundary between terrestrial permafrost and marine systems. Based on two ∼30 m long sediment cores from two contrasting thermokarst lagoon systems, we examined the depositional history, the impact of marine inundation and the amount of stored organic carbon for different stages of thermokarst lagoons. We identified several different phases of sediment deposition in the semi-closed Uomullyakh Lagoon and the nearly closed Polar Fox Lagoon. The sedimentological analysis of both lagoon cores shows that fluvial deposition from a paleoriver system dominated at the sites from >360 ka to the late Pleistocene, when silty to sandy ice-rich Yedoma deposits started to accumulate in this region. These former Yedoma deposits are now present as thawed talik sediments. With the beginning of the warmer Late Glacial to Holocene period, the formation of thermokarst lakes began and resulted in lacustrine sediment deposition at both lagoons. At the more open Uomullyakh Lagoon, marine sediments overlying lacustrine deposits indicate the onset of the lagoon phase starting in the early to mid-Holocene. Inundation with seawater and formation of seasonal bedfast ice caused the saline enrichment of the pore water, which led to cryotic talik formation.
Our case study from the Bykovsky Peninsula provides a first estimate of organic carbon quantity and quality based on deep sediment cores below thermokarst lagoons. We have identified large organic carbon stocks in these special Arctic lagoon systems. We measured a mean C-density of 15.3 ± 2 kg/m 3 for the 30 m thick deposits below the lagoons, with higher values found below the semi-closed lagoon. The measured carbon density for the Yedoma deposits (24.72 ± 4 kg/m 3 ) is in the upper range of the estimate of terrestrial Yedoma (19 + 13/-11 kg/m 3 ), but lower range of terrestrial thermokarst deposits (33 + 25/-19 kg/m 3 ) in the Yedoma region. Based on the mean C-density and the remote sensing-based spatial extent of thermokarst lagoons, we calculate a total C storage of 5.7 Mt C in the upper 30 m of sediment in five thermokarst lagoons on the southern Bykovsky Peninsula.
Our data indicate that open thermokarst lagoons under high seawater influence seem to better preserve organic carbon more then freshwater dominated systems. This shows the potential importance of the role and contribution of thermokarst lagoons in the transformation of terrestrial permafrost carbon under increasingly marine influence.

DATA AVAILABILITY STATEMENT
The datasets presented in this study can be found in online repositories. The names of the repository/repositories and accession number(s) can be found below: https://doi.pangaea. de/10.1594/PANGAEA.922169.

AUTHOR CONTRIBUTIONS
MJ and JS designed this study. MJ led the writing of the first draft of the manuscript. JS, GG, and MG developed the overall coring plans for the Bykovsky Peninsula field campaign. JS, MG, LS, BB, PO, MA, SL, and GG conducted the field work. MA and PO did the subsampling for the Polar Fox sediment core and MJ for the Uomullyakh core. MJ and PO performed laboratory analyses. MJ led the data interpretation for both cores. MJ and IN conducted the spatial analysis and mapping of the lagoons. AM performed optical luminescence analyses and data interpretation. All co-authors contributed within their specific expertise to data interpretation.

FUNDING
This study was carried out within the NERC-BMBF project CACOON (Changing Arctic Carbon cycle in the cOastal Ocean Near-shore, grant no. 03F0806A), the ERC Project PETA-CARB (#338335), and the BMBF project KoPF (03F0764B, 03F0764F). The Alfred Wegener Institute Helmholtz Centre for Polar and Marine Research (AWI) provided baseline funding and together with the Germany Research Centre for Geosciences (GFZ) AWI provided funding for the expedition.