The deposits below both lakes were drilled from the lake ice in March 2015 during a joint German-Russian field expedition. The drilling depths reached down to about 20–21 m below the lake-ice surface. We used a URB2-47 drilling rig without a core catching system. Due to partly unfrozen sediment layers, parts of the cores were lost during the drilling process. The retrieved core segments were pushed out with air pressure, described and stored in plastic bags and kept frozen until laboratory analysis. From here on, all core depths are given in centimeters as mean depth of a sampling interval below lake-ice surface (cm bls).

The frozen cores were split lengthwise using a band saw and were subsequently subsampled in a freezer room at −5°C. After cleaning, the cores were described visually and photographed (Figure 3). The sampling half was cut into subsamples and divided for further analysis. We did a regular sampling on both cores according to visual stratigraphy. Hence, the alas lake core YU-L7 was sampled in intervals between 20 and 70 cm (n = 19) and the Yedoma lake core YU-L15 at intervals between 10 and 40 cm (n = 49).

For the biogeochemical and sedimentological analyses, the gravimetric water and/or ice content was determined as the difference between wet mass of a sample (as ice for in situ frozen parts) and the dry sample mass and is expressed in weight percent (wt%). Unless stated otherwise, we will use the term gravimetric water/ice content (GW/IC) onwards, due to changing unfrozen and frozen conditions within the cores.

We radiocarbon-dated seven samples from YU-L7 and eight from YU-L15 using Accelerator Mass Spectrometry (AMS) with the Mini Carbon Dating System (MICADAS) at Alfred Wegener Institute Helmholtz Centre for Polar and Marine Research (AWI) Bremerhaven, Germany (Mollenhauer et al., 2021). Macrofossils such as plant and wooden remains (i.e. roots, leaf fragments, twigs, Larix spp. (larch) needles) were picked and dated for nine samples. For the other six samples, bulk sediment was dated. The radiocarbon ages were calibrated using the OxCal4.4 software (Bronk Ramsey, 2009) with the IntCal20 calibration curve after Reimer et al. (2020), and are expressed in calibrated years before present (cal. yrs. BP).

All subsamples of the YU-L7 alas lake core (n = 19) were analyzed for sedimentology. Of the Yedoma lake core YU-L15, all samples above 1,362 cm and below 2064 cm bls were analyzed but only every second subsample from the middle part (n = 31).

A Bartington Magnetic Susceptibility Meter Model MS2 (Bartington Instruments, United Kingdom) was used to measure the mass-specific magnetic susceptibility (MS). The MS reflects differences in magnetizability of clastic material. Measurements were done at low frequency of 0.465 kHz in an external magnetic field in 10⁻⁸ m³/kg.

Grain-size analyses were carried out on organic-free (treated with 35% H₂O₂) sub-samples using a Malvern Mastersizer 3000 (Malvern Instruments, United Kingdom) with an attached Malvern Hydro LV wet-sample dispersion unit. The proportions of sand, silt and clay fractions are given as sums between 2 mm and 63 µm, 63 and 2 µm and <2 µm, respectively. Grain-size parameters were calculated with the software Gradistat (Version 8.0; Blott and Pye, 2001).

We analyzed the total carbon (TC) content using an Elementar Vario EL III (Elementar Analysensysteme GmbH) (detection limit of 0.05 wt%). Furthermore, total organic carbon (TOC) measurements were performed using the Elementar Vario Max C (detection limit of 0.1 wt%). No external removal of the total inorganic carbon (TIC) is needed for this method. The TIC content was calculated as the difference of TOC and TC. All values of TC, TOC, TIC are calculated in wt% (n = 68).

The stable carbon isotope ratio (δ¹³C) was determined using a ThermoFisher Scientific Delta-V-Advantage gas mass spectrometer equipped with a FLASH elemental analyzer EA 2000 and a CONFLO IV gas mixing system in the Stable Isotope Laboratory at the AWI Potsdam, Germany. An aliquot per sample of 15 samples per core with TOC values above detection limit were measured. The δ¹³C measurements have been calibrated by three international reference materials. The results are compared to the common Vienna Pee Dee Belemnite (VPDB) standard and expressed in per mill (‰ vs. VPDB; n = 30; measurement accuracy: equal or smaller than 0.15‰). The δ¹³C values give information about plant metabolism, as they are generally very similar to those of the former vegetation. We interpret more negative δ¹³C values as an indication of less-degraded OM, while less negative δ¹³C values are linked to stronger decomposition processes (Schirrmeister et al., 2011; Strauss et al., 2015).

Pore water was extracted from sediment in order to measure stable water isotope ratios. The samples were thawed at 4°C and pore water was extracted using Rhizone soil moisture sampler (Eijkelkamp Soil & Water). The pore size of the membrane was between 0.12 and 0.18 µm. Several sediment subsamples did not hold enough water for pore water extraction. Thus, stable water isotope measurements could be carried out for 42 subsamples of YU-L15 and five subsamples for YU-L7.

Stable water isotopes, i.e. the ratios of hydrogen (δD) and oxygen (δ¹⁸O), were measured using a Finnigan MAT Delta-S mass spectrometer following standard procedures described in Meyer et al. (2000). The isotopic ratios of δD and δ¹⁸O were calculated using the software ISODAT and expressed in per mill (‰) relative to the Vienna Standard Mean Ocean Water (VSMOW) value. The analytical precision is better than ±0.1‰ for δ¹⁸O and ±0.1‰ for δD (Meyer et al., 2000). We furthermore calculated the second-order parameter deuterium excess d (d = δD-8*δ¹⁸O, Dansgaard, 1964). While δ¹⁸O and δD in precipitation are considered as temperature proxies (lower values reflect lower temperatures, higher values reflect higher temperatures), d provides information on moisture generation processes, mainly relative humidity and sea surface temperature. Freeze-thaw processes in the active layer cause secondary fractionation processes, which may overprint the original isotopic signatures.

All samples were analyzed for their element concentrations by non-destructive X-ray fluorescence (XRF) spectrometry. Therefore, bulk samples were homogenized, mixed with a wax binder (CEREOX Licowax at a ratio of 4 to 1), and pressed into 32 mm pellets. Measurements were carried out using the Energy-Dispersive Polarisation XRF (EDPXRF) SPECTRO XEPOS (SPECTRO Analytical Instruments Ltd.) analyzer in a helium gas atmosphere. The contents of all elements from sodium to uranium were simultaneously determined and adjusted to sample weight. After intensive literature research, analyzing the measured single elements, and calculating a large number of ratios, we concentrated on six elemental ratios (S/Cl, Ca/Ti, Sr/Ca, Fe₂O₃/MnO, Si/Al, Zr/Al) for detailed discussion of sediment properties and permafrost conditions during Yedoma deposition as well as thermokarst and talik processes. We applied the oxide representation of iron and manganese because both, iron and manganese are characterized by several oxidation states compared to the other element ratios. Thus, the ratio in the oxidized form would be more understandable and easier to compare with literature data. The nominal values of the ratio differ only slightly but would lead to the same conclusions. All measured single element contents for both cores are presented in the Supplementary Tables 1, 2.

Multi-element weathering indices provide, relative to the parent material, semi-quantitative measures of mineral weathering intensity, which is directly dependent on, e.g., water availability and temperature conditions. They are traditionally calculated from the concentrations of several mobile (e.g., Ca, K, Na, Mg) and immobile elements (e.g., Al, Fe, Ti, Zr) (Schatz et al., 2015). From the large number of weathering indices (Schatz et al., 2015, cf. Table 1), we ultimately selected three indices that best correspond to geochemical conditions in frozen and thawed Yedoma deposits: the FENG, MIA, and ICV indices.

Feng (1997) developed the FENG index in particular for calcareous substrates. In contrast to many other indices, Ca is not taken into account in this index. Since the solution and relocation history of Ca is almost exclusively traced in the case of clearly calcareous initial substrates, misinterpretations can result (Buggle et al., 2011). By placing easily soluble elements in the denominator of the formula, an increase in the FENG index indicates increased weathering. In accordance to Feng (1997), the index is calculated in molar proportions using Eq. 1. FENG = ( Al 2 O 3 + Fe 2 O 3 ) / Na 2 O + K 2 O + MgO + P 2 O 5 (1)

Babechuk et al. (2014) proposed the mafic index of weathering (MIA) especially for Mg and Fe-containing minerals. The specific arrangement of the MIA_(R) was used in consideration of predominantly reducing soil conditions. This takes into account the redox-dependent weathering behavior of iron. Under reducing conditions, Fe²⁺ responds as a mobile element and is leached along with Mg during the mafic weathering process (Babechuk et al., 2014). Hence, the MIA_(R) index increases with increased weathering. In molar proportion and considering only the silicate-bound Ca (CaO^*), the MIA_(R) calculation is: MIA ( R ) = 100 x [ Al 2 O 3 / ( Al 2 O 3 + Fe 2 O 3 ( T ) + MgO + CaO * + Na 2 O + K 2 O ) ] (2)

The index of composition variability (ICV) was defined by Cox et al. (1995) and can be applied as a measure of compositional maturity. The index measures the abundance of alumina in relation to other major cations (in molar proportions; Cox et al., 1995). ICV = ( Fe 2 O 3 + K 2 O + Na 2 O + CaO + MgO + MnO + TiO 2 ) / Al 2 O 3 (3)

An index value close to 1 shows hardly any mineral conversions, while an index value less than 1 indicates several cycles of weathering, erosion and accumulation; i.e. the smaller the index, the stronger the weathering. This index is therefore inversely proportional to the FENG and the MIA_(R) index.

For the examination of relationships between sedimentology as well as organic and inorganic geochemistry of frozen and thawed Yedoma deposits, statistical analyses were performed. Descriptive statistics were carried out using the Python software packages Pandas (McKinney, 2010) and NumPy (Oliphant, 2006). For calculation of Pearson correlation matrices including p-values and to compare the two sediment cores by Mann–Whitney U nonparametric tests the Pingouin statistical package in Python 3 was used (Vallat, 2018). Principal component analysis (PCA) was performed on all pre-selected sedimentological, organic, and inorganic geochemical parameters. Ca, S and Cl were included, additionally, in order to evaluate the influence of these elements on individual elemental ratios. The PCA was carried out using the Python package Scikit-learn (Pedregosa et al., 2011). Since missing values are incompatible with Scikit-learn estimators, the missing values of the grain size parameters from the Yedoma lake core YU-L15 (see above) were initially replaced using the k-Nearest Neighbors approach. By default, an Euclidean distance metric that supports missing values was used to find the nearest neighbors. For more information on methodology see Troyanskaya et al. (2001). For running the PCA, all variables were standardized by subtracting the mean and scaling the actual data range to unit variance (i.e., z-transform). The resulting ordination diagram presents standardized metric scores (i.e. scale_i = i × 1.0/(i.max – i.min)) and expresses the relationship among metrics as correlations. Additionally, all core samples were projected onto the ordination graph for interpretation purposes only, because they do not affect calculations. Finally, the sample scores on the first and second principal component (PC1 and PC2, respectively) were plotted against core depth.

Dating results are shown in Table 1. The mean calibrated ages of the YU-L7 alas core showed a general trend over depth ranging from 3,740 cal. yrs. BP at 520 cm bls, over 13,750 cal. yrs. BP at 758 cm bls to 43,020 cal. yrs. BP at 1,998 cm bls. The dating of plant remains at 1,560 cm bls showed an inverse age (49,070 ± 970 cal. yrs. BP). While the youngest calibrated age from the uppermost core part of YU-L7 came from plant remains dating, the other two dates resulted from bulk sediment dating. Three infinite ages at 1,278 cm bls (>24,570 yrs. BP), 1998 cm bls (>39,600 yrs. BP), and 1746 cm bls (>49,300 yrs. BP) generally fitted into the late Pleistocene age sequence.

The mean calibrated radiocarbon ages of the Yedoma lake core YU-L15 showed a relatively short age sequence during the interstadial Marine Isotope Stage (MIS) 3 with a bulk sediment ages of 41,870 cal. yrs. BP at 783 cm bls and 43,540 cal. yrs. BP at 1,917 cm bls. The dating of the plant remains at 1,917 cm bls (43,180 cal. yrs. BP) suggested an age similar to that of the bulk sediment from the same depth, thus confirming the period in which sediments were deposited in this core depth. Plant remains from the uppermost core sample were dated to a nearly modern age (140 ± 80 cal. yrs. BP). The inverse mean age at 1,395 cm bls (51,420 cal. yrs. BP) showed an age near the limit of the ¹⁴C calibration curves (Reimer et al., 2020). Three infinite radiocarbon dates of the YU-L15 core also point to late Pleistocene age of the core (Table 1).

Stable water isotopes were plotted in relation to the Global Meteoric Water Line (GMWL) (Craig, 1961) for each core (Figure 8A) and we used the season-specific Local Meteoric Water Line (LMWL) for Yakutsk from Papina et al. (2017). While the LMWL for the warm season exhibited a lower slope of 7.2 and a negative intercept of −18.9 in comparison to the GMWL, the LMWL for the cold season ran parallel but above the GMWL (Figure 8A). The Local Evaporation Line of Central Yakutia (LEL; Wetterich et al., 2008) reflects regional evaporation effects on the surface waters in the study area.

The δ¹⁸O and δD values of the alas lake core YU-L7 strongly differed from the GMWL and LMWLs, respectively, but their slope (5.7) roughly corresponded to the LEL (slope of 5). The δ¹⁸O values ranged from −13.7‰ at 1,302 cm bls to −12.9‰ at 1796 cm bls. The δD values ranged from −134‰ (1,649 cm bls) to −129‰ (1,328 cm bls). The d values varied between −27.6 and −24.9‰ and slightly decreased towards the bottom of the core.

The Yedoma lake core YU-L15 showed a much steeper slope (11.5) than GMWL, LMWLs, and LEL, respectively, with δ¹⁸O values between −29.5‰ at 1,316.5 cm bls and −25.3‰ at 1,931.5 cm bls (Figure 8B). The δD values range between −230‰ (1,316.5 cm bls) and −184‰ (1931.5 cm bls). Notably, the isotopic values from the uppermost frozen deposits below the talik (1,265–1,483 cm bls) were highest, while the samples between 1820 and 1992 cm bls showed the lowest and most depleted isotopic values. Generally, the δ¹⁸O and δD values decreased with depth until 1931.5 cm bls to −29.5‰ and to 230.4‰, but then increased again to the bottom layer to values that were more similar to the uppermost layers. The d values of the Yedoma lake core showed a similar course and lie between 5.1‰ (1957.5 cm bls) and 22.6‰ (1,377.5 cm bls).

Due to their inverse proportionality, a positive correlation (r) of the FENG or MIA_(R) index with a variable corresponds to a negative correlation with the ICV. High negative (respectively positive for ICV) and statistically significant (p < 0.05) correlations are obvious for both cores between the weathering indices and TOC, TIC, TC, and the S/Cl ratio (see Supplementary Tables 3, 4). Further positive correlations for both sediment sequences exist between TOC and the clay and silt grain size fraction. For the YU-L15 Yedoma lake core, there are particularly high positive (respectively negative for ICV) and statistically high significant (p < 0.001) correlations between the three weathering indices and the sand content, the mean grain size, the Si/Al ratio, and the Sr/Ca ratio. A negative and high significant correlation between GW/IC and the sand content can only be seen in YU-L15 (r = −0.64, p < 0.001). In addition, also all three weathering indices show significant negative (resp. positive) correlation with GW/IC in YU-L15. In contrast to YU-L15, all weathering indices show significant (p < 0.05) positive (respectively negative) correlations with depth in the YU-L7 alas lake core.

The PCA of both cores show two major principal components (PC1 and PC2) that explain 68.8% (YU-L7) and 69.2% (YU-L15) of the total variance of all sedimentological and biogeochemical parameters. The sample scores on PC1 and PC2 are plotted against depth for both cores in Figures 5, 7, illustrating their variations against the ordination results within the different sediment units. However, the ordination biplots (Figure 9) show different combinations of variables on the two axes. While the variables within the PCA for the Yedoma lake core YU-L15 are clustered quite clearly, for the alas lake core YU-L7 they show a broad distribution in small groups across the entire biplot. The vector lengths (i.e. PCA loadings) within the biplots indicate how strongly the individual variable is related to the displayed ordination. Hence, for both sediment cores, the highest positive loadings on PC1 are shown for Si/Al, Sr/Ca, Sand, the FENG index, and the MS. The MIA_(R) index also shows a positive correlation with PC1 for YU-L15. High negative loadings on the first axis for both cores are shown for TC, clay, silt, sorting, TIC, and TOC. In addition, negative correlations with PC1 are also shown by the variables ICV and GW/IC for YU-L15 and the variables S/Cl, S, Ca and Cl for YU-L7. Variables with high positive loadings on PC2 for YU-L7 are Ca/Ti, ICV, and GW/IC and for YU-L15 S and S/Cl. Negative correlations with PC2 for YU-L7 are shown for MIA_(R), Fe₂O₃/MnO, and depth, while for YU-L15, the variables depth, Zr/Al, Ca, Ca/Ti, and Cl are oriented along the negative range of PC2.

For our studied cores we found varying amounts of pore water/ice which is linked to the different host sediments and their respective water storage capacity (Figure 8). With a mean GW/IC value of 26.5 ± 13.8 wt%, the Yedoma lake core contains slightly more water/ice in comparison to the alas lake core (mean: 23.2 ± 9.5 wt%). The deposits of the alas lake core were completely thawed already. Thus, no further ground subsidence is expected for the ground below the alas lake. The higher GW/IC in the Yedoma lake core results from the perennial frozen conditions, particularly of the lower parts of the YU-L15 core. This together with further vertical and lateral expansion of the lake and the talik, would lead to more thawing of the sediments and further ground subsidence below the Yedoma lake (Fedorov et al., 2014).

In general, the composition of stable water isotopes in thermokarst lakes is mostly influenced by precipitation, evaporation and occasionally by meltwater from the surrounding catchment and the underlying permafrost (e.g., Turner et al., 2014). Earlier measurements of stable water isotopes of lake water from different thermokarst lakes at the Yukechi study site show evaporation by enriched δ¹⁸O (−13.8 to −8.8‰) and δD (−133.8 to −110.9‰) values, which is also seen in the LEL for Central Yakutia of δD = 5*δ¹⁸O - 61.4 (sample nos. Yak-01 to Yak-05; Wetterich et al., 2008).

Because of ground thawing (i.e., talik development) and lake water infiltration, the lake water has a strong influence on the pore water isotope composition of the talik sediments. In comparison to Northern Yakutia, the influence of local air temperatures and resulting evaporation effects on lake water δ¹⁸O are stronger due to much higher summer temperatures in Central Yakutia (Wetterich et al., 2008).

As shown in Figure 8, the δ¹⁸O and δD values from the alas lake core YU-L7 (mean: −13.2 ± 0.4‰ and −131.7 ± 2.2‰, respectively) are distinctly less depleted than those from the Yedoma lake core YU-L15 (−27.4 ± 1.3‰ and −205.0 ± 15.4‰, respectively). Furthermore, the alas lake core data are clearly independent from the GMWL but align with the LEL. Overall, they are similar but slightly more depleted than stable isotope values of the nearby drilled alas sediment core presented in Windirsch et al. (2020). Despite the low number of samples (n = 5), both the slope of 5.73 (R ² = 0.97) and the enriched isotopic composition with extremely negative d values suggest that the pore waters of the alas core originated from water infiltration from the alas lake above the thawed sediments. The alas lake water itself can be interpreted to represent a rather recent water source signal (i.e. precipitation), which is strongly affected by evaporative fractionation and further depends on cyclical water balance changes typical for the alas lakes in Central Yakutia (Ulrich et al., 2017b). The relatively homogenous isotope signals of YU-L7 result from the saturation of the sediments with surface water. This suggests subsurface water flow after talik development, which probably has also implications for local hydrology in a way that new flow paths through former impermeable layers are possible under warmer climate conditions (Johansson et al., 2015).

The Yedoma lake core (YU-L15) stable isotope data (n = 43) are highly correlated (R ² = 0.97) with a δ¹⁸O-δD slope of 11.53. The low δ¹⁸O and δD and high d values indicate that the water/ice originates from intra-sedimental pore and segregation ice of the still frozen, i.e. intact, Yedoma IC below the lake. Hence, they represent the original isotopic composition during permafrost formation. An overprinting by isotopically different lake water can be ruled out. Furthermore, the high slope and the high d values indicate that neither evaporation nor freezing fractionation have played a major role before fixation of the isotope signals in permafrost. While the upper part (above 1,800 cm bls) of the Yedoma lake core is dominated by structureless non-visible pore ice, the sediments below 1,800 cm bls represent a larger variety of intra-sedimental cryostructures: from structureless non-visible ice to varying amounts of micro and macro ice lenses to almost pure ice lenses (Figure 10). Hence, distinct minima in δ¹⁸O, δD, and d between 1,800 and 2,000 cm bls might be related to the higher amount of segregation ice even though the isotope profile below and above 1,800 cm bls do not differ considerably in terms of the δ¹⁸O-δD slopes. Furthermore, changes in climate and hydrological conditions, their seasonal variations as well as active-layer dynamics may have contributed to this pattern.

Overall, the mean δ¹⁸O and δD values of the Yedoma lake core are very similar to those of the nearby drilled on-land Yedoma core presented by Windirsch et al. (2020), whereas the mean d value of our core is lower, mainly driven by the distinct minimum between 1,800 and 2,000 cm bls. Nevertheless, our δ¹⁸O-δD slope (11.53) differs distinctly from that of the Windirsch et al. (2020) core (6.61) and is much higher than season-specific LMWLs for Yakutsk (Figure 8). Generally, the low δ¹⁸O and δD and the high d values are close to those of a nearby ice wedge (Windirsch et al., 2020) as well as other Central and Interior Yakutian MIS 3 ice wedges (Popp et al., 2006; Opel et al., 2019). This suggests a substantial contribution of winter precipitation (characterized by lower δ and higher d values) to the formation of the intra-sedimental ice. Comparable high pore ice d values have also been reported for some stratigraphic units exposed at the Batagay megaslump (Opel et al., 2019) and might therefore represent an over-regional pattern for Central and Interior Yakutia related to the extreme continentality of climate during MIS 3.

Quantifying the OM content is important to estimate the relevance of the study site regarding potential future GHG release from Yedoma deposits (Turetsky et al., 2020). The TOC for the YU-L7 core has a mean value of 0.72 wt%, whereas the mean value for YU-L15 core lies at only 0.37 wt%. These values are in the marginal range of averaged TOC values reported by Strauss et al. (2013) for Siberian Yedoma deposits (3.0 + 1.6/−2.2 wt%), but fit well with results of previous studies carried in Central Yakutia (Ulrich et al., 2019; Windirsch et al., 2020). In both cores, TOC correlates negatively with grain size (YU-L7: r = −0.48, p < 0.05; YU-L15: r = −0.43, p < 0.01; Supplementary Tables 1, 2). TOC, TIC, and TC values show local peaks in the fine-grained sediment units; most prominent at the silt shift at 1,672 cm bls (Y-SU-III; Figure 7). Extreme low TOC values are rather linked to the sand-dominated units.

One explanation for the low TOC content can be a strong post-sedimentary OC decomposition during thermokarst processes (Strauss et al., 2015; Weiss et al., 2016). But this could only apply for the thawed YU-L7 core sediments. The cryostratigraphy and stable water isotope values for the middle and lower parts of core YU-L15 indicate that at least this part stayed frozen since accumulation. Therefore, a deep OM degradation is unlikely here. However, short-term and rather shallow OM degradation within the active layer after different sedimentation phases of YU-L15 could be an explanation. A subsequent slowly rising permafrost table would leave enough time for the carbon mobilization before the sediment would freeze.

Another explanation for the low TOC values is a low initial OC content. Windirsch et al. (2020) argued that the fluvial and alluvial deposition of coarse-grained (sandy), organic-poor sediments led to the low TOC content. In addition, plant growth and organic input could have been hampered by substrate quality and disturbance (i.e. accumulation) frequency.

Further OM decomposition processes within the sediments below both thermokarst lakes are possible and traceable by δ13C isotopes (Schirrmeister et al., 2011; Strauss et al., 2015). For the talik unit Y-SU-V of the Yedoma lake core, δ13C indicate lower quality and higher carbon decomposition rates (Strauss et al., 2015; Weiss et al., 2016). The high quality signal in the upper part of the YU-L7 core (−26.12‰) likely resulted from modern input by the alas lake above. However, Jongejans et al. (2021a) found higher C/N ratio in the Yedoma lake core and in the bottom of the alas lake core compared to the top and middle parts of YU-L7, which suggest that OM was at least better preserved in frozen sediments below the Yedoma lake compared to large parts of YU-L7.

Finally, the conditions that promote the weathering processes also lead to the gradual degradation of the OM, as can be indicated by the high significant correlations of the weathering indices with the TOC values (e.g., r = −0.74, p < 0.001 for FENG versus TOC, see Supplementary Tables 3, 4). However, analysis of close-by sediment cores drilled on dry land at the Yukechi site point against strong post-sedimentary decomposition of OM but rather for the input of organically poor and pre-decomposed material (Windirsch et al., 2020). Based on our data, a combination of the two can be assumed. Relatively few OM was deposited, which could then be decomposed comparatively quickly. Despite the low OC content of the Yukechi Yedoma IC sediments, Jongejans et al. (2021a) found substantial greenhouse gas production from sediments upon thawing. Their findings showed that the OM quality and turnover history are the main driver for GHG production (Jongejans et al., 2021a).

When comparing the two cores, it must be taken into account that the Yedoma lake core YU-L15 is significantly sandier than the alas lake core YU-L7 (Table 2; p < 0.001). Usually, there are more adsorption surfaces in silty substrates, depending on the element considered, but this could explain the generally higher element concentrations, and at least higher TOC content, in YU-L7 (p < 0.01). Thus, the comparison of weathering indices and elemental ratios is more suitable to compare different substrates than pure single element contents.

Previous studies from NE Siberia and Central Yakutia showed that Yedoma IC deposits are not only rich in carbonates (Schirrmeister et al., 2011) but also sulfides (Siegert, 1979; Siegert, 1987; Biskaborn et al., 2012; Ulrich et al., 2017a; Ulrich et al., 2019), the latter are common in the form of pyrite (FeS₂). Pyrite forms under reducing conditions and oxidizes relatively easily under oxygen-rich conditions. When the surface thaws and comes into contact with oxygen, sulfide oxidizes to sulfate and can then be easily removed, together with Fe, which is dissolved. Thus, individual surfaces can be detected that represent these short-term thawing processes and subsequent re-freezing (in case of the YU-L15 core) by increasing S and Fe values. Those are in particular prominent in the lower and upper parts of both cores (Supplementary Tables 1, 2). Furthermore, S naturally shows a close relation to the organic matrix, so that high S contents could be detected in the core sections with high TOC concentrations. Kokelj et al. (2013) reported S releases in the form of SO₄ ^2- during permafrost thaw and carbonate weathering from Late Pleistocene ice-rich glaciogenic sediments under oxidizing conditions on the Peel Plateau of northwestern Canada. The carbonate weathering includes sulfide oxidation, which leads to an increase in solute fluxes (Zolkos et al., 2018). With Cl⁻ as an ion that showed only minor increase with permafrost degradation, Kokelj et al. (2013) have shown that the SO₄ ²⁻/Cl⁻ ratios have more than doubled since the 1960s in a Canadian Arctic river. This was resulting from an increase in retrogressive thaw slump (RTS) activities within its watershed due to increasing air temperatures. The following statements are possible if these findings are applied to our results. Since chloride is very easily soluble, it is removed faster than sulfate, which is only generated from sulfide through oxidation processes. The S/Cl ratio in the thawed core YU-L7 and within the talik of YU-L15 is thus higher than in the frozen part of YU-L15 (Table 2), where Cl is bound by the frozen state and has not been removed, yet. Therefore, individual maxima of the S/Cl ratio within YU-L15 suggest also brief thawing processes during pausing sedimentation. The syngenetically rising permafrost table has preserved the resulting S/Cl ratio.

The Ti/Ca ratio is often used as the primary weathering index, which increases as the degree of weathering increases, since the immobile Ti is more stable against weathering than the Ca, which is often bound in CaCO₃ (Fischer et al., 2012). Since we used the reverse Ca/Ti ratio, it drops when the weathering processes increases. When comparing the two cores, it becomes obvious that YU-L7 shows overall lower Ca/Ti values than YU-L15 (p < 0.01), which indicates a higher overall degree of mineral weathering in the alas lake core than in the Yedoma lake core. This is likely caused by the fact that YU-L7 is consistently thawed and is therefore subject to stronger and longer weathering than the frozen core part of YU-L15. Moreover, Ti is generally known to reflect detrital input (usually in lakes) as it is released from source rocks or deposits by physical weathering and minerals containing Ti are not sensitive to dissolution (Demory et al., 2005). Haberzettl et al. (2007) used the Ca/Ti ratio to interpret hydrological variability in lacustrine deposits from Southern Patagonia with the immobile Ti as a proxy for allochthonous clastic input. Using this method, we can differentiate between tendencies of lithogenic to carbonate deposition, with lower Ca/Ti values reflecting an intensified allochthonous input and higher values for autochthonous mineral precipitation.

The ratio of (Sr/Ca)*10⁴ is used to separate lithogenic and biogenic Ca, whereby high ratios indicate lithogenic carbonate, while small ratios reflect the biogenic influence (mussel and snail shells). Overall, the Sr/Ca ratio of the alas lake core was showing much smaller values compared to the Yedoma lake core (Table 2; p < 0.1). This is suggesting a higher biogenic influence on the carbonate amount in the sediments below the alas lake and a predominantly lithogenic (clastic) carbonate input in the Yedoma deposits (Wien et al., 2005; Demina et al., 2019). Beside the predominantly lithogenic carbonate input in the Yedoma deposits, as supposed above, this ratio drops in Y-SU-III (1,672 cm bls) suggesting a short-term biogenic Ca input in this sediment unit. Moreover, if the (Sr/Ca)*10⁴ ratio is >100, it is suggested to be an indication of the presence of aragonite (Sr-rich) rather than calcite (Wegwerth et al., 2014). Aragonite is mostly formed secondary by chemical precipitation from calcareous water. The prerequisite for this chemical precipitation is an increase in the CaCO₃ concentration until a supersaturated solution is reached (Wegwerth et al., 2014; Lampe et al., 2016). Aragonite precipitation is generally more characteristic of marine environments, but Roeser et al. (2016) were able to detect the accumulation of aragonite in lake sediments from NW Turkey in connection with very shallow lake levels and phases of increasing temperatures during the late Pleistocene (i.e. MIS3) and the early Holocene. The CaCO₃ precipitation from calcareous (ground-)water can only occur in thawed subsurface sediments as only then, a subsequent weathering process can be initiated. Thus, the Sr/Ca ratio, if compared with the FENG index, shows high correlations for YU-L7 and YU-L15 of r = 0.83 and r = 0.91, respectively. In this way, different phases of weathering of different degrees can be identified.

The Si/Al and Zr/Al ratios provide information about grain-size changes and can be used to estimate transport energy in sedimentary catchments (Lopez et al., 2006; Urbat, 2009; Bertrand et al., 2012). The amount of Si and Zr in surface sediments is controlled by their association with heavy minerals and are interpreted as signals of detrital input. However, Al is relatively insensitive to changes in the nature of the sediment sources, the size of the catchment area and hydrodynamic processes. Therefore, it represents the ideal element to interpret the lithogenic fraction of the sediment and to use it to normalize other lithophile elements (Bertrand et al., 2012). Increasing Si/Al and Zr/Al ratios during cold to moderate climate conditions may, therefore, reflect higher transport energy and terrestrial sediment supply, since silicate minerals cannot be weathered and removed as much under such conditions (Lopez et al., 2006; Urbat, 2009). The perennially frozen conditions below 1,250 cm bls in the YU-L15 core are reflected in overall high Si/Al ratios (Table 2; Figures 5, 7). The Zr/Al ratio is understood in a similar context indicating cooler phases by high values. Bertrand et al. (2012) demonstrate that the Zr/Al ratio is in particular sensitive in proximal sediment environments of fjords in Northern Chilean Patagonia and increases with increasing transport energy. Especially in lower core parts of YU-L15 (below 1700 cm bls), the Zr/Al values are higher on average (almost twice as high), than in YU-L7 and the talik section of YU-L15. The relatively immobile Zr was only released and removed during increasing weathering processes (Urbat, 2009; Bertrand et al., 2012). Within both cores, the Zr/Al ratio decreases upwards, which suggests warmer conditions with longer lasting thawing cycles. Lower Zr/Al values (in particular within YU-L15), can be interpreted as old sediment surfaces, where higher summer temperatures could have favored deeper active layers.

We consider similar chemical composition and hardly any post-sedimentary modification of the Yedoma deposits parent material, which have, as commonly acknowledged, rather local sources from where different transport mechanisms enabled the polygenetic Yedoma formation (e.g., Schirrmeister et al., 2013). This has been shown in Schirrmeister et al. (2020) by comparing grain-size distributions and endmembers of Yedoma deposits, which revealed a localized pattern and short-distance transport. As we studied a single Yedoma site, i.e. Yukechi, we assume the same general Yedoma formation conditions here and thus a rather homogeneous chemical composition. We further assume that alas deposits formed within and from Yedoma deposits, and are largely composed of taberal and re-located Yedoma deposits, although authigenic deposition during thermokarst lake stages is expected. Thus, any further diagenetic evolution of the material as reflected by increased weathering and expressed by the values of the weathering indices refers to the same base. Although variations between the single units of each core occur (see Figures 5, 7) that allowed us to differentiate our proxies between the sediment units, with the still frozen sediments of YU-L15 as reference.

The applied weathering indices clearly indicate alternating freezing and thawing processes in the study area. Particularly during thawing phases, intensive chemical weathering processes took place due to the presence of (melt) water. Thus, all weathering indices indicated specific periods within both sediment sequences when weathering processes could have been more active. The good accordance of the indices is shown by their significant and high correlation coefficients in both cores (p < 0.001; see Supplementary Tables 3, 4). The correlations between the indices and the sand content, the mean grain size, the Si/Al, and the Sr/Ca ratio, suggests that weathering processes have mainly affected the lithogenic fraction of the sediments and thus the least pre-weathered portion. This would confirm the assumption that the parental material of the Yedoma IC deposits were hardly or not at all pre-weathered. Comparatively short transport routes and primary sedimentation under cold climatic conditions, i.e. conditions that are unfavorable for mineral weathering, could be reasons for this (Schirrmeister et al., 2020). However, these correlations are lower for the YU-L7 core, in particular with the MIA_(R) index. This is probably due to the redox sensitivity of the MIA_(R) index (Babechuk et al., 2014) and sustained reducing conditions in the talik zone. Contemporary sustained reducing conditions are found especially for the upper core sediments that are located immediately below the lakes. The significant negative correlation of the MIA_(R) and FENG indices (positive for ICV; p < 0.001) with the GW/IC in YU-L15 confirms that (chemical) weathering processes stagnate under frozen conditions. The application of weathering indices can therefore be transferred very well to periglacial conditions in order to estimate past freezing processes.

The sampling of YU-L7 was not as detailed as for YU-L15 due to core loss. Nevertheless, homogeneous weathering conditions can be inferred from the relatively constant values of the indices. The uniform ICV values between 0.6 and 0.8 suggest only moderate mineral conversion for the YU-L7 sediments (Cox et al., 1995). Because these sediments were thawed completely, similar weathering processes likely took place throughout the YU-L7 sequence. At 1,762 cm bls and below about 1,900 cm bls, clear inhibited weathering processes are indicated by decisive indices. One explanation for this is the fluctuating expansion of the talik below the lake in connection with changing lake levels. Strongly changing lake levels including phases of complete drying out are typical for Central Yakutian alas lakes (Soloviev, 1959; Soloviev, 1973; Ulrich et al., 2017b). It is generally assumed that with loss of water bodies, taliks within the permafrost re-freeze relatively quickly from the sides (Grosse et al., 2013).

In general, the core YU-L15 is showing a much more differentiated picture for all utilized indices. Since we assume a largely homogeneous chemical composition of the parent material of both sediment sequences, we see the frozen sediments as a reference in our study on which we test our proxies comparatively. Most noticeable within YU-L15 is the significant jump in the index values within Y-SU-III. The frozen core section can clearly be structured from bottom up into three phases of weathering of different degrees (first: below about 1,700 cm bls; second: about 1,700–1,600 cm bls; and third: above about 1,600 cm bls). Starting from the lowest first phase, the weathering tendency increases to the third phase with a low during second phase. This suggests that weathering conditions have increased gradually from the bottom to top of the sediment sequence, which is associated with frequent changes of long-lasting thawing and subsequent freezing processes. Furthermore, we assume a deepening active layer caused by increasing ground temperatures, since the detected chemical weathering can take place more intensively only at higher positive temperatures and thawing permafrost (Zolkos et al., 2018).

Finally, our findings are supported by the Fe₂O₃/MnO ratio. The ratio of iron to manganese in sediments is known to change with the degree of oxidation (usually in lacustrine environments) during accumulation (Davison, 1993; Biskaborn et al., 2013). During reducing soil/sediment conditions, the Fe₂O₃/MnO ratio increases, since Mn is more soluble than Fe and is therefore faster mobilized. For YU-L7, the ratio is always above 50, which suggests sustained reducing conditions throughout the core (Lampe et al., 2016). For the frozen part of core YU-L15, a clear difference was obvious between the first and third weathering phases. Below about 1,700 cm bls, the ratio is significantly smaller than above, which in turn confirms the discussions of the FENG and MIA_(R) value changes. It finally suggests that due to the long-lasting frozen state of the sediment units Y-SU-I and Y-SU-II of the Yedoma lake core, Mn was not reduced and discharged, which is reflected in lower Fe₂O₃/MnO ratios. The upper part of the frozen sediments has already been subject to several thawing phases and higher water availability due to rising temperatures, which is reflected in higher ratios. Since reducing anaerobic soil conditions significantly slows the degradation of organic matter (Knoblauch et al., 2012), higher TOC values for the areas of higher Fe₂O₃/MnO ratio could be expected. However, this correlation is not significant for the sediment cores presented here and ultimately cannot answer the question of the generally low TOC contents. In addition, higher TOC levels in the active layer can also be a sign of increasing plant production during warmer periods in the past (Schirrmeister et al., 2011).

Our dating results do not allow a detailed age-depth model. Nevertheless, they suggest that the sediment sequence under the Yedoma lake was deposited quickly during MIS3. In contrast, the calibrated radiocarbon age range of the sediments under the alas lake points to a much longer sedimentation history, which extends from the MIS3 into the middle Holocene (Table 1). However, similar geochemical properties of both cores allow us to conclude that both sediment sequences originate from the same parent material and a primarily terrestrial sedimentation history. Evidence of higher terrestrial plants found throughout both cores by Jongejans et al. (2021b), indicated by the dominance of specific long-chained alkanes, corroborates the lack of lacustrine sediments even under the alas lake. The lenticular, layered and suspended cryostructures in the middle and lower section of the YU-L15 core (Figure 10) point to a gradual syngenetic freezing of the active layer of subaerial deposits, crack infilling and different rates of sedimentation (French and Shur, 2010).

Our PCA results are particularly helpful to discuss the general sedimentation environment of the Yedoma IC deposits at the Yukechi study site. In both cores, PC1 is generally determined by organic enrichment in connection with clayey-silty grain sizes on the one hand (negative sample scores) and sand-rich grain sizes, higher transport energy of sili-clastic material and conditions promoting weathering on the other hand (positive sample scores). The different distribution of the variables along PC2, however, reflects the differences between the two cores with regard to their depth-dependent thawed and/or frozen states and related redox conditions (Figures 5, 7, 9). We found that the sedimentation conditions gradually change from bottom to top, but return to conditions in the top most recent stratigraphic units comparable to those of the bottom sediment units.

For the East Siberian lowland region, Schirrmeister et al. (2013), Schirrmeister et al. (2020) developed the concept of a polygenetic Yedoma IC accumulation (including alluvial, fluvial and partly aeolian transportation) from nearby sources and syngenetic ground ice segregation, ice-wedge growth, peat aggradation, cryosol formation and cryoturbation. Furthermore, many paleo-environmental studies from NE Siberia discuss the MIS3 interstadial (57–29 ka BP), i.e. the main Yedoma IC formation period, as a climate optimum stage that was punctuated by strong climate fluctuations (Anderson and Lozhkin, 2001; Zech et al., 2008; Wetterich et al., 2014; Diekmann et al., 2017). For instance, plant macrofossils studied by Kienast et al. (2005) on Bykovsky Peninsula in the Central Laptev Sea region suggest relatively warm and moist summers even for Northern Yakutian regions with mean July temperatures above 12°C. Comparable variations in summer climate conditions in Central Yakutia could be the reason for the strong fluctuations in our proxy data, especially in the YU-L15 core. Warm and moist phases caused increased fluvial-alluvial sedimentation of sand-rich sediments, deeper active layer depths, increased mineral weathering and higher OM decomposition. In contrast, cooler and drier phases are dominated by aeolian sedimentation, inhibited weathering processes, and hampered decomposition of OM (Strauss et al., 2013; Jongejans et al., 2021b). The clearly unimodal GSDs with a high peak in the coarse silt to fine sand fraction (Figures 4, 6) and the lack of stratification within both sediment sequences indicate the high aeolian component in our sediment cores (Péwé and Journaux, 1983; Murton et al., 2015). However, Schirrmeister et al. (2020) point out that the formation of finer-grained sediments in Yedoma deposits can also be the result of severe in situ frost weathering. But, high proportion of sand in the middle sediment units of both cores and especially in Y-SU-IV of the YU-L15 core (maximum mean grain size of 260 µm) can only be explained by fluvial-alluvial transport processes (Schirrmeister et al., 2020; Windirsch et al., 2020). These were intense short events of fluvial-alluvial deposition of large sediment amounts, which we suggest to be comparable to the short time frame covered by the Yedoma lake core and the proxy cyclicity represented by the sample score variations of PC1. A representative exception in the sedimentation history of the YU-L15 core is the unit Y-SU-III, which was defined as the second weathering phase, with low indices but exceptionally high TOC values and Fe₂O₃/MnO ratios. All proxies in this unit indicate that short-term thermokarst processes could have taken effect here (Figure 7) and that a shallow thermokarst lake may have developed and disappeared during a very short period (Ulrich et al., 2017a). A similar proxy pattern can be seen in the A-SU-II of the alas lake core at a depth of 1,763 cm bls. However, it cannot be conclusively clarified whether these short-term processes are the same in terms of time and thaw intensity.

For the transition to the MIS2 and the Last Glacial Maximum (LGM, 26.5–19 ka BP; Clark et al., 2009), significant summer and winter cooling and arid conditions are often described for NE Siberia (Wetterich et al., 2011). Such conditions are supposed to be typical for the formation of the NE Siberian Yedoma IC deposits (Schirrmeister et al., 2013; Schirrmeister et al., 2020). We cannot identify the MIS2 based on our dating and proxy data. A Holocene cover of the Yedoma IC deposits at the Yukechi study site is also not at all visible in the YU-L15 core and even difficult to identify in the YU-L7 core. Wetterich et al. (2020) describe several hiatuses for Northeast Siberian Yedoma deposits during MIS3-2 as typical regional disturbance leading to either low accumulation and/or fluvially-triggered or climate-triggered erosion of deposited IC material. We assume that for the Yukechi study site, climate-related thawing of the permafrost and successive changes of accumulation areas during the Lateglacial to Holocene transition and/or the HTM are the reason for the lack of younger deposits.

Various studies suggest that large alas basins in Central Yakutia formed during the pronounced climate warming after the LGM and during the HTM as a result of extensive thawing processes, ground-ice melting, and subsidence of Yedoma deposits (Soloviev, 1973; Ulrich et al., 2017a; Ulrich et al., 2019). During the formation of these thermokarst basins, younger OM likely was also deposited (e.g., mobilized from thermo-erosion processes and basin expansion) as seen by the two uppermost early and mid-Holocene dates from the YU-L7 core. In contrast, the inverse age at 1,560 cm bls of YU-L7 (Table 2) is interpreted as the result of old OM incorporation already during Yedoma deposition or thermokarst lake and talik development (Gaglioti et al., 2014). The calibrated radiocarbon dates from the uppermost core sample in YU-L15 (142 ± 49 cal. yrs. BP) suggest the deposition of modern sediments during the formation of the young thermokarst lake above, but could also reflect fresh input prior to the lake phase and cryoturbation.

The questions regarding the significantly lower sand proportion in the YU-L7 core remains. Biskaborn et al. (2013) and Ulrich et al. (2019) discussed that depositional activities and grain size fractionation during thermokarst (lake) basin development might have overprinted the original sediment characteristics of the Yedoma source sediments. Concerning this, distal sediment transport after lake bank erosion or a low-energy sedimentation milieu within an expanding thermokarst lake basin could have led to a homogeneous fractionation of finer-grained sediments in a basin center as represented by the YU-L7 core site.