Thawing Siberian permafrost stabilizes organic carbon from recent plant litter inputs

Greenhouse gas release due to microbial decomposition of thawing permafrost organic matter receives ample attention but the other side of the permafrost soil carbon budget, the stabilization of organic matter due to rising plant litter input in a greening Arctic has hardly been addressed. Here we explore whether thawing permafrost material may act as a long-term sink of fresh plant litter carbon. To identify the magnitude and drivers of litter carbon stabilization in thawing permafrost material, we incubated samples from the permafrost layer under oxic and anoxic conditions with ¹³C-labelled plant litter. Subsequently, we used the microbial CO₂ and CH₄ production from the added litter carbon (litter-C) and from the carbon in the thawed permafrost material (permafrost-C) to calibrate a carbon decomposition model with a fast and a slow carbon pool. Beside the size of the different pools, their mean residence times (MRT) were calculated as an indicator for carbon stabilization in these soils. Finally, we fractionated the remaining organic matter into a dissolved, a mineral-associated and a particulate fraction. At the end of the experiment, after nine years, on average 40 % to 60 % of the added litter-C persisted in the thawed permafrost material. The MRT of the slow litter-C pool of 18 years (oxic) and 52 years (anoxic) indicate a substantial stabilization of fresh litter-C over the course of the experiment. More than 80 % of the remaining litter-C was part of the mineral-associated fraction, but in contrast to current understanding, litter decomposability was positively correlated with the size of the mineral bound litter-C pool. Although the fraction of mineral-bound permafrost-C (64 % to 68 %) was significantly smaller than of litter-C, the MRT of the slow permafrost-C pools was more than 10-fold higher. Hence, the size of the mineral bound carbon pool alone may not be a suitable measure of carbon stabilization. We furthermore identified interactions between new litter carbon and pre-existing mineral-bound carbon from the thawed permafrost material as an important driver of litter-C stabilization. Such interactions could reduce net carbon emissions from thawing permafrost and add complexity to the permafrost carbon climate feedback.

1 Introduction

About 15 %–18 % of the land area of the northern hemisphere is underlain by permafrost (Zhang et al., 2000; Obu et al., 2019). Soils in these circum-Arctic terrestrial landscapes store about 1300 Pg organic carbon, 800 Pg of which is preserved in the permafrost (Hugelius et al., 2014). The storage of such large amounts of organic matter (OM) is due to an imbalance between plant net primary productivity and microbial plant litter mineralization. The unique environmental conditions in the permafrost region facilitate the accumulation of OM, as low temperatures slow down OM decomposition. Furthermore, water-saturated, anoxic soils are widespread in the permafrost region since permafrost impedes water drainage. In these cold, waterlogged soils, microbial carbon decomposition is particularly low and OM accumulation is high (Heikkinen et al., 2002; Eckhardt et al., 2019).

The carbon dynamics in permafrost-affected soils differ fundamentally between the active layer, i.e., the surface soil layer that thaws every year, and the permafrost, i.e., the earth material that remains frozen for at least two consecutive years (van Everdingen, 1976), underneath. If vegetated, the active layer receives regular inputs of fresh OM from root exudates and decaying plant material, which are decomposed by microorganisms. Soil organic matter (SOM) is transported by cryoturbation and advection (Beer et al., 2022) into deeper soil layers where decomposition slows down (Kaiser et al., 2007) and SOM is ultimately stabilized through incorporation into the permafrost (Ping et al., 2015). Although microbial activity may proceed at a low level even at sub-zero temperatures (Mikan et al., 2002; Natali et al., 2019), freezing conditions in the permafrost largely protects SOM from decomposition, causing the persistence of labile OM for thousands of years without contact to fresh plant litter.

The strong warming of the Arctic, which is up to four times higher than the global average (Rantanen et al., 2022), causes the active layer to gradually deepen, but also induces abrupt thaw processes, particularly in ice-rich permafrost deposits (Turetsky et al., 2020). Such abrupt thaw processes induce fundamental changes in the land surface through thermo-erosion and exposes material from the deep permafrost layer to the surface (Lewkowicz and Way, 2019). Thawing permafrost liberates OM, which is decomposed by microorganisms to the greenhouse gases carbon dioxide (CO₂) and methane (CH₄) (Schädel et al., 2016; Knoblauch et al., 2018). But importantly, Arctic warming also causes longer growing seasons (Euskirchen et al., 2006; Arndt et al., 2019) with higher plant productivity (Natali et al., 2012; Berner et al., 2020) and increased inputs of plant litter (Elmendorf et al., 2012) as source for SOM formation.

The permafrost region has been a net carbon sink in the past, documented by the huge accumulation of organic carbon. Whether it will remain a carbon sink or turn into a carbon source depends on how primary productivity and OM decomposition respond to climate change. Current studies based on field flux measurements and modelling approaches indicate that the permafrost region may still be a weak carbon sink, at least in the boreal forest, but the tundra regions seem close to the threshold or have already switched to a net carbon source (Belshe et al., 2013; Virkkala et al., 2021; Hugelius et al., 2024). However, uncertainties remain high, also due to an incomplete understanding of long-term processes.

Although the decomposition of OM in thawing permafrost and the subsequent release of carbon has received ample attention (Elberling et al., 2013; Schädel et al., 2014; Drake et al., 2015; Faucherre et al., 2018; Gentsch et al., 2018; Beer et al., 2022; Guo et al., 2024), little is known about the stabilization of fresh plant litter in thawing permafrost, despite its potential impact on the long-term carbon budget. Thawing removes the freeze-induced stabilization of OM in the former permafrost, increases the volume of the active layer, and allows fresh plant litter to be incorporated into SOM. Microbial SOM decomposition not only produces CO₂ and CH₄ but also promotes the stabilization of litter decay products and microbial necromass through different mechanisms such as physical protection or sorption onto mineral surfaces (Schmidt et al., 2011; Cotrufo et al., 2013; Begill et al., 2023; Qin et al., 2024). These processes and the factors determining the long-term persistence of OM from recent plant litter in thawing permafrost are poorly understood.

Soil organic matter is a mixture of organic material of different age, molecular composition and chemical properties. To simulate future carbon release from soils, organic matter decomposition models need to differentiate among several pools with different decomposability. Such process-based models are calibrated using experimentally determined OM decomposition rates (Tuomi et al., 2009; Schädel et al., 2014; Knoblauch et al., 2018). However, these models do not allow characterizing the composition of SOM or identifying the mechanisms of SOM stabilization. Therefore, physical SOM fractionation methods followed by the chemical analysis of the derived fractions are applied (Schrumpf et al., 2013; Torn et al., 2013; Prater et al., 2020; Martens et al., 2023). It has been shown that the OM associated with the mineral phase (MAOM) is generally older and has a higher $\mathrm{C}/\mathrm{N}$ ratio than the particulate OM fraction (POM), indicating a higher decomposability or lability of the POM fraction (Torn et al., 2013; Prater et al., 2020; Martens et al., 2023). However, only a few studies have measured the microbial decomposability of the different OM fractions in permafrost soils (Jagadamma et al., 2014; Gentsch et al., 2015a) and no clear difference in OM lability between these fractions was detected. On the other hand, carbon from plant litter may rapidly be incorporated into the mineral fraction (Swanston et al., 2005; Vogel et al., 2014; Haddix et al., 2016) but the lability of this recently bound OM remains largely unknown.

Our study, based on a nine years lasting incubation experiment with ¹³C-labelled plant litter, aims to improve our understanding on two contrasting processes and their interaction after permafrost thaw: (1) the long-term stabilization of OM from decaying fresh plant litter in thawing permafrost and (2) the mineralization of old OM in the thawing permafrost. We hypothesize that (h1) fresh plant litter is stabilized on mineral surfaces thereby contributing to the SOM build-up in thawing permafrost soils and that (h2) the persistence of mineral bound plant litter is lower than the persistence of mineral bound OM in the thawed permafrost.

2 Materials and Methods

2.1 Sampling sites

Samples from the permafrost layer were collected from the uppermost 4.3 m of deposits from two islands in the Lena River Delta, Russia as described in Knoblauch et al. (2013). The island Kurungnakh (72.327503° N, 126.277597° E) is comprised of Pleistocene Yedoma deposits, which are covered by silty, cryoturbated Holocene sediments. Samples were taken at from a core drilled from the surface into the permafrost layer and from a headwall of a retrogressive thaw slump (Knoblauch et al., 2013; Bischoff et al., 2013) In contrast, Samoylov Island (72.369444° N, 126.475° E) is comprised mainly of sandy, fluviatile sediments that were deposited during the late Holocene. There, a core was drilled from the surface to a depth of 4.5 m. The tundra vegetation at both sites is composed of a moss/lichen layer and a grass/sedge layer. The dominant vascular plants in the wet tundra on Samoylov and Kurungnakh are sedge species of the Carex aquatilis aggregate (Kutzbach et al., 2004; Lashchinskiy et al., 2020).

Samples from Kurungnakh are characterized by relatively high organic carbon concentrations (33.9–121 mg C g⁻¹), intermediate $\mathrm{C}/\mathrm{N}$ ratios between 12 and 19 and increasing pH values ranging from extremely acidic at the surface of the permafrost layer to neutral in the deepest layers (Knoblauch et al., 2013), (Appendix, Table A1). The texture of the Kurungnakh samples is dominated by the silt fraction with generally more than 50 % and a relatively constant clay fraction slightly above 20 %. The organic carbon concentrations in the samples from Samoylov (6–47 mg C g⁻¹) are lower than on Kurungnakh while the $\mathrm{C}/\mathrm{N}$ values are higher (15–25). The pH values are relatively stable and in the slightly acidic range (Knoblauch et al., 2013), (Appendix, Table A1). The sand content of the Samoylov permafrost samples is decreasing from 85 % at the surface of the permafrost layer to 37 % in the deepest layer while the clay content ranges from 3 % to 13 %.

2.2 Incubation experiment

The incubation experiment has been described in detail previously (Knoblauch et al., 2018). Briefly, twelve samples from the permafrost layer in Kurungnakh and Samoylov (Appendix, Table A1) were incubated in triplicate at 4 °C under oxic and anoxic conditions in 120 mL injection vials that were closed with butyl rubber stoppers and were not opened during the experiment to keep moisture constant (Appendix, Fig. A1). Anoxic conditions were established by repeatedly evacuating and flushing the headspace of the bottles with molecular nitrogen at the beginning of the experiment. The headspace of the oxic incubations was repeatedly flushed with synthetic air (20 % O₂, 80 % N₂) when CO₂ concentrations exceeded 3 %. During the first four years, labile OM was decomposed and an active methanogenic community developed in the anoxic incubations (Knoblauch et al., 2013). Data from this pre-incubation experiment are not considered in this manuscript. After four years, the amount of organic carbon that was mineralized to CO₂ and CH₄ in each of the samples was replenished with organic carbon from ¹³C-labelled Carex aquatilis (Appendix, Table A1) that was grown on Samoylov Island. Such prepared incubations of permafrost material contained two isotopically distinct carbon sources, the carbon from the thawed permafrost samples (permafrost-C, −26.5 ‰ to −29.1 ‰ VPDB) and the carbon from the added Carex litter (litter-C, +774 ‰ VPDB). The samples were then incubated for a further nine years. Based on the different δ¹³C-signatures of permafrost-C and litter-C, a two-endmember model was applied to quantify the amount of CO₂ and CH₄ produced from the two carbon sources permafrost-C and litter-C (Amundson and Baisden, 2000). Organic matter decomposition rates were calculated from three successional measurements at the beginning (initial rates) and end (final rates) of the incubations, and were related to the concentration of the respective carbon source present at the respective time interval.

2.3 Gas quantification and carbon stable isotope analysis

The production of CO₂ and CH₄ was quantified by measuring the gas pressure (LEO1, Keller Druckmesstechnik, Switzerland), and gas concentrations in the headspace of the incubation flasks with a gas chromatograph (GC 7890, Agilent Technologies, USA) as described before (Knoblauch et al., 2018).

The δ¹³C-values of CO₂ and CH₄ were determined with an isotope ratio mass spectrometer (Delta V, Thermo Fischer Scientific, Bremen, Germany) equipped with a GC Isolink (Thermo Fischer Scientific, Milano, Italy). CH₄ was measured against the external standards NIST 8561 (−73.27 ‰ VPDB) and USGS HCG-1 (−1.51 ‰ VPDB) while CO₂ was measured against LSVEC (−46.60 ‰ VPDB), CO-8 (−5.76 ‰ VPDB) and IAEA 303A (+93.3 ‰ VPDB). Stable carbon isotope analysis of particulate organic carbon and DOC were carried out with a Delta VPlus isotope ratio mass spectrometer and a ConFlo IV (ThermoFisher Scientific, Bremen, Germany) that was connected to a Flash EA (ThermoFisher Scientific, Bremen, Germany) for δ¹³C-analysis of carbon in the bulk samples, the POM and the MAOM fraction, and with an Aurora 1030 (OI Analytical, Yellow Springs, USA) for δ¹³C-DOC analyses. ¹³C-analysis of DOC and solid samples were calibrated against the isotope standards IAEA-CH-7 (−10.45 ‰ VPDB), USGS40 (−26.39 ‰ VPDB) and USGS41 (+37.63 ‰ VPDB).

2.4 Soil analysis

Total carbon and nitrogen in the bulk samples and in the POM and MAOM fractions were determined with an elemental analyser (VarioMAX Elementar Analysensysteme GmbH, Hanau, Germany). If the available sample amount was too small, carbon analysis was conducted with the Flash EA in parallel to δ¹³C-analysis. Carbon and nitrogen concentrations were calibrated against external standards with carbon concentrations between 48.3 % and 2.2 % and nitrogen concentrations between 2.30 % and 0.20 % (IVA Analysentechnik, Meerbusch, Germany). Soil texture was determined by sieving and sedimentation (ISO 11277:2020, 2020). Therefore, air dried soil (<2 mm) was suspended in a 0.1 M Na₄P₂O₇ solution and sieved through a 63 µm mesh sieve. The retained fraction was dried at 105 °C and subsequently fractionated by sieving through 630, 200, 125, and 63 µm mesh sieves. The fraction of <63 µm was separated by sedimentation using a Sedimat 4-12 (UGT, Müncheberg, Germany). Organic matter was destroyed with a 30 % H₂O₂ solution prior to texture analysis.

2.5 Organic matter fractionation

At the end of the incubation, the dissolved organic matter (DOM) was extracted from each of the replicates by mixing air dried soil with 0.01 M CaCl₂ solution (1:2; w:v), overhead shaking for 15 min and centrifugation for 30 min at 3000 g. The supernatant was filtered through a GFF filter, acidified to pH 2, and stored at 4 °C until analysis. Subsequently, the OM from each replicate of each permafrost sample was separated by density fractionation. Due to the low amount of sample (7–15 g) we separated only two fractions, the POM and the MAOM. Dry permafrost samples were mixed with a Na-polytungstate solution (TC Tungsten Compounds, Grub am Forst, Germany) with a density of 1.6 g cm⁻³ (Cerli et al., 2012) at a ratio of 1:5 (w:v). To break up aggregates we used a sonication with 200 J mL⁻¹ (Sonoplus HD 200/Sonotrode MS 70, Bandelin, Berlin, Germany) in an ice bath. Subsequently, the samples were centrifuged at 2000 g for 30 min and the floating material (POM) was carefully separated from the sedimented material (MAOM) with a vacuum pump. The POM fraction was washed with deionized water by pressure filtration until the electrical conductivity in the washing solution decreased below 50 µS. Subsequently the POM fraction was detached from the filter (0.22 µm PVDF membrane, Berrytec, Harthausen, Germany) by ultrasonication and subsequently dried at 60 °C. The MAOM fraction was washed by centrifugation (6000 g) with deionized water until the electrical conductivity in the washing solution decreased below 50 µS. The sediment with the MAOM fraction was then dried at 60 °C. The median carbon recovery, i.e. the sum of carbon in the DOM, POM and MAOM fraction divided by the total organic carbon in the unfractionated sample was 96.2 % (IQR 10.4). The share of carbon in the three different OM fractions (DOM, POM, MAOM) to the total permafrost-C and litter-C was calculated by dividing the amount of carbon in the respective fraction through the total amount of carbon recovered.

2.6 Organic matter dynamic decomposition model

The cumulative CO₂ and CH₄ production from thawed permafrost-C and litter-C from each sample and replicate were used to calibrate an organic carbon decomposition model as described before (Knoblauch et al., 2013). The model follows the principles of the Introductory Carbon Balance Model (Andrén and Kätterer, 1997). A first-order kinetics equation represents the change of organic carbon content in time. This equation is applied to two carbon pools with a high and a low rate constant, respectively. The initial fraction of the fast pool is treated as a parameter and the initial fraction of the slow organic carbon pool is then calculated as the difference to the total organic carbon content. Using a nonlinear least-squares approach with a trust-region-reflective algorithm in MATLAB R2023b (MathWorks Inc., USA), the following four model parameters were optimized: two decomposition rate constants (fast and slow pool), initial fast carbon pool fraction, and the stabilization coefficient. Mean residence time of each pool (MRT) is defined as the reciprocal of the decomposition rate constant. In case of oxic incubations, the model was fitted with the cumulative CO₂ production, in case of anoxic incubations we used the total carbon mineralization, i.e., the sum of CO₂ and CH₄ production. We had to exclude replicates from the analysis for which the two-pool model could not be fitted to the data because the data was already too linear after the long pre-incubation period. In this case the MRT of the slow pool reached 1 million years, the upper limit of the model. The total number of replicates analysed for the oxic incubations was n=34 (litter-C pools) and n=31 (permafrost-C pools) and for the anoxic incubations n=28 (litter-C pools) and n=22 (permafrost-C pools).

2.7 Statistics

Prior to all analyses, data were inspected visually for their distribution via histograms, QQ-plots and by a Shapiro-Wilk test. Mean values were compared using a two-sided independent t-test. In the case of non-normal distribution, the medians were tested using a Mann-Whitney rank sum test or a Kruskal-Wallis ANOVA. Significance was tested on a level of p<0.05, using the Bonferroni correction in case of multiple tests. If mean values are presented, they are followed by the standard deviation (±SD), in case of medians they are followed by the interquartile range (IQR). Since most of the data did not show a normal distribution, we used for the correlation analysis a Spearman rank order correlation. Due to the correlation between the sand, silt, and clay fractions, the clay fraction was included as the single variable to represent soil texture. These statistical analyses were conducted with IBM SPSS Statistics 29.0 (IBM Corporation, Armonk, USA).All further statistical analyses were performed in R version 4.3.2. To identify the most important controls on organic matter stabilization we applied a random forest (RF) machine learning approach. Input variables for the RF analyses were selected based on stepwise linear regression using R-package MuMIn (Bartón, 2023). Values were standardized using z-score scaling [(x − mean)$/$standard deviation]. The best linear model was selected based on the lowest Akaike Information Criterion, a Delta value <2 and ≥6° of freedom to obtain a sufficient number of variables for the subsequent RF analysis. Stepwise linear regression analysis was conducted for eight separate models using the following dependent variables: Fraction of litter-C decomposed, fraction of permafrost-C decomposed, MRT stable litter-C pool, MRT stable permafrost-C pool. The most important predictors for these dependent variables were identified separately for oxic and anoxic samples. Results of the stepwise-linear regression are listed in Appendix, Table A2. For the subsequent RF analysis, models for oxic and anoxic samples were merged to include the same set of predictor variables. The clay fraction was included as the single variable to represent soil texture. RF models were constructed for the same models created in the stepwise linear regression using the R-package randomForest (Breiman, 2001). We used 500 trees (ntree = 500) to create the random forests, and the function `tuneRF' was applied to determine the number of variables tried at each split individually for each of the eight models. Highly-correlated variables (Spearman's correlation coefficient > 0.9) were excluded from the RF models. Variable importance was assessed using the average increase in node purity of the regression trees based on splitting on the various environmental variables (Louppe et al., 2013). Due to the lower number of replicates from the anoxic incubations that were used for the two-pool model analysis we excluded the MRT of the anoxic incubations from the RF analysis. The predictive power in the RF analysis was high (64 %–92 % explained variance). Due to the relatively small number of observations in our dataset, we used partial dependence plots (R-package pdp) to assess the validity of the assumed responses (Appendix Fig. A2).

3 Results

3.1 Mineralization of litter carbon versus permafrost carbon

At the beginning of the incubation experiment, litter-derived CO₂ and CH₄ dominated gaseous carbon production. Initial litter-C mineralization rates were 35 (IQR 41) and 165 (168) times higher than initial permafrost-C mineralization rates under oxic and anoxic conditions, respectively, with about 60 % of total litter-C mineralization taking place in the first year (Appendix, Fig. A1a, Table A3). After one year, litter-C mineralization rates decreased significantly more rapidly than those of permafrost-C (p<0.001). By year nine, median litter-C mineralization rates dropped to 9.3 (12.6) % (oxic) and 1.5 (3.0) % (anoxic) of their initial values. However, the median OM mineralization rates of litter-C at the end of the experiment were still about ten times higher than those of the remaining permafrost-C under both oxic and anoxic conditions (Appendix, Table A3).

The fraction of total litter-C that was decomposed during the incubation period ranged from 23 % to 97 % (mean 61 ± 21 %) under oxic conditions and was significantly higher (p<0.001) than under anoxic conditions (mean 41 ± 12 %, Fig. 1a). The fraction of permafrost-C that was decomposed was significantly lower than that of litter-C, ranging between 1.3 % and 16.3 % (median 4.8 %) under oxic conditions and between 0.5 % and 4.4 % (median 1.0 %) under anoxic conditions (Fig. 1b).

Figure 1 Total amount of organic carbon decomposed into CO₂ and CH₄ from litter-C (a) and permafrost-C (b) over nine years in samples incubated under oxic and anoxic conditions. Values are relative to the total amount of litter-C and permafrost-C at the beginning of the experiment. Error bars represent one standard deviation of the mean (n=3).

Under oxic conditions, the median size of the fast litter-C pool, as determined by the two-pool carbon decomposition model, was slightly higher than of the slow litter-C pool, while the opposite was the case under anoxic conditions (Fig. 2a, c). In contrast, permafrost-C pools were clearly dominated by the slow pool with 92.7 (12.9) % under oxic and 98.8 (2.2) % under anoxic conditions (Fig. 2c).

Figure 2 Size of the fast (a) and slow (c) litter-C and permafrost-C pools and mean residence times (MRT) of the respective carbon pools (b, d) determined by a two-pool carbon decomposition model. Significant differences are indicated by * = p<0.05, ** = p<0.01 and *** = p<0.001 (Kruskal-Wallis ANOVA). The Box shows the $\mathrm{25}/\mathrm{75}$ percentiles, the whiskers the $\mathrm{5}/\mathrm{95}$ percentiles, the horizontal line the median and the closed square the mean. Outliers are presented as open circles.

The MRT of the fast carbon pools of litter-C and permafrost-C were similar, with median values below 0.5 years, except for the fast permafrost-C pool under anoxic conditions with a significantly higher median MRT (Kruskal-Wallis ANOVA, p=0.004; Fig. 2b). However, the median MRT of the slow permafrost-C pool was more than ten times higher than that of the slow litter-C pool, both under oxic conditions (195 year vs. 17.6 year) and anoxic conditions (1572 year vs 52.4 year) (Fig. 2d).

The median ratio of initial oxic to anoxic litter mineralization rates was 3.6 (7.2) but this increased substantially to 8.0 (4.7) towards the end of the experiment. A similar increase was observed for the oxic/anoxic ratios for permafrost-C decomposition, rising from an initial ratio of 3.3 (1.3) to a final ratio of 6.9 (4.9).

3.2 Fractionation of permafrost carbon and litter carbon

Most organic carbon from both carbon sources remaining at the end of the incubations was found in the mineral-associated (MA) fraction (Fig. 3a, b). The median $\mathrm{C}/\mathrm{N}$ ratio of 15.7 in the total MAOM fraction (sum of MA-permafrost-C and MA-litter-C) was significantly lower than the ratio of 29.5 in the total POM fraction.

The median values for mineral-associated permafrost-C under oxic and anoxic conditions were 64.4 (20.0) % and 68.0 (17.0) %, respectively (Fig. 3b). Surprisingly, we found a significantly higher proportion of litter-C remaining after nine incubation years in the mineral-associated fraction: 82.5 (35.3) % under oxic conditions and 83.8 (21.4) % under anoxic conditions (Fig. 3a). The median contribution of DOC to total permafrost-C and litter-C was generally below 1 %, except for the dissolved litter-C under anoxic conditions, which contributed 4.6 % to the total litter-C (Appendix Fig. A3).

Figure 3 Contribution of litter-C (a) and permafrost-C (b) to the particulate (P-carbon), the mineral associated (MA-carbon) and the dissolved (D-carbon) carbon fractions in permafrost samples from Kurungnakh and Samoylov. Error bars representing one standard deviation of the mean (n=3) are presented to the right of the respective bar.

3.3 Factors affecting organic carbon stabilization

The final permafrost-C decomposition rate, which was calculated based on the amount of carbon remaining at the end of the incubations, was the most important predictor of the parameters used to describe permafrost-C stabilization, i.e. the fraction of permafrost-C mineralized and the MRT of the slow pool of permafrost-C (Fig. 4). Texture, represented by the size of the clay fraction, was also an important predictor of not only permafrost-C and total nitrogen concentrations, but also of the fraction of permafrost-C decomposed and the MRT of the slow permafrost-C pool. Furthermore, the size of the clay fraction was negatively correlated with the final permafrost-C decomposition rates. However, the size of the MA-permafrost-C pool did not correlate with any of the parameters related to the decomposition of permafrost-C such as permafrost-C decomposition rates, the fraction of permafrost-C decomposed or the MRT of the stable permafrost-C pool (Fig. 4a). In the RF model (Fig. 4d, e), the size of the mineral-associated permafrost-C was also one of the least important parameters explaining permafrost-C stabilization.

Figure 4 (a) Correlation table showing Spearman's correlation coefficients between soil properties and the fraction of permafrost-C and litter-C decomposed and the MRT of the slow pool of litter-C and permafrost-C. Results from samples incubated under oxic conditions are presented in the top right-hand corner, and results from anoxic incubations are presented in the bottom left-hand corner. The colours and size of the circles code the size of the correlation coefficient, the asterisk indicates the probability values with * = p<0.05, ** = p<0.01, *** = p<0.001. Predictors: (b) relative importance of predictors on the fraction of litter-C decomposed, (c) the mean residence time (MRT) of slow litter-C pool, (d) the fraction of permafrost-C decomposed and (e) the MRT of the slow permafrost-C pool . The relative importance was determined with a random forest model for the individual treatments. The MRT of anoxic incubation was excluded from the RF analysis since the two-pool model could not be fitted to the data of several anoxic replicates. For the direction of responses see Appendix Fig. A2a–d. Numbers in panels (b)–(e) give the percentage of variance in the data explained by the model.

The final litter-C decomposition rates were the most important predictor of both the fraction of litter-C decomposed and the MRT of the slow litter-C pool (Fig 4a, b, c). But surprisingly, and contrary to permafrost-C decomposition, the fraction of litter decomposed correlated positively with the proportion of MA-litter-C, which was also the parameter in the RF-model explaining most of the variance of the data under oxic conditions. Additionally, the proportion of MA-litter-C correlated negatively with the MRT of the slow litter-C pool; that is to say, the more litter-C that was bound to the mineral fraction, the lower the MRT of the slow C-pool. Texture was only a minor predictor of litter stabilization, since the size of the clay fraction only correlated with the parameters of litter stabilization under anoxic conditions, and only under anoxic conditions did the size of the clay fraction explain a substantial fraction of the data variance in the RF-model (Fig. 4 b). Interestingly, the proportion of MA-litter-C showed a strong positive correlation with the proportion of MA-permafrost-C under both oxic and anoxic conditions.

4 Discussion

4.1 Organic matter decomposability

The high initial carbon mineralization rates of the added plant litter demonstrate the high share of labile OM in the Carex aquatilis litter that may readily be mineralized, as well as a microbial community adapted to the substrate and the low incubation temperature of 4 °C. Litter-C loss during the nine-year incubation period is similar to the carbon loss from different sedges (44.7 %–61.0 %) incubated for two years under in situ conditions in Arctic soils (Cornelissen et al., 2007). However, the 90 % decrease of litter-C decomposition rates towards the end of the incubation indicates that the remaining 40 %–60 % of the added litter-C may resist a rapid decomposition.

In contrast, the carbon loss from permafrost-C of our Siberian permafrost samples (0.5 %–16 %) is at the lower end of the range (9 %–75 %) reported for a similar 12-year oxic incubation experiment with three permafrost soil samples from Greenland (Elberling et al., 2013). Prior to adding litter-C to our permafrost samples, they were pre-incubated for four years until the most labile permafrost-C had been mineralized (Appendix, Fig. A1) (Knoblauch et al., 2013). Including the carbon loss during this pre-incubation period increases the permafrost-C loss over the 13-year incubation period to a median of 11 % (oxic) and 2.3 % (anoxic) in our study. Decade-long incubation experiments with permafrost samples are extremely rare but simulations of oxic carbon decomposition over ten incubation years based on data from over 120 samples from the circum-Arctic permafrost region indicate a wide range of carbon loss between 0.7 % and 75 % with a median value of 7.5 % of the initial soil organic carbon content (Schädel et al., 2014). Therefore, our permafrost samples seem to represent the median permafrost OM decomposability rather than the reported extremes.

The about tenfold difference between final litter-C and permafrost-C mineralization rates and between the MRT of slow carbon pools indicate different properties of OM from permafrost and of the remaining plant litter. The MRT of the slow litter-C is in the order of decades, consistent with the carbon turnover time of tundra ecosystems (Carvalhais et al., 2014). In contrast, the MRTs of centuries to millennia of the slow permafrost-C pool, which represents more than 90 % of total permafrost-C, disagree with reports of high amounts of labile, rapidly decomposable OM in permafrost deposits, which were identified using the chemical composition of permafrost OM (Waldrop et al., 2010; Mueller et al., 2015), differences in carbon contents between frozen and unfrozen sediments (Vonk et al., 2012) or incubation experiments over relatively short time periods (Waldrop et al., 2010). Initial high OM decomposition rates in thawing permafrost may be due to the elevated amount of DOC in permafrost, which is rapidly decomposed upon thaw (Vonk et al., 2013; Drake et al., 2015). However, data from oxic long-term (≥1 year) incubation studies from more than 20 permafrost affected circum-Arctic study sites show a significantly higher OM decomposability in the uppermost 1 m surface soil, comprising the active layer, than from the permafrost located below 1 m (Schädel et al., 2014). Similar results were obtained by incubating over 100 samples from permafrost soil profiles (Faucherre et al., 2018). Therefore, current long-term OM decomposition measurements are consistent with our findings of a relatively low long-term decomposability of deep permafrost OM.

In our experiments, the initial ratio of oxic to anoxic carbon-based litter mineralization rates was similar to the ratio of 3.4 reported for a wide range of permafrost soil incubations (Schädel et al., 2016). The substantial increase of this ratio towards the end of experiment proves a stronger decrease of anoxic than oxic OM decomposition rates over time. The lower decomposability of OM under anoxic versus oxic conditions is well established (Bridgham et al., 1998; Schädel et al., 2016) and likely due to the lower energy yield of microbial OM decomposition in the absence of oxygen (LaRowe and Van Cappellen, 2011). However, similar decomposition rates under oxic and anoxic conditions have repeatedly been observed in short-term incubations and during the initial phase of OM decomposition (Bastviken et al., 2001; Lin et al., 2021; Knoblauch et al., 2021) indicating that labile OM may support high anoxic decomposition rates. Such labile OM may originate from fresh plant litter. Furthermore, anoxic conditions followed by microbial iron reduction causes the dissolution of iron(III) (oxoyhydr)oxides and the subsequent release of co-precipitated OM, i.e. labile OM formerly bound to the iron minerals is released (Patzner et al., 2020; Dong et al., 2023). A similar decomposability of labile litter OM under oxic and anoxic conditions is supported by similar MRT of the fast litter-C pool in our incubations in the presence and absence of oxygen. Nevertheless, over longer time scales, decomposition rates under anoxic conditions are lower than those under oxic conditions as observed in our incubation study, causing significantly higher MRT and the highest carbon accumulation in water saturated soils (Köchy et al., 2015).

We selected samples from the permafrost layer on Kurungnakh and Samoylov with highly variable soil properties (pH, organic carbon and nitrogen content, texture, $\mathrm{C}/\mathrm{N}$ ratio, etc.) to represent a wide range of environmental conditions at the time these permafrost soils formed during the Holocene. Since all these soil properties affect OM decomposability and OM stabilisation, we also found a high variability in the observed carbon decomposition data and the results from the two-pool carbon decomposition model. A further reason for the observed variability in OM decomposability may be differences in the composition of clay minerals and the content in iron and aluminium oxides, which all play a crucial role for the interaction between the mineral surfaces in soils and OM stabilization (Lim et al., 2022; Jia et al., 2024; Singh et al., 2018). However, due to insufficient sample material, we were unable to investigate the extent to which the variability in our carbon decomposition and stabilisation data is due to the mineralogy of the clay and iron oxides.

4.2 Properties of persistent organic carbon

The significantly higher $\mathrm{C}/\mathrm{N}$ ratio of the POM in comparison to the MAOM remaining at the end of the incubations indicate different properties of the OM in these two fractions. The $\mathrm{C}/\mathrm{N}$ ratio has been used as an indicator of OM decomposition (Schädel et al., 2014; Kuhry and Vitt, 1996) as $\mathrm{C}/\mathrm{N}$ ratios of SOM decrease during on-going decomposition (Bonanomi et al., 2013; Gentsch et al., 2015b), which would indicate that the MAOM fraction contains more decomposed OM. However, lower $\mathrm{C}/\mathrm{N}$ ratios in the MAOM than the POM fraction may also be caused by different sources of OM of the two fractions (Haddix et al., 2016), with a preferential accumulation of microbial necromass, with generally low $\mathrm{C}/\mathrm{N}$ ratios, in the MAOM fraction (Wang et al., 2020) and a relative accumulation of refractory structural plant components such as lignin with a high $\mathrm{C}/\mathrm{N}$ ratio in the POM fraction (Cheng et al., 2023). As we lack data on the sources of the OM in the different fractions, we cannot substantiate either explanation for the higher $\mathrm{C}/\mathrm{N}$ ratio in the POM fraction.

The proportion of mineral-associated permafrost-C under oxic and anoxic conditions, ranging from 60 % to 70 %, falls within the values previously reported for Siberian permafrost samples (Beer et al., 2022; Martens et al., 2023; Gentsch et al., 2015b) as well as for deeper permafrost (>∼ 80 cm) in drained thermokarst lake basins in Alaska (Mueller et al., 2015). In contrast, the active layer may be dominated by the POM fraction (Prater et al., 2020; Mueller et al., 2015), probably due to repeated input of POM through decaying plant litter. The surprisingly high share of more than 80 % of litter-C remaining after nine incubation years in the mineral-associated fraction indicates an efficient incorporation of plant-litter decomposition products and microbial necromass into the mineral fraction over the course of the experiment (Swanston et al., 2005; Haddix et al., 2016).

Although the DOC pool is substantially smaller than the particulate and the mineral-associated carbon pools in permafrost soils, as shown in this and previous studies (Xu et al., 2009), it plays a central role in the decomposition and stabilization of organic matter since it is highly mobile, fastest decomposable (Vonk et al., 2013; Drake et al., 2015) and both the source for greenhouse gases when decomposed by microbes but also for OM stabilized at mineral surfaces (Cotrufo et al., 2013). The significantly higher concentrations of DOC from litter than from permafrost-C provide further evidence of the relatively higher lability of the remaining litter-C compared to permafrost-C.

4.3 Control of organic carbon stabilization

We use the MRT of the slow carbon pools and the fraction of litter-C or permafrost-C decomposed as indicators of organic carbon stabilization. The fact that the final mineralization rates of litter-C and permafrost-C were the most important predictors of carbon stabilization, while the initial mineralization rates provided little additional information, highlights the importance of long-term incubations when studying OM stabilization. Among the soil properties, it was mainly the size of the clay fraction that affected permafrost-C stabilization. The central role of clay and iron minerals for organic matter stabilization is well established (Kögel-Knabner et al., 2008; Dong et al., 2023), primarily because they provide a large mineral surface area for binding organic matter (Gentsch et al., 2015a). The MAOM is generally considered to be the least decomposable and therefore most stable fraction of organic matter in soils (Kögel-Knabner et al., 2008; Mueller et al., 2015; García-Palacios et al., 2024). This interpretation is mainly based on observations of lower $\mathrm{C}/\mathrm{N}$ ratios, lower ¹⁴C-content and lower abundance of functional groups associated with labile OM in the MAOM fraction (Mueller et al., 2015; Beer et al., 2022; Prater et al., 2020). Although we found a significantly lower $\mathrm{C}/\mathrm{N}$ ratio in the MAOM than in the POM fraction, indicative for a higher grade of degradation (Schädel et al., 2014; Kuhry and Vitt, 1996) and a contribution of microbial necromass to the MAOM fraction (Wang et al., 2020), the size of the MA-permafrost-C pool did not contribute to the explanation of permafrost-C stabilization. Therefore, our long-term decomposition experiment provided no evidence that the MAOM pool in permafrost is the most resistant or that the size of the POM pool can be used as an indicator of the amount of labile, fast-decomposable OM in permafrost. The few existing, direct measurements of the decomposability of MAOM in comparison to POM or bulk OM in permafrost soils showed that MAOM in the surface active layer is significantly less decomposable than POM or bulk OM but that this is not the case for MAOM in permafrost (Gentsch et al., 2015a, 2018). The reason for the lack of difference between permafrost MAOM and POM degradability may lie in the generally lower decomposability of permafrost OM in comparison to OM from the active layer (Schädel et al., 2014). Furthermore, POM in permafrost is preferentially occluded in aggregates (Martens et al., 2023; Beer et al., 2022), which makes the latter POM fraction less accessible to microbial decomposition.

The factors affecting litter-C stabilization differed from those affecting permafrost-C stabilization, with the exception that the final litter decomposition rate was also one of the most important predictors for both the fraction of litter decomposed and the MRT of the stable litter-C pool. Contrary to previous findings showing that the mineral fraction contains the most stable carbon in soils (Kögel-Knabner et al., 2008; Mueller et al., 2015; García-Palacios et al., 2024; Heckman et al., 2022; Schrumpf et al., 2013; Qin et al., 2021) our study revealed a decrease in litter-C stabilization with an increasing proportion of litter-C in the mineral fraction. The more litter-C was associated with the mineral fraction at the end of the experiment the higher the final litter-C decomposition rates. One possible reason for this unexpected result is that we fractionated the litter-C at the end of our incubation experiment, after the most labile OM had already been mineralized leaving the least decomposable OM in the POM fraction behind. In this context is has to be considered that we initially added litter-C as POM. This means that all litter-C that was present in the MAOM fraction at the end of the experiment was bound to the mineral fraction during the nine years lasting incubation.

The strong positive correlation between MA-litter-C and MA-permafrost-C, and the absence of a correlation between the clay fraction and litter stabilization under oxic conditions indicates that mineral bound permafrost-C is more important for incorporating fresh plant litter into the mineral fraction than the size of the clay fraction itself. These results are supported by previous findings showing that fresh plant litter binds to pre-existing organo-mineral clusters rather than to free mineral surfaces (Vogel et al., 2014; Kang et al., 2024). The presence of C-pools with different degradability in the MAOM fraction has previously been proposed, for example based on the rapid incorporation of fresh plant litter into the MAOM fraction (Swanston et al., 2005). Kleber et al. (2007) proposed a model of organo-mineral interactions based on multiple layers, characterized by an inner contact zone between the mineral surface and OM with the highest strength of association, and an outer zone where OM is more loosely bound to pre-existing OM on the mineral surface and can more easily be decomposed. Evidence for a preferential binding of fresh OM to pre-existing OM at mineral surfaces (OM-OM interactions) and a relatively labile mineral associated OM fraction have been repeatedly described (Murphy et al., 1990; Possinger et al., 2020; Kang et al., 2024; Jilling et al., 2025). The nature of these OM-OM interaction is not well understood but hydrophobic associations (Coward et al., 2019) and polysaccharide-rich microbial necromass (Kang et al., 2024) are likely involved. Such OM-OM interactions align with our observations of a significantly lower MRT of the slow litter-C pool compared to the slow permafrost-C pool. Assuming that the decomposition products of added litter OM are bound to the outer layer of pre-existing permafrost OM mineral clusters, the litter-C would be more readily available for microbial decomposition. This model could also explain why we found, contrary to our expectations, a positive correlation between litter-C decomposition parameters and the proportion of litter-C bound to the mineral fraction, which indicates that also relatively labile litter-C is present in the mineral-associated OM fraction.

4.4 Significance of laboratory data for in situ carbon turnover

Our nine-year lasting incubation study at constant laboratory conditions aimed at elucidating some of the basic processes of fresh plant litter stabilization in thawing permafrost material. Our Study was not designed to mimic in situ environmental conditions, which constantly change in the Arctic tundra, or to consider changes in vegetation composition, which are expected to occur with Arctic warming. We chose a constant incubation temperature of 4 °C to represent summer temperatures in the active layer of the Siberian sampling sites. However, the temperature in the active layer on Samoylov may fluctuate between about +10 and −20 °C over the year with less than 4 months of non-freezing conditions (Boike et al., 2019). Hence; carbon decomposition will likely be higher under constantly unfrozen conditions in the laboratory then during the long freezing period in situ. On the other hand, annual freeze-thaw cycles in the tundra cause the liberation of labile OM resulting in peak carbon decomposition rates after such thaw-events (Walz et al., 2017), which was not considered in our incubations and may cause an underestimation of in situ OM decomposition. Furthermore, we simulated the most extreme conditions in terms of oxygen supply, either permanently fully oxic or fully anoxic while the oxygen concentrations and correspondingly the redox potential fluctuate in tundra soils during the year, e.g., after snow melt or during the deepening of the active layer (Liebmann et al., 2025). Such changes in oxygen concentrations will directly affect microbial OM decomposition. Additionally, decreasing redox potentials support the microbial reduction of iron(III)(oxoyhydr)oxides and the subsequent release of labile OM bound to the iron minerals (Patzner et al., 2020; Dong et al., 2023). In contrast, at rising redox potentials dissolved OM may adsorb and co-precipitate with newly formed iron(III)(oxoyhydr)oxides (Jia et al., 2024) thereby reducing the availability of degradable OM. We also did not consider the dispersal of microbial communities into thawing permafrost (Monteux et al., 2020) or priming effects, although priming seems of minor importance in our permafrost samples (Knoblauch et al., 2018). Considering these limitations of our laboratory experiments, we cannot provide robust estimates on the turnover time of fresh plant litter in thawing permafrost, i.e., the average time that litter-C remains in the soils before it is decomposed to CO₂ or CH₄. Particularly the model results from the anoxic incubations are less robust since the observational data from several replicates could not be used for model calibration.

Despite these limitations of our controlled laboratory experiment, our study sheds new light and mechanistic understanding on the fate of fresh plant litter in thawing permafrost soils. We provide first evidence, that thawing permafrost not only contributes to the release of greenhouse gases after thaw but may also function as a long-term carbon sink for fresh plant litter. This mechanism is particularly relevant in hillslope abrupt thaw features, such as those observed at the Kurungnakh study site, where abrupt permafrost thaw causes the active layer to erode, exposing deep permafrost material to the surface within a short period of time. Such hillslope erosional features may become strong sources of greenhouse gases while they are active (Knoblauch et al., 2021; Cassidy et al., 2016; Mu et al., 2017) and contribute substantially to the greenhouse gas release form the permafrost zone (Turetsky et al., 2020). However, the regrowth of vegetation stabilizes these thaw features (Huebner et al., 2022), and the stabilization of fresh plant litter in these soils may turn them from a carbon source into a carbon sink (Wickland et al., 2020).

5 Conclusions

The current study demonstrates that thawing permafrost material is not only a source of organic carbon that may be decomposed to the greenhouse gases CO₂ and CH₄ but that it may also stabilize a substantial fraction of recent plant litter over decades mainly in the mineral associated fraction. By the end of the experiment, a higher fraction of OM from added litter was bound to the mineral fraction than of OM from the thawed permafrost material. Nevertheless, the decomposability of OM from recent plant litter is one order of magnitude higher than that of permafrost carbon. These data suggest that the mineral-associated OM in our incubations comprises OM with low decomposability from the original thawed permafrost and OM with higher decomposability from the added plant litter. The higher decomposability of litter derived OM in the mineral fraction may be explained by a preferential adsorption of fresh litter derived OM on pre-existing OM on the mineral surfaces of the thawed permafrost material, thereby forming an outer layer of OM most available to microbial decomposition.

Furthermore, we found no evidence that particulate OM from the thawed permafrost may be faster decomposed than the mineral bound OM, which challenges the widely used classification of particulate OM as most labile and the mineral bound OM as stable. At least in case of permafrost organic matter, further parameters apart from mineral interactions contribute to organic matter stabilization. Our study emphasizes that the carbon balance of thawing permafrost soils depends not only on the release of old organic carbon from former permafrost, but also on its capacity to stabilize fresh organic matter from increasing plant litter input due to higher primary productivity in a warmer Arctic. This is particularly relevant in landscapes affected by abrupt thaw features, where surface vegetation is eroded and deep permafrost is exposed, allowing new vegetation to establish over time and stabilize the soil surface. However, to decipher the carbon balance of thawing permafrost landscapes, a combination of field observations, more complex laboratory experiments and process models is required, which consider the capacity of thawing permafrost to stabilize fresh plant litter.

Appendix A

Figure A1 Timeline of the incubation experiment. Left panel: Pre-incubation of permafrost samples for four years under oxic and anoxic conditions at 4 °C until constant CO₂ and CH₄ production rates have been established. Data of this pre-incubation experiment are not considered in this manuscript. At the end of the pre-incubation ¹³C-labelled Carex aquatilis litter was added. Right panel: (a) Main incubation experiment. Cumulative CO₂ and CH₄ production under oxic and anoxic conditions from permafrost-C and litter-C in a permafrost sample from the Siberian Island Kurungnakh (0.9–1.1 m depth) incubated at 4 °C for nine years. The left axis shows CO₂ and CH₄ produced from the two carbon sources (litter-C, permafrost-C). Black line represents the simulation of the calibrated two-pool carbon decomposition model for CO₂ (solid) and CH₄ (dashed). (b) Density fractionation of the organic carbon into the dissolved (DOM), the particulate (POM) and the mineral associated organic matter (MAOM) fraction at the end of the experiment. (c) Partitioning of the remaining carbon in the three fractions into the two different carbon sources (permafrost-C, litter-C) using the δ¹³C-signatures in the respective pools.

Figure A2 Partial dependence plots for Random Forest models, showing the responses of predictors used in Fig. 4. Panels show the fraction of litter-C decomposed (a), the mean residence time (MRT) of stable litter-C pool (b), the fraction of permafrost-C decomposed (c) and the MRT of the stable permafrost-C pool (d) under oxic and anoxic conditions. Note that variables were scaled prior to analyses using z-score centring.

Figure A3 Contribution of the dissolved carbon fraction (DOC) to total litter-C and total permafrost-C recovered at the end of the experiment under oxic and anoxic incubation conditions. Significant differences are indicated by *** = p<0.001 (Kruskal-Wallis ANOVA). The Box shows the $\mathrm{25}/\mathrm{75}$ percentiles, the whiskers the $\mathrm{5}/\mathrm{95}$ percentiles, the horizontal line the median and the white square the mean. Outliers are presented as circles.

Table A1 Characteristics of sampled permafrost deposits.

* data from Knoblauch et al. (2018).

Table A2 Selected model parameters based on stepwise linear regression. Df = degrees of freedom, AIC = Akaike Information Criterion. For further analysis, models for aerobic and anaerobic samples were merged to include the same set of predictors, and due to the correlation between the different soil fractions, the clay fraction was included in all models to represent soil texture.

Table A3 Fraction of total carbon decomposed during the first incubation year and final carbon decomposition rates standardized to the amount of carbon of the respective pool remaining at the end of the incubations.