Carbon Dioxide and Methane Release Following Abrupt Thaw of Pleistocene Permafrost Deposits in Arctic Siberia

The decomposition of thawing permafrost organic matter (OM) to the greenhouse gases (GHG) carbon dioxide (CO2) and methane forms a positive feedback to global climate change. Data on in situ GHG fluxes from thawing permafrost OM are scarce and OM degradability is largely unknown, causing high uncertainties in the permafrost‐carbon climate feedback. We combined in situ CO2 and methane flux measurements at an abrupt permafrost thaw feature with laboratory incubations and dynamic modeling to quantify annual CO2 release from thawing permafrost OM, estimate its in situ degradability and evaluate the explanatory power of incubation experiments. In July 2016 and 2019, CO2 fluxes ranged between 0.24 and 2.6 g CO2‐C m−2 d−1. Methane fluxes were low, which coincided with the absence of active methanogens in the Pleistocene permafrost. CO2 fluxes were lower three years after initial thaw after normalizing these fluxes to thawed carbon, indicating the depletion of labile carbon. Higher CO2 fluxes from thawing Pleistocene permafrost than from Holocene permafrost indicate OM preservation for millennia and give evidence that microbial activity in the permafrost was not substantial. Short‐term incubations overestimated in situ CO2 fluxes but underestimated methane fluxes. Two independent models simulated median annual CO2 fluxes of 160 and 184 g CO2‐C m−2 from the thaw slump, which include 25%–31% CO2 emissions during winter. Annual CO2 fluxes represent 0.8% of the carbon pool thawed in the surface soil. Our results demonstrate the potential of abrupt thaw processes to transform the tundra from carbon neutral into a substantial GHG source.

This large range of model predictions is related to fundamental uncertainties about the decomposability of permafrost OM. Current estimates on how fast permafrost OM may be mineralized to CO 2 after thaw range between 15% in 100 years and 66% in one summer thaw season (Elberling et al., 2013;Knoblauch et al., 2013;Plaza et al., 2019;Schädel et al., 2014;Schuur et al., 2015;Vonk et al., 2012). While some biomarker and fractionation studies indicate its rapid decomposability (Mueller et al., 2015;Strauss et al., 2015), longterm laboratory incubations demonstrate the dominance of slowly degradable OM with turnover times of hundreds of years Schädel et al., 2014). Results from laboratory incubations are subsequently used to constrain simulations of future GHG release from permafrost landscapes (Koven et al., 2015;Schneider von Deimling et al., 2015). However, laboratory incubations may not consider all of the environmental parameters that strongly affect in situ OM turnover, such as gradients of temperature, oxygen and moisture, freezing and thawing or the input of fresh OM Schädel et al., 2016;Walz et al., 2017;Wang et al., 2020;Wild et al., 2016). Thus, it still remains unclear how far organic carbon decomposition rates from laboratory studies represent those under in situ conditions.
Most studies on GHG fluxes from thawing permafrost landscapes report ecosystem respiration (R eco ) fluxes from vegetated soils affected by thawing permafrost (Abbott & Jones, 2015;Fouche et al., 2014;Mauritz et al., 2017;Natali et al., 2015;Pegoraro et al., 2021;Ravn et al., 2020). These CO 2 fluxes originate from both microbial decomposition of soil OM [heterotrophic respiration (R h )] and plant respiration [autotrophic respiration (R a )]. Furthermore, R h fluxes include both CO 2 from OM decomposition of recently fixed carbon in the surface active layer and from older permafrost carbon thawing at the soil's bottom . A direct measurement of in situ turnover of recently thawed permafrost OM is challenging since permafrost thaws in typically 0.3-1 m depth at the bottom of the active layer.
The few studies that estimated GHG fluxes from permafrost OM used a combination of field manipulations, incubation studies and 14 C analysis (Cooper et al., 2017;Estop-Aragonés et al., 2018;Hicks Pries et al., 2016;Schuur et al., 2009) but in situ measurements of permafrost OM decomposition are still rare. Since the degradability of deep permafrost OM may substantially differ from those of the surface active layer (Cooper et al., 2017;De Baets et al., 2016;Schädel et al., 2014;Song et al., 2020), knowledge on the degradability of permafrost OM under in situ conditions is needed to constrain future GHG fluxes from thawing permafrost landscapes.
The aims of the current study were to (a) quantify the in situ decomposability of Pleistocene and Holocene permafrost OM, (b) evaluate how far data from short-and long-term laboratory incubations may be used to predict in situ CO 2 and methane fluxes, and (c) quantify the amount of permafrost OM released as CO 2 during one year. Therefore, we measured in situ R h fluxes of CO 2 and methane in July 2016 and July 2019 from different soils in an active thaw slump of the northern Siberian Lena River Delta and related these to the thawed carbon pools in the active layer. In parallel, we incubated soil samples from the same sites and used the resulting CO 2 and methane production rates to calculate potential GHG fluxes. Finally, we calibrated a regression model with in situ flux observations, and a dynamic model with long-term incubation data from samples of the same permafrost deposits, to estimate annual CO 2 production from the thawed permafrost and to determine the fraction of permafrost OM decomposed over a whole year.

Study Sites
Field measurements and sampling were conducted at an active thaw slump on the island Kurungnakh (72.339°N, 126.292°E, Figure 1) in the Lena River Delta, North Siberia. The investigation area is situated in the zone of continuous permafrost with a trans-arctic, continental climate. The annual mean air temperatures on the neighboring island Samoylov was −12.3°C with annual soil temperatures ranging between about 10°C and −28°C at 6 cm depth (Boike et al., 2019). The soils are frozen for about eight months of the year. Only during the short summer season, a shallow surface layer thaws to a depth of less than 70 cm. Kurungnakh is comprised of Yedoma deposits (17-34 m above river level) that developed during the middle and late Weichselian (MIS 3 and MIS 2) and contain huge syngenetic ice-wedges. The OM frozen in these deposits originate mainly from a grass-dominated tundra-steppe vegetation. The Yedoma is overlain by an about 5 m thick cover of Holocene deposits. The mineral soil is mainly comprised of sandy to silty loam with intermediate organic carbon contents (2%-7%) Wetterich et al., 2008). Due to their high ice content, Yedoma deposits are strongly affected by thermokarst processes and particularly vulnerable to abrupt thaw (Czudek & Demek, 2017). After thawing of massive ice-wedges exposed at the bluffs of the Lena River channels, retrogressive thaw slumps develop. Such retrogressive thaw slumps are a widespread form of thermokarst and one of the most rapidly developing erosional features (Bröder et al., 2021;Cassidy et al., 2016;Tanski et al., 2017). In front of the thaw slump's headwall, soil material is relocated and forms a mud-slurry, which is then transported downslope. Due to rapid erosion, the surface of an active thaw slump is bare of recent vegetation and only R h contributes to in situ CO 2 and methane fluxes, which makes these features ideal systems for studying the in situ decomposition of recently thawed permafrost OM.
We measured in situ gas fluxes in July 2016 at seven field sites in such an active thaw slump (Table 1 and Figure S1 in Supporting Information S1). The sites were located in an area of about 0.5 ha with a distance of about 30 m to each other and were grouped according to their parent soil material. Five sites were sampled at the thaw slump floor (SF) and consisted of relocated soil material, which mainly comprised of recently thawed permafrost but also some surface mineral soil material and remnants of surface vegetation from active layer detachment at the front of the headwall (Figures 1c-1e). Further two sites were sampled at two thermokarst mounds (TM) that were not affected by relocated surface material and consisted of mineral soil of the former polygonal centers between the ice wedges (Figures 1b and 1f). In July 2019, two SF sites (SF3 and SF6) and the two TM sites were resampled.

Field Measurements and Sampling
CO 2 and methane fluxes were measured about every other day between July 12, 2016 and July 28, 2016 and between July 10, 2019 and July 23, 2019 with opaque emission chambers (diameter 24 cm, 20 cm height). Three soil collars were pre-installed at each site at least 24 hr before the first measurement. At each sampling day, three emission measurements were conducted above these three different soil collars per site. CO 2 concentrations inside the chambers were quantified over a period of 5 min by portable GHG analyzers  using a Li-840a (LI-COR Biosciences, Lincoln, USA) during the first week of measurements and an UGGA (Los Gatos Research, San Jose, USA) during the second week in 2016. Methane emissions were measured during the first week by collecting gas samples over a period of 30 min from the closed chambers that were analyzed the same day by gas chromatography (Preuss et al., 2013). During the second week, methane concentrations were quantified with the UGGA in parallel to CO 2 concentration measurements. Both analyzers were calibrated at the field station on Samoylov with an external calibration gas (403 ppm CO 2 , 3.1 ppm CH 4 ). In parallel to each gas flux measurement, thaw depth, air temperature and soil temperature gradients (5-cm increments between surface and thaw depth) were measured. In July 2019, CO 2 and methane fluxes were measured with an UGGA, air and soil temperatures and thaw depth as described above.
Soil samples were collected in 2016 and 2019 between the soil surface and the thaw depth using steel tubes (inner diameter 5.8 cm). In 2016, one profile was sampled at each sampling site at the beginning of the measurement period (July 10, 2016). Two soil samples were retrieved per site, the first between the surface and 20 cm depth and the second between 20 cm and the thaw depth (Table 1). In 2019, three different profiles were sampled per sampling site with a 10 cm depth increment between the surface and the thaw depth in the middle of the measurement period (July 18, 2019). All soil samples were quantitatively retrieved from the sample tube and the bulk density was calculated from sample volume and dry weight.

Soil Analysis
Total soil carbon and nitrogen were quantified with an elemental analyzer (VarioMAX Elementar Analysensysteme GmbH, Hanau, Germany) after the soil had been sieved (<2 mm), milled, and dried at 105°C. Total inorganic carbon (TIC) was quantified from the amount of CO 2 released after sample treatment with phosphoric acid. Total organic carbon (TOC) was calculated as the difference between total carbon and TIC. Soil pH was measured in a suspension of 10 g of dried soil in 25 mL of demineralized water. Soil water content was calculated from the weight difference between fresh samples and those dried at 105°C. The 14 C-age of bulk TOC was analyzed with a MICADAS AMS system (Ionplus, Dietikon, Switzerland) (Steinhof et al., 2017). AMS 14 C data were calibrated using Oxcal 4.3 and the INtCal13 calibration curve (Bronk Ramsey, 2009).

Incubation Experiments
Short-term laboratory incubations were conducted with samples from sites SF1, SF3, SF4, SF6 and TM2 (Table 1) in July 2016. Samples were processed at the day of sample collection and incubated at the research station on Samoylov Island. For anaerobic incubations, samples were weighted into glass flasks (100 mL) in a glove bag under nitrogen. The bottles were closed with rubber stoppers and the headspace was repeatedly evacuated and flushed with nitrogen. Aerobic incubations were prepared under ambient air in 250-mL bottles. All samples were incubated in three replicates for 18 days at 4°C. The amount of CO 2 and methane was measured every day (aerobic incubations) or every other day (anaerobic incubations) with a gas chromatograph . CO 2 and methane production rates were calculated from cumulative gas production during the last two weeks of the incubations, to omit initial disturbances during sample preparation. Samples from the surface soil of TM2 were incubated at 4°C for two years and reanalyzed in July 2018.

Calculation of Carbon Dioxide and Methane Fluxes
Methane and CO 2 emissions from bare soils were calculated by a linear regression of gas concentrations inside the emission chambers as described by Eckhardt and Kutzbach (2016). Gas emissions were normalized to surface area (f A ).
Field CO 2 fluxes were also normalized to the amount of organic carbon thawed in the active layer at the day of gas flux measurements, using the carbon concentrations in the different sampled soil layers according to: with C A = Carbon pool in the active layer, ρ i = dry bulk density of soil layer i, d i = depth of soil layer i and [TOC] i = TOC concentration in layer i.
The CO 2 or methane flux normalized to thawed carbon during each of the emission measurement (f C ) was calculated according to: To evaluate how far results from short-term laboratory incubations represent GHG production from thawing permafrost OM under in situ conditions, CO 2 and methane production rates from laboratory incubations were upscaled to a potential areal flux (f p ), which was calculated for each day of field flux measurements according to: with r i = aerobic or anaerobic CO 2 or methane production rate determined at 4°C in samples incubated from layer i, T i = mean temperature of layer i during the day of field CO 2 and methane flux measurement. Q 10 is the factor by which a reaction rates increases after a temperature rise of 10°C. A Q 10 value of 2 was applied, which is based on permafrost incubation studies (Schädel et al., 2016;Vaughn & Torn, 2019).

Soil Temperature Modeling
To extrapolate summer CO 2 emissions to an annual CO 2 flux, a year-round time series of soil profile temperature data is required. This time series was derived for 5 cm layers until 1 m soil depth using the land surface model JSBACH Porada et al., 2016) with a specific setup for the Lena Delta (Chadburn et al., 2017;Ekici et al., 2015) but without considering any vegetation insulation layer (e.g., mosses) as the sampling sites were not vegetated. Only for model evaluation, we additionally compared a model run with total moss cover to the soil temperature data from Samoylov Island nearby. The meteorological forcing data is based on the CRUNCEP dataset (Viovy, 2018). This time series was bias-corrected to observations from the Samoylov station (Boike et al., 2019) using an overlap period 2003-2016 and an additive method except for precipitation and wind speed, for which a multiplicative technique was applied (Beer et al., 2014). Then, climate data were replaced by observed meteorological data from the Samoylov station (Boike et al., 2019) during 2002-2016. Winter precipitation was approximated by positive day-to-day snow depth changes. Soil water, ice content and temperature were initialized using a classical spin-up procedure following Beer et al. (2018), followed by a transient model run from 1901 to 2016. Soil temperature results with and without moss layer were evaluated with observations at the Samoylov station (with moss layer) and by direct thaw slump measurements during July 2016 (without moss layer).

Annual Carbon Dioxide Fluxes of the Thaw Slump
The annual CO 2 fluxes from heterotrophic respiration averaged over the different thaw slump sites were estimated using two different modeling strategies and by using different observations. First, a Q 10 model was calibrated against July 2016 chamber measurements of CO 2 emissions and observed average topsoil (0-20 cm) temperature. Subsequently, the Q 10 model was applied to year-round topsoil temperature estimates from a site-level run of the land surface scheme JSBACH (Section 2.6). Second, a first-order kinetics model (ICBM) with two connected reservoirs (Andrén & Kätterer, 1997) was applied. The parameters at a reference temperature of 4°C were calibrated using long-term incubation data from either Holocene or Pleistocene deposits nearby . The decomposition rate constants are modified by a Q 10 model using a Q 10 value of 2 (Schädel et al., 2016;Vaughn & Torn, 2019). The ICBM was driven by soil temperature data (20 five-cm-layers) from the site-level JSBACH model run (Section 2.6) but with observed soil OM content (Table 1) for topsoil and subsoil layers from each site. This approach led to several results per thaw slump site due to a variation of calibrated decomposition rate constants and results in a distribution of CO 2 fluxes from which we present the median and the upper and lower quartiles. Model outputs of the ICBM were considered at temperatures down to −10°C since below this temperature, microbial activity was shown to be practically negligible (e.g., Mikan et al., 2002;Natali et al., 2019).

Quantification of Methanogenic Archaea
Total genomic DNA was extracted in duplicates using the DNeasy Power Soil Kit (Qiagen, Hilden, Germany) according to the manufacturer's protocol. Quality and quantity of the DNA was assessed through gel electrophoresis and photometry using the Qubit2 system (Invitrogen, HS-quant DNA). DNA was purified using the 15HiYield PCR Clean-Up & Gel-Extraction Kit (SLG, Gauting, Germany) to reduce PCR inhibitors prior to PCR applications. The enumeration of methanogenic gene copies was realized through quantitative PCR (qPCR) targeting the mcrA gene of methanogens as described by Liebner et al. (2015). Each of the duplicate DNA extractions per sample was amplified in triplicates.

Statistics
Data were tested for normal distribution by the Shapiro-Wilk test. Since data transformation (log10 and square root) of non-normally distributed data did not improve data distribution, statistical tests were performed with non-transformed data. In case of normal distribution, differences between mean values were compared by a two-sided independent t-test and correlations between variables by Pearson's correlation coefficient. In case of nonnormal distribution, non-parametric tests were performed, for two groups a Mann-Whitney U test and for more than two groups an independent samples Kruskal-Wallis test followed by a Dunn-Bonferoni post-hoc test. A Spearman's rho correlation analysis with a two tailed significance test was applied to test correlation between non-normal distributed data. Significance was tested on a level of p < 0.05. All statistical analyses were conducted with IBM SPSS Statistics 27 (IBM Corporation, Armonk, USA).

Soil Properties
All soils were classified as Protic Cryosols due to the presence of permafrost and a lack in soil horizon development (WRB, 2014) ( Figure 1). Neutral pH values dominated with extreme values between moderately acidic at SF5 to moderately alkaline at TM2 (Table 1). Organic carbon concentrations ranged between 1.2% and 6.1% and C/N ratios between 11 and 21 and both variables showed no clear change with depth. The average TOC concentration of SF soils (3.6 ± 1.2% (mean ± SD) in 2016, 3.8 ± 1.2% in 2019) was substantially higher than in TM soils (2.4 ± 0.4% in 2016, 2.6 ± 1.0% in 2019) but this difference was not significant (t-test, p > 0.07). Except for one site (TM2), TOC pools in the uppermost 30 cm did not differ between 2016 and 2019 (Table 1). However, due to substantially higher air temperatures during the thaw period, the soil temperature was higher in 2019 than in 2016 (Table S1 in Supporting Information S2) with significant differences at SF6 and TM1 (t-test, p < 0.01). Consequently, both thaw depth and thawed carbon pools were significantly higher (Mann-Whitney U test, p < 0.005) at all sampling sites in 2019 in comparison to 2016. The organic carbon originated from the Pleistocene/Holocene transition (10-13 kyr) at SF1, the site most distant from the headwall, while TOC of the remaining SF sites were of mid to late Holocene age (Table S2 in Supporting Information S2). In contrast, OM at the TM sites originated from the last glacial maximum (TM2) and the Kargin interstadial (TM1) with a maximum age of 33 kyr BP. Only SF soils contained mixedin organic rich material of modern 14 C-age (see SF4 Table S2 in Supporting Information S2). Soil moisture did not clearly differ in July 2019 and July 2016 (Table S3 in Supporting Information S2).

Abundance of Methanogenic Archaea
The abundance of methanogenic archaea was significantly higher (Mann-Whitney U test, p < 0.01) in the SF soils than in TMs (Table 1). Highest abundance of methanogens was found in the bottom samples of SF3 and lowest in the bottom layer of TM2. The number of mcrA genes was below the detection limit of about thousand copies per gram in 61% of the TM soil amplifications while the qPCR was successful in all of the replicates of the SF soil samples. mcrA gene copy numbers correlated significantly with CH 4 production (Spearman's rho = 0.79, p = 0.007) and anaerobic CO 2 production (Spearman's rho = 0.69, p = 0.029) in laboratory incubations but with none of the soil parameters tested (TOC, N, C:N, pH, water content). Furthermore, methanogen abundance in the active layer did not correlate with methane emissions from the respective sampling sites (Spearman's rho = −0.185, p = 0.265). The undisturbed active layer soil above the headwall contained in its bottom mineral horizon (8-20 cm depth) 26.0 ± 20.9 × 10 5 mcrA gene copies per gram.

Carbon Dioxide and Methane Fluxes
Median CO 2 fluxes from the different thaw slump soils in 2016 and 2019 ranged between 0.24 and 2.6 g CO 2 -C m −2 d −1 (Figure 2a) In 2016, CO 2 fluxes per square meter were significantly correlated (Spearman's rho 0.32 to 0.38, p < 0.01) with the average active layer temperature, TOC ( Figure S2a in Supporting Information S1) and total nitrogen in the active layer, while in 2019, CO 2 fluxes significantly correlated only with the amount of TOC and nitrogen in the active layer during the day of measurement.
When normalizing CO 2 fluxes to the amount of unfrozen TOC in the active layer during flux measurements, median fluxes at the different sites ranged between 0.020 and 0.118 g CO 2 -C (kgC) −1 d −1 and 0.015 and 0.067 g CO 2 -C (kgC) −1 d −1 in 2016 and 2019, respectively, with highest CO 2 fluxes from TM1 in both years (Figure 3). Despite substantially higher soil temperatures in 2019 (Table S1 in Supporting Information S2), median fluxes from all sites decreased (Mann-Whitney U test, p < 0.001) from 0.061 g CO 2 -C (kgC) −1 d −1 in 2016 to 0.033 g CO 2 -C (kgC) −1 d −1 in 2019 (Table  S4 in Supporting Information S2). Also the site wise comparison showed  a decrease of CO 2 fluxes at each site with a significant decrease at SF3, TM1 and TM2 (Mann-Whitney U test, p < 0.005). Furthermore, median CO 2 fluxes normalized to thawed carbon did significantly differ (Mann-Whitney U test, p < 0.001) between SF sites (0.053 g CO 2 -C (kgC) −1 d −1 ) and TM sites (0.094 g CO 2 -C (kgC) −1 d −1 ) in 2016 but not in 2019.
In 2016, median methane fluxes ranged between 3.7 and 5.1 mg CH 4 -C m −2 d −1 at the four sites that showed significant methane emissions (SF3, SF4, SF5, and SF6) (Figure 4a). Methane fluxes from the three replicate plots of the same site could differ substantially, e.g., between 132 and 0.8 mg CH 4 -C m −2 d −1 at the same day at SF4. Methane fluxes at SF3 correlated significantly with thaw depth in 2016, and thawed carbon and nitrogen during measurement (Spearman's rho = −0.717, p < 0.001) but no significant correlation was found at any other site with any environmental parameter tested (thaw depth, surface and mean soil temperature, amount of thawed TOC and nitrogen). In 2019, methane fluxes from SF6, the only site releasing considerable amounts of methane, were positively correlated with thaw depth, and thawed carbon and nitrogen (Spearman's rho = 0.488, p = 0.001). Median methane fluxes at SF6 were higher in 2016 (4.8 mg CH 4 -C m −2 d −1 ) than in 2019 (2.8 mg CH 4 -C m −2 d −1 ) and this difference was close to significant (Mann-Whitney U test, p = 0.067) (Figure 4a). At sites with detectable methane emissions, methane contributed between 0.02% and 6.5% (median 0.41% and 0.11% in 2016 and 2019, respectively) to total carbon fluxes to the atmosphere.

Carbon Dioxide and Methane Production in Laboratory Incubations
Microbial CO 2 production rates under aerobic conditions were highest in surface samples of SF3 and lowest in the surface sample of SF1 ( Figure S3 in Supporting Information S1). The median value of aerobic CO 2 production rates from all sampling sites was 0.73 (0.82, 0.41 upper and lower quartile) μmol CO 2 g −1 d −1 . Anaerobic conditions slowed down CO 2 production in all except for one soil layer (TM2, 20-43 cm) by an average factor of 3.8 ± 2.6 (n = 9) ( Figure S3 in Supporting Information S1). Aerobic CO 2 production was most strongly correlated with TOC and total nitrogen (Spearman's rho = 0.69, p < 0.001) and the C:N ratio (Spearman's rho = 0.68, p < 0.001). Anaerobic CO 2 production rates did not correlate with any of these parameters.
Anaerobic methane production was observed in all samples except for those from TM2, but the rates were about three orders of magnitude lower than anaerobic CO 2 production rates (median CO 2 :CH 4 ratio 657 (215, 4,707) (upper and lower quartile)) with highest methane production in the surface samples of SF1 and SF3 and lowest in SF6 ( Figure S2c in Supporting Information S1). However, after two years of incubation at 4°C on Samoylov, methane production could also be detected in the TM2 surface soil with a final mean ratio between total CO 2 and methane production of 5.9 and 6.9 (n = 2).

Potential Gas Fluxes From Laboratory Incubations
Potential CO 2 fluxes (f p ), calculated from CO 2 production rates in aerobic short-term incubations significantly (Kruskall-Wallis, p < 0.01) overestimated in situ CO 2 fluxes by a factor of 1.3-6.4 (mean 4.4 ± 1.4 SD) (Figure 2). Since anaerobic conditions substantially reduced CO 2 production rates, the difference between in situ and potential CO 2 fluxes decreased with significantly higher potential fluxes at sites SF1 and TM2 (Kruskall-Wallis, p < 0.02), similar fluxes at SF4 and significantly lower potential fluxes at SF3 and SF6 (p < 0.01). The mean ratio between in situ fluxes and potential anaerobic fluxes was 1.9 ± 1.5. In situ CO 2 fluxes from all sites correlated both with potential aerobic CO 2 fluxes (Spearman's rho = 0.62, p < 0.001, Figure S2b in Supporting Information S1) and with potential anaerobic CO 2 fluxes (Spearman's rho = 0.28, p < 0.05).
Potential methane fluxes from laboratory incubations were significantly lower (Mann-Whitney-U test, p < 0.002) than in situ emissions at sites SF4 and SF6 (Figure 4b). Since no significant methane emissions could be detected at SF1, potential methane fluxes overestimated in situ emissions at this site while potential and in situ fluxes did not significantly differ at site SF4. At TM2, we detected no methane production in short-term laboratory incubations and also no methane emissions from the soils.

Q 10 Model Calibration
The Q 10 regression model could only be calibrated at sites SF1, SF3, and SF4 ( Figure S4 in Supporting Information S1) since the in situ soil temperature range during flux measurements was too small at the other sites. Q 10 values were very similar at SF1 and SF3 (2.7 and 2.6, respectively) but substantially higher at SF4 (4.1). In contrast, CO 2 fluxes at 8°C were most similar at SF1 and SF4 (0.92 and 0.94 g CO 2 -C m −2 d −1 ) and more than twice as high (2.58 g CO 2 -C m −2 d −1 ) at SF3, although sites SF3 and SF4 have similar OM content, which were higher than at site SF1.

Model Evaluation
The estimated year-round soil temperature dynamics simulated by the site-level run of the JSBACH land surface scheme without a moss layer shows a high agreement with observations at the thaw slump in July 2016 ( Figure 5). Additionally, the JSBACH model results considering a 90% moss cover show a very good agreement with the observations from a vegetated polygon at nearby island Samoylov ( Figure 5). There is a positive temperature bias in December 2016 when comparing the model simulation including a moss layer to the temperature data of the Samoylov soil station. The resulting positive bias of estimated CO 2 emission during December 2016 is about 3% (Q 10 model) and 5% (ICBM).
The CO 2 emissions from the thaw slump in 2016 simulated with the Q 10 model and the ICBM show a high agreement with in situ chamber measurements ( Figure 6). The ICBM results, based on both topsoil and subsoil temperature, show less day-to-day variation than the Q 10 model that purely relates emissions to temperature and therefore translate all  temperature variation immediately into emission variation. However, both model results are similar on average after July 15 (bias = −0.03 g CO 2 -C m −2 d −1 and RMSE = 0.15 g CO 2 -C m −2 d −1 ), although both models are calibrated using completely independent datasets (chamber measurements vs. long-term incubation experiments) and are driven by independent data (topsoil temperature vs. soil OM content and total soil temperature column dynamics).

Annual Carbon Dioxide Emissions From the Thaw Slump
The median cumulative annual carbon release as CO 2 from the thaw slump in 2016 was estimated by the two models as 184 (269,116) gCO 2 -C m −2 (ICBM) and 160 (342, 139) gCO 2 -C m −2 (Q 10 model) (Figure 7a). These are the median and (75, 25)-percentiles from a range of model results using information from the different thaw slump sites (Q 10 model) or a calibrated parameter distribution (ICBM). Since simulated surface soil temperatures were below −10°C between January 1, 2016 and May 24, 2016 ( Figure 5), zero CO 2 production was assumed in this time period by the ICBM. Soil temperatures below 0°C were simulated for the topsoil layer by the JSBACH land surface scheme from of January 1 until May 28 and after September 24, 2016 and the CO 2 fluxes during this period accounted for 25% (ICBM) to 31% (Q 10 model) of annual CO 2 fluxes.
When relating these emissions to the available soil OM thawed in the active layer at the respective day of the year, the ICBM estimates that 0.78 (1.16, 0.52)% of thawed soil OM was mineralized to CO 2 over the whole year 2016 (Figure 7b).

In Situ Carbon Dioxide and Methane Fluxes
Abrupt permafrost thaw is affecting only a minor fraction of the northern permafrost region but GHG fluxes from these areas are expected to almost equal those from the remaining permafrost region until 2300 (Turetsky et al., 2020). Since the studied thaw slump soils were not vegetated, GHG emissions originated only from recently thawed permafrost and do not consider plant autotrophic CO 2 respiration like most previous studies (Abbott & Jones, 2015;Natali et al., 2015;Ravn et al., 2020). The heterotrophic CO 2 fluxes at Kurungnakh (0.24-2.6 g CO 2 -C m −2 d −1 ) were relatively high in comparison to the median CO 2 emissions from soils of the same region, e.g., from similar thawing Yedoma deposits (0.07-0.43 g CO 2 -C m −2 d −1 ) (Vonk et al., 2012), or different polygonal tundra soils (0.24-0.48 g CO 2 -C m −2 d −1 ) (Eckhardt et al., 2019). However, since soil temperature and active layer carbon pools were the main driving parameters for CO 2 exchange fluxes, comparing gas fluxes from different tundra soils without knowledge of latter parameters is problematic. Similar heterotrophic CO 2 fluxes as on Kurungnakh were only measured from active thaw slump soils in north-western Alaska (Jensen et al., 2014) but at substantially higher soil temperatures and active layer depths. The relatively high heterotrophic respiration fluxes from the Kurungnakh thaw slump soils indicate the presence of highly labile OM but laboratory and field incubations indicate that this labile carbon pool comprises only a minor fraction of TOC Moni et al., 2015;Schädel et al., 2014;Zimov et al., 2006). Therefore, it was expected, that CO 2 fluxes would decline over time. This was also the case after normalizing in situ CO 2 fluxes to the amount of organic carbon in the active layer. Despite substantially higher soil temperatures in 2019, significantly lower CO 2 fluxes per thawed carbon were observed than in 2016, which indicates a decreasing degradability of the remaining organic carbon in the active layer. Also surface erosion, causing a loss of labile carbon, might have contributed to lower CO 2 fluxes in 2019 than in 2016.
Mean CO 2 fluxes normalized to thawed organic carbon were lower at the SF soils than at the TM soils, both under in situ conditions and in laboratory incubations. This is surprising, since the SF soils contained fresh OM from the surface vegetation of the eroding headwall, and it has been shown that fresh OM may increase OM decomposition rates in thawing permafrost (positive priming) (Walz et al., 2017;Wild et al., 2016). The reason for lower decomposition rates in the SF soils might be a higher disturbance due to strong erosion, which was less pronounced at the TM sites. The higher decomposability of the OM at the TM sites is further striking considering its substantially higher age (20-33 kyr BP) than the OM at the SF sites (2-13 kyr BP). These results support data from laboratory incubations that show no effect of organic carbon age on decomposition rates and a particularly high decomposability of organic carbon in the Siberian deposits of the MIS3 Kargin interstadial Walz et al., 2018). Thus, we present evidence that the permafrost OM in the Pleistocene sediments was well preserved during its storage for thousands of years under frozen conditions, and that microbial activity in the permafrost is likely of minor importance at least at the low average temperature of −10°C in the deep permaforst (Boike et al., 2019). The latter conclusion is supported by multi-proxy studies indicating a low decomposition state of the Pleistocene Yedoma deposits in the study region .
Methane emissions from the sampling sites were about two orders of magnitude lower than CO 2 emissions and contributed only a minor amount to total carbon fluxes. Since methane is produced only in the absence of oxygen, its formation in soils is generally restricted to water saturated conditions. Furthermore, methane may be oxidized to CO 2 by aerobic methane oxidizers when migrating through aerobic soils (Knoblauch et al., 2008;Zheng et al., 2018). Therefore, methane emissions depend on a variety of environmental parameters with lowest or even negative fluxes in well drained soils and highest fluxes from water saturated, organic soils vegetated by vascular plants (Christensen et al., 1995;Kwon et al., 2017;Natali et al., 2015;Olefeldt et al., 2013). As expected, mean methane fluxes from the unsaturated, bare soils of the thaw slump were in the lower range of those reported for wet arctic tundra soils in Siberia, Greenland and Alaska (20-250 mg CH 4 -C m −2 d −1 ) (Kutzbach et al., 2004;Kwon et al., 2017;Ström et al., 2015;Vaughn et al., 2016) but, at least in 2016, higher than those from other drained arctic tundra soils (Davidson et al., 2016;Jørgensen et al., 2015;Kutzbach et al., 2004;Kwon et al., 2017). Furthermore, the large variability in methane fluxes with maximum values of up to 132 mg CH 4 -C m −2 d −1 demonstrate substantial heterogeneity of soil conditions but also a high degradability of the thawing permafrost OM even under anaerobic conditions.
The absence of methane emissions at SF1 and SF3 (2019) is striking since these sites had the highest mrcA gene copy numbers and produced methane in laboratory incubations. Hence, the most likely explanation for lacking methane fluxes from these sites is the complete oxidation of produced methane in the unsaturated surface soil (Knoblauch et al., 2008;Zheng et al., 2018). In contrast, the abundance of methanogens in TM1 and TM2 was lowest among all soils and in part even below the detection limit and no methanogenesis was detected in TM2 soil incubations. It has been shown that permafrost deposits, including those from the studied Kurungnakh sites, may contain low numbers of methanogens (Holm et al., 2020;Knoblauch et al., 2018), which may result in the observed very low or absent methane production in recently thawed permafrost (Chowdhury et al., 2015;Knoblauch et al., 2013;Waldrop et al., 2010;Walz et al., 2018). Not only methanogenesis but also other microbial functions, like nitrification, may got lost in Yedoma deposits over the millennia under frozen conditions, which constrain carbon and nitrogen cycling in these soils after thaw (Monteux et al., 2020). In the current thaw slump, the dispersal of methanogens from the active layer 10.1029/2021JG006543 13 of 18 above the headwall, which contained high numbers of methanogens, into the soils of the thaw slump seems essential for active methane production and methane emission, at least in the initial phase after thaw.
The aerobic CO 2 production rates in laboratory incubations of Kurungnakh soils are in the middle of the range of rates reported for mineral permafrost material of the Eurasian Arctic, quantified by short-term incubation studies (≤6 months) at temperatures below 10°C (Faucherre et al., 2018;Gentsch et al., 2018;Zimov et al., 2006). Anaerobic CO 2 production rates are rarely reported from low temperature (<10°C) incubations of thawed permafrost and ranged in Kurungnakh and Samoylov deposits between 0.02 and 0.20 μmol CO 2 g −1 d −1 Walz et al., 2017), which covers the range of anaerobic CO 2 production rates in the SF soils studied here.
It is still unclear to which extent laboratory incubation data represent in situ carbon decomposition rates in thawing permafrost. The current aerobic short-term incubations substantially overestimated in situ CO 2 fluxes. However, even well drained upland soils may contain anaerobic microsites enabling anaerobic microbial processes such as fermentation, denitrification and methane production (Sexstone et al., 1988;von Fischer & Hedin, 2002). The distribution of anaerobic microsites is strongly regulated by water filled pore space and microbial respiration in the soil (Keiluweit et al., 2016;Wagner, 2017). Field methane fluxes indicated anaerobic conditions at least in some parts of the soils. Since a fraction of CO 2 emissions originate from anaerobic OM decomposition, which is about three times slower than aerobic respiration Schädel et al., 2016), aerobic CO 2 production rates will overestimate in situ decomposition rates but may be used as the upper benchmark for potential permafrost OM decomposition. On the other hand, anaerobic decomposition rates from laboratory incubations should underestimate in situ CO 2 fluxes. However, this was only the case for two sites. The reason is likely the increased availability of labile OM after mixing the soil during the preparation of the incubation experiments. This known bias may be avoided by longer incubation studies whose data may then be used for constraining dynamic carbon decomposition models (Elberling et al., 2013;Knoblauch et al., 2013;Schädel et al., 2014). Such models may provide carbon decomposition rates similar to those under in situ conditions, which was demonstrated by using a two pool decomposition model (ICBM), calibrated with data from a three years incubation study with permafrost samples from the study site ( Figure 6).
Surprisingly, in situ methane emissions were generally higher than potential methane emissions. The opposite was expected, since a fraction of the produced methane will be oxidized before emitted under in situ conditions but not in anaerobic incubations. Initially low methane production rates in incubations indicate that methane production is not accelerated by labile OM, as was observed for CO 2 production Lee et al., 2012;Waldrop et al., 2010), but that the conditions for methanogenesis were not optimal at the beginning of the experiments (Chowdhury et al., 2015;Høj et al., 2007;Knoblauch et al., 2013;van Hulzen et al., 1999). Methanogenesis is a syntrophic process requiring a close cooperation with fermenting organisms under specific micro-environmental conditions (Schink, 1997). Disturbances during sample preparation, e.g., by mixing of the soil samples, may disrupt such consortia, change micro-environmental conditions and hence reduce methane production until stable conditions are re-established. Therefore, methane fluxes may generally be underestimated in short-term incubations.

Total Annual Carbon Release
Soil temperature and the amount of thawed organic carbon were the main driving factors of heterotrophic respiration at the study sites. Therefore, the two independent models to extrapolate summer field observations to annual CO 2 fluxes were based on these two factors. The range of soil temperatures was large enough to be used for the Q 10 model calibration only at three sites and we could only use positive temperatures. However, our Q 10 model gives similar results below 0°C as an exponential model fitted to a wide range of winter CO 2 fluxes (Natali et al., 2019). These results give confidence that our Q 10 model can be used for simulating CO 2 fluxes at soil temperatures below those used for calibration. Furthermore, the obtained Q 10 values (2.6-4.1) are in the typical range of Q 10 values reported for the temperature response of CO 2 fluxes from permafrost ecosystems (Fouché et al., 2017;Natali et al., 2019;Wang et al., 2020). The simulation of active layer soil temperatures in 2016 with the JSBACH land surface model was very close to the observations at the nearby tundra on Samoylov when including an insulating moss cover into the model ( Figure 5). Since no vegetation was present at our thaw slump sites, we used the simulated temperature data without a vegetation cover, which were very close to the observations during the field flux measurements. There was only a positive bias of JSBACH concerning the length of the zero-curtain period in autumn, which might result in higher CO 2 emissions simulated with the model than in the field.
Data on annual CO 2 emission fluxes from R h of tundra soils are currently not available and only few studies report R h fluxes during the thaw season, which strongly differ between about 32 and 240 g CO 2 -C m −2 (Eckhardt et al., 2019;Nobrega & Grogan, 2008;Schuur et al., 2009). Heterotrophic respiration strongly depend on ecosystem conditions with highest CO 2 fluxes reported form moist organic tundra soils with a relatively high mean annual temperature (MAT) of −1°C (Schuur et al., 2009) and lowest at wet, mineral tundra soils with low MAT (≤−9°C) (Eckhardt et al., 2019;Nobrega & Grogan, 2008).
Ecosystem respiration in the arctic tundra continues during the non-growing season and may contribute a considerable fraction to the annual CO 2 and methane budget (Euskirchen et al., 2012;Fahnestock et al., 1998;Zimov et al., 1996;Zona et al., 2016). The two models applied in the current study behave differently at subzero temperatures. The Q 10 model simulates CO 2 fluxes at any temperature below 0°C, since in situ CO 2 fluxes were observed even at very low winter temperatures below −15°C (Natali et al., 2019). In contrast, the ICBM is simulating microbial CO 2 production only at temperatures above −10°C. Between about −2°C and −5°C the water availability in the soil sharply decreases (Romanovsky & Osterkamp, 2000), resulting in a drop in microbial activity (Tilston et al., 2010), which becomes practically negligible below −10°C (e.g., Mikan et al., 2002;Natali et al., 2019;Schaefer & Jafarov, 2016). CO 2 fluxes may continue even below −10°C without microbial CO 2 production, e.g., by CO 2 diffusion along the concentration gradient between the soil and the atmosphere. The simulated CO 2 fluxes at sub-zero temperatures are substantial and accounted for 25%-31% of total annual fluxes (160-184 g CO 2 -C m −2 ) from heterotrophic respiration.
Observations of R h fluxes in tundra soils during the winter season are not available but ecosystem respiration was shown to contribute a variable fraction of 15%-55% to annual CO 2 fluxes in Alaskan tundra soils (Euskirchen et al., 2012;Schuur et al., 2009). Since soil temperatures at the Siberian study sites are below −10°C for five to six months of the year (Boike et al., 2019), the minor contribution of winter fluxes to total annual fluxes becomes plausible.
One of the central research questions concerning the permafrost carbon climate feedback is how fast permafrost OM will be mineralized under in situ conditions to CO 2 and methane after thaw. However, estimates on permafrost OM decomposability based on in situ CO 2 fluxes are still not available. Based on our calibrated ICBM, about 0.8% of thawed organic carbon is mineralized to CO 2 in one year. This carbon mineralization rate is much higher than the simulated loss of 15% thawed permafrost carbon until 2100, which was derived from laboratory incubations and dynamic models Schuur et al., 2015). However, our field observations show a decrease of OM decomposition rates over time, which was also observed in laboratory incubations (Lee et al., 2012;Schädel et al., 2014), and is most likely caused by a decrease in the fast decomposable organic carbon pool. Hence, the initial decomposition of almost 1% of thawed permafrost organic carbon during one year will likely further slow down.
Recent upscaling of CO 2 flux measurements indicate that the circum-arctic permafrost region in general and the upland tundra in particular are small annual carbon sources although uncertainties are still high, particularly due to missing data during the long winter period (Belshe et al., 2013;Virkkala et al., 2021). Also multi annual CO 2 flux measurements on Samoylov show that the annual NEE of the polygonal tundra close to Kurungnakh fluctuates over the years (n = 8) between a carbon source and a carbon sink (+26 to −24 g CO 2 -C m −2 ) and it was a slight carbon source in 2016 (+7 g CO 2 -C m −2 ) (data from Holl et al., 2019). Hence, abrupt permafrost thaw, which erodes the vegetation surface and liberates old permafrost OM for microbial decomposition, has a substantial effect on the annual carbon balance turning the tundra from carbon neutral to a CO 2 source. However, the models used in the current study were calibrated only with CO 2 flux measurements during the summer thaw period and with decomposition rates from permafrost deposits of the Lena Delta. Hence, while they improve our process understanding on in situ CO 2 emissions from thawing permafrost substantially, the simulated CO 2 fluxes are more accurate for the summer thaw period and most representative for permafrost deposits of the Lena Delta. For constraining the carbon mass balance of permafrost thawing in the northern hemisphere, pan-arctic field studies are required including the long winter period.

Conclusions
Abrupt permafrost thaw of ice-rich Pleistocene Yedoma deposits liberates OM that is initially rapidly decomposed into mainly CO 2 , turning the carbon-neutral tundra into a substantial greenhouse gas source but the depletion of the small labile carbon pool causes a decrease in CO 2 fluxes over time. The formation of methane from the Pleistocene OM not only requires anaerobic conditions but also the introduction of active methanogenic communities since these are not abundant in the studied Yedoma deposits. Short-term laboratory incubations may provide an upper benchmark for in situ CO 2 fluxes but give only very limited information on in situ methane fluxes. In contrast, dynamic carbon decomposition models, calibrated with long-term incubation data and CO 2 field fluxes provide meaningful estimates of annual in situ CO 2 release of which 25%-31% took place outside the summer thaw season. The bare soils are carbon sources initially after abrupt permafrost thaw, but the development of a vegetation cover will provide an additional carbon input. Therefore, multi-annual observations are required that include carbon fluxes from the upcoming vegetation and lateral fluxes of particulate and dissolved organic carbon to evaluate if these soils will stay a net carbon source, as observed during the first years after thaw, or develop to a net carbon sink in future.