Moisture and temperature effects on the radiocarbon signature of respired carbon dioxide to assess stability of soil carbon in the Tibetan Plateau

Microbial release of CO₂ from soils to the atmosphere reflects how environmental conditions affect the stability of soil organic matter (SOM), especially in massive organic-rich ecosystems like the peatlands and grasslands of the Qinghai–Tibetan Plateau (QTP). Radiocarbon (¹⁴C) is an important tracer of the global carbon cycle and can be used to understand SOM dynamics through the estimation of time lags between carbon fixation and respiration, often assessed with metrics such as age and transit time. In this study, we incubated peatland and grassland soils at four temperature (5, 10, 15 and 20 °C) and two water-filled pore space (WFPS) levels (60 % and 95 %) and measured the ¹⁴C signature of bulk soil and heterotrophic respired CO₂. We compared the relation between the Δ¹⁴C of the bulk soil and the Δ¹⁴CO₂ of respired carbon as a function of temperature and WFPS for the two soils. To better interpret our results, we used a mathematical model to analyse how the calculated number of pools, decomposition rates of carbon (k), transfer (α) and partitioning (γ) coefficients affect the Δ¹⁴C bulk and Δ¹⁴CO₂ relation, with their respective mean age and mean transit time. From our incubations, we found that ¹⁴C values in bulk and CO₂ from peatland were significantly more depleted (old) than from grassland soil. Our results showed that changes in temperature did not affect the Δ¹⁴C values of heterotrophic respired CO₂ in either soil. However, changes in WFPS had a small effect on the ¹⁴CO₂ in grassland soils and a significant influence in peatland soils, where higher WFPS levels led to more depleted Δ¹⁴CO₂. In our models, the correspondence between Δ¹⁴C, age and transit time highly depended on the internal dynamics of the soil (k, α, γ and number of pools) as well as on model structure. We observed large differences between slow and fast cycling systems, where low values of decomposition rates modified the Δ¹⁴C values in a non-linear pattern due to the incorporation of modern carbon (¹⁴C bomb) in the soil. We concluded that the stability of carbon in the peatland and grassland soils of the QTP depends strongly on the direction of change in moisture and how it affects the rates of SOM decomposition, while temperature regulates the number of fluxes. Current land cover modification (desiccation) in Zoigê peatlands and climate change occurring on the QTP might largely increase CO₂ fluxes along with the release of old carbon to the atmosphere potentially shifting carbon sinks into sources.

1 Introduction

Studying soil organic matter (SOM) stability and persistence in a globally changing environment is of fundamental importance to understand temporal variations of carbon cycling in the earth–climate system. Soil constitutes the largest carbon (C) stock in the terrestrial biosphere (Chen et al., 2021), making it a key component in global climate models (McGuire et al., 2001; Wieder et al., 2013; Xu et al., 2016). Physical, chemical and biological properties determine SOM decomposition and in consequence its persistence over decades to millennia (Schmidt et al., 2011). Temperature and moisture are two of the most important abiotic variables controlling the rates of SOM decomposition (Sierra et al., 2015) and its transit time across the different ecosystem pools. Changes in these environmental conditions are, however, occurring simultaneously, highlighting the necessity to conduct multifactorial experiments to disentangle the dominant mechanisms on such cycling times. Hence, discerning the relationship between SOM persistence and carbon cycling times in terrestrial ecosystems is imperative to improve climate change models and to inform decisions on CO₂ mitigation and management strategies (Bradford et al., 2016; Mesfin et al., 2021).

The Qinghai–Tibetan Plateau (QTP) with an area of 2.5×10⁶ km² is a macroregion of global importance for the cycling of water, carbon and other biogeochemical elements (Anslan et al., 2020). Land cover is mainly dominated by alpine grasslands (44 %) (Scurlock and Hall, 1998; Gao et al., 2014), which store 23.4 % of China’s total organic carbon and 2.5 % of the global soil carbon (Genxu et al., 2002). The QTP also hosts one of the largest high-mountain marshes in the world (Xiang et al., 2009; Chen et al., 2014; Ma et al., 2016), the Zoigê peatlands, considered to be the most important carbon stock in peatlands for China (Liu et al., 2018). However, climate change and land cover change are currently modulating the net carbon balance through CO₂ and CH₄ effluxes from the grasslands (Piao et al., 2012; Chen et al., 2017; Du and Gao, 2020) and the peatlands at both daily and interannual timescales (Hao et al., 2011; Chen et al., 2014, 2021; Kang et al., 2018; Liu et al., 2019b). Anthropogenic activities in the Zoigê peatlands (drainage, peat mining and overgrazing) have caused a degradation of approximately 30 % of wet and dry meadows (Chen et al., 2014; Zhou et al., 2021). Additionally, the QTP has been facing an air temperature increase of 0.2 °C per decade over the past 50 years (Zhang et al., 2013; Chen et al., 2014; Yang et al., 2014; Ganjurjav et al., 2016), equivalent to 2–3 times faster than the world average (Yao et al., 2019; Nieberding et al., 2020), along with a moderate increase in precipitation (Dong et al., 2018). According to the IPCC (Arias et al., 2021), climate models predict an increment in heavy precipitation events and high temperature extremes for this region.

Such changes in temperature and soil moisture control the magnitude in which stabilization and destabilization mechanisms enable carbon storage in the QTP (Xiang et al., 2009; Ma et al., 2016). It has been observed that warming increases soil respiration (Rustad et al., 2001; Lu et al., 2013; Pold et al., 2015) and that soil moisture modulates ecosystem and soil respiration (Geng et al., 2012; Piao et al., 2012; Moyano et al., 2013; Pan et al., 2022; Azizi-Rad et al., 2022). Although the influence of destabilization mechanisms on SOM decomposition has been already studied in the Zoigê peatlands (Zhao et al., 2011; Wang et al., 2015; Liu et al., 2016, 2018; Li et al., 2018; Liu et al., 2019a), temporal scales of SOM persistence under temperature and soil moisture variations are not clearly understood yet. Therefore, a quantification of the SOM persistence is of vital importance to predict climate change feedback magnitudes and pathways.

A useful approach to quantify SOM persistence is through the theory of compartmental dynamical systems, where soil is understood as a set of interconnected pools with transformations of carbon to different forms (gas, dissolved or solid) (Emanuel et al., 1984; Schimel, 1995; Sierra et al., 2018a). Therefore, carbon can be characterized by the time it remains inside a compartment or the entire system (Eriksson, 1971; Bolin and Rodhe, 1973) through the calculation of age and transit time. Here, age is defined as the time elapsed since the carbon entered the system until the time of observation in the bulk soil, while transit time is defined as the time spent by the carbon between the entry to the system and its exit as respired CO₂ (Eriksson, 1971; Bolin and Rodhe, 1973; Manzoni et al., 2009; Sierra et al., 2018b). These timescale metrics are, however, not the same for all the carbon atoms in a soil, since physical, chemical and biological stabilization and destabilization processes of SOM modify the rate at which mixing and storage occurs.

Radiocarbon measurements in bulk soil and in respired CO₂ are a powerful tool to approximate ages and transit times of carbon in soils as they trace back the trajectory of carbon through the different stocks on decadal to millennia timescales (Trumbore, 2000; Sierra et al., 2014; Schuur et al., 2016; Estop-Aragones et al., 2020; Chen et al., 2021). The radiocarbon content reflects the time a carbon atom has been in the soil since it was fixed by photosynthesis from the atmosphere (Trumbore, 2000). Carbon fixed since the early 1960s has higher ¹⁴C than carbon fixed previously due to the enrichment by thermonuclear weapon testing. By comparing the ¹⁴C relative abundance in the carbon pools with the atmospheric ¹⁴CO₂ concentrations, it is possible to model soil–atmosphere carbon cycling times (Trumbore, 2000).

¹⁴C values can be altered by SOM destabilization processes and soil characteristics such as soil organic carbon (SOC) content, age and diagenetic state, thaw depth, redox state, seasonality, etc. (Gaudinski et al., 2000; Trumbore, 2000; von Lützow et al., 2008; Sierra et al., 2018b; Estop-Aragones et al., 2020). In recent decades, it has been widely observed that temperature plays a major role on the dynamics of SOM (Knorr et al., 2005; Davidson and Janssens, 2006; Feng and Simpson, 2008) and the mean age of respired CO₂ (Hopkins et al., 2012; Chen et al., 2021) by increasing the decomposition rates from fast-cycling pools (Trumbore et al., 1996) and mobilizing old carbon in warming conditions (Dutta et al., 2006; Briones et al., 2010). In contrast, other studies have suggested that warming does not lead to release of old carbon (Briones et al., 2021; Dioumaeva et al., 2002) and that non-labile carbon decomposition is insensitive to temperature increase (Liski et al., 1999). By comparison, drying phenomena increased the release of modern ¹⁴CO₂ from shallow soil layers but preserved the old soil carbon pools in deeper layers (Kwon et al., 2019). Observations in rich-SOC soils of the Arctic indicated the release of old carbon from deep soil layers after thaw and drainage processes (Schuur et al., 2009; Estop-Aragones et al., 2020; Pegoraro et al., 2021). Such a process might be occurring in the Zoigê peatland soils due to the presence of seasonal frozen layers (Liu et al., 2021; Yang et al., 2022). Nonetheless, the influence of temperature and soil moisture on SOM persistence is a matter of current debate.

Changes in temperature and moisture contribute to the destabilization of carbon in soils from the QTP. Hence, we hypothesize that higher temperature would increase the age of respired CO₂, and changes in soil moisture would increase or decrease (depending on the direction of moisture change) the age of respired CO₂ in soils subjected to controlled manipulations. Greater ages of respired CO₂ would imply that previously stabilized carbon is destabilized by changes in the manipulated environmental factors. For this purpose, we posed three specific research questions:

1. Do specific changes in temperature and moisture result in the release of old carbon in the respiration flux?

2. Are there differences in the age of heterotrophic respired CO₂ between a grassland and a peatland soil from the QTP?

3. How can radiocarbon data in bulk soil and heterotrophic respired CO₂ be interpreted to understand the effect of changes in decomposition rates on ages and transit times?

To address these questions, we conducted a controlled multifactorial experiment with soils from two different ecosystems in the QTP and measured the radiocarbon signature in the bulk soil and the heterotrophic respired CO₂. In addition, we used a mathematical model to better interpret the interaction between decomposition rates change, expressed through internal dynamics of the soil (k, α, γ and number of pools) and the Δ¹⁴C values by targeting the range found in the incubations. Thus, our observations and models strengthen each other to gain a deeper comprehension of the relationship between soil carbon stability, Δ¹⁴C, age and transit time.

2 Materials and methods

To evaluate the influence of soil moisture and temperature on the ¹⁴C values from bulk soil and heterotrophic respired CO₂ of grasslands and peatlands soils, we incubated soil samples from a high-elevation grassland at the Nam Co catchment collected in September 2018 and a peatland from the Zoigê region sampled in July 2021, both located in the QTP.

2.1 Site description and sampling

2.1.1 Nam Co grassland

The Nam Co site (Fig. 1) is located in the central Tibetan Plateau (altitude of the Nam Co lake 4726 m a.s.l., 30°46^′ N, 90°59^′ E) and represents a frequent study location for monitoring and tracking of environmental changes over different timescales (Anslan et al., 2020). The dominant climate at Nam Co consists of cold winters and short and moist summers. Mean annual temperature (MAT) measured at the NAMORS research station was −0.6 °C (2006–2017), and mean annual precipitation (MAP) was 406 mm, occurring mostly during the monsoon season between May and October (Anslan et al., 2020). The sampling area was mainly covered by Kobresia pygmaea and has been grazed intensively by yaks and sheep. A total amount of 40 kg of soil was collected from randomly selected spots at depths between 5 and 15 cm within an area of about 40 m². The soil was then mixed homogeneously to produce one single sample which consisted of a sandy loam with a bulk density of 1.3 g cm⁻³, pH of 7.5, 3.2 % of total organic carbon, 0.05 % of inorganic carbon and 0.3 % of total nitrogen, cation exchange capacity of 89 mmol_c kg⁻¹ soil, and electrical conductivity of 89.6 µS cm⁻¹.

Figure 1 (a) Overview of the QTP: (b) location of the Nam Co lake with land cover dominated by grasslands and (c) Zoigê peatlands highlighting the location of the city of Zoigê. The sampling site is located in the southeast of the peatland (Ruokeba, Hongyuan County). SRTM elevation data from Jarvis et al. (2008).

2.1.2 Zoigê peatlands

The Zoigê peatlands are located in the northeast of the QTP (average elevation of 3400 m a.s.l.) and cover an area of about 4605 km² on the headwaters of the Yellow River basin. These peatlands have been recognized as one of the largest high-mountain marshes in the world (Xiang et al., 2009; Ma et al., 2016) and store 88 % of the carbon in the QTP (Chen et al., 2014). Based on multi-decadal records, the MAP was 720 mm and the MAT was 1.5 °C, where the warmest monthly temperature (11 °C) is recorded in July and the coldest in January (−10.1 °C), showing a slight drying trend and a warming of 0.4 °C per decade since 1970 (Chen et al., 2014; Yang et al., 2014). The soils are unfrozen from April to October, while the layer between 0 and 50 cm depth is seasonally frozen between November and March. Vegetation cover consists mainly of Potentilla anserine, Blysmus sinocompressus, Kobresia myosuroides and Scirpus triqueter. The samples (0–35 cm depth) were collected with a spade from an area of 50 m² at a long-term monitoring site (33°4^′5^′′ N, 102°33^′52^′′ E) and then thoroughly mixed to produce a single sample. The soil bulk density was 0.3 g cm⁻³ and contained on average 27.6 % organic carbon, 0.06 % inorganic carbon and 1.8 % total nitrogen.

2.2 Incubation of peatland and grassland soils

We conducted two sets of incubation experiments without pre-incubation period: one set with the grassland soil and a second set with the peatland soil. Samples were incubated with their original roots to minimize disturbance and allow comparisons with field conditions; however stones were removed. Briefly, 50 mL vials were filled with 12 g of soil (±1.5 g) and placed inside glass flasks along with 0.2 mL of water at the bottom of the flask (away from contact with the sample) to avoid possible drying (Dioumaeva et al., 2002). Thereafter the flasks were sealed with rubber plugs and screwed with plastic caps. Each of the sets was placed at two different water-filled pore space (WFPS) levels (60 % and 95 %), which were selected in order to reassemble the thaw and consequent water saturation of seasonally frozen soils in grasslands and for peatland soils the process of drying (through artificial desiccation) after high water saturation. Despite the high WFPS, the soil samples still had contact with air inside the vials, which guaranteed microbial decomposition of the organic matter and accumulation of heterotrophic CO₂ respiration. Soil moisture treatments were combined with four temperature levels (5, 10, 15 and 20 °C) for a total of 69 samples: 33 for grassland and 36 for peatland (2 sites × 4 temperature levels × 2 moisture levels × 3–6 analytical replicates, to prevent scarcity of data due to eventual failure in CO₂ extraction). Flasks with samples were flushed with synthetic air (CO₂ free) to remove atmospheric CO₂. This flushing marked the starting day of the incubations.

Incubations for each subset concluded concurrently once all the samples reached a C concentration (from CO₂) in the headspace estimated to be equal to or exceeding 2 mg, sufficient for subsequent radiocarbon analysis. This approach was not possible in 2 of the 16 subsets due to lab material limitations; therefore, grassland samples were incubated between 15 and 67 d, while peatland samples were incubated for 13 d (Table A1). The headspace volume of the incubation flasks was measured as 587 mL, which was corrected after adding the soil. On average, the final headspace was 575.6 and 533 mL for the incubated grassland and peatland soils, respectively. These values were used to calculate the fluxes of heterotrophic CO₂ respiration.

The rate of accumulation of CO₂ in the headspace of our incubations represents the diffusion rate of heterotrophic CO₂ respiration released from the incubated soils. We refer throughout the article to heterotrophic respiration, but we acknowledge that our measurements better capture how this heterotrophic CO₂ respiration flux diffuses out of the soil. Rates were measured at intervals of 1 to 2 weeks using a CO₂ analyser LI-COR 6262 for every treatment, and mean heterotrophic CO₂ respiration rates (mg CO₂ g soil⁻¹ d⁻¹) were calculated through the division of CO₂ concentration in the headspace by the product of the accumulation duration (days) (Table A1) and the mass of the introduced soil (g).

2.3 Radiocarbon analysis of incubated soils

Radiocarbon analyses were conducted in the bulk soil of each sample after the incubation. Soil inorganic carbon was eliminated through decalcification following the preparation protocol by Steinhof et al. (2017). Also, heterotrophic CO₂ respiration accumulated in the headspace of incubation flasks was extracted and purified on a vacuum line, graphitized by Fe reduction in H₂ and measured for Δ¹⁴C by an accelerator mass spectrometer (MICADAS, Ionplus, Switzerland) in the radiocarbon laboratory of the Max Planck Institute for Biogeochemistry in Jena, Germany (Steinhof et al., 2017).

Radiocarbon data are expressed as Δ¹⁴C (the deviation in ‰ from 0.95 times the oxalic acid standard in 1950). The values were corrected to a δ¹³C value of −25 ‰ for differences in biological mass-dependent fractionation (Stuiver and Polach, 1977). These data are presented as percentage of modern carbon (pMC) that can be converted to F¹⁴C by dividing it by 100 and later to Δ¹⁴C using Eq. (1) (Stuiver and Polach, 1977).

$$\begin{array}{}\text{(1)}& {\mathrm{\Delta }}^{\mathrm{14}}\text{C}=\left[{{\mathrm{F}}^{\mathrm{14}}\mathrm{Ce}}^{{\mathit{\lambda }}_{\text{C}}(\mathrm{1950}-t)}-\mathrm{1}\right]\times \mathrm{1000}\left[\mathrm{‰}\right],\end{array}$$

where F¹⁴C is the fraction modern, i.e. the ratio of the measured sample normalized to a δ¹³C value of −25 ‰, divided by 0.95 times the measured ratio of the oxalic acid I standard (OX-I) (Schuur et al., 2016), λ_C is the updated radiocarbon decay constant (equals $\mathrm{1}/\mathrm{8267}$ [yr⁻¹]), and t is the year of sampling.

The effects of soil moisture and temperature manipulation on Δ¹⁴C were evaluated for the bulk and CO₂ fractions in each ecosystem separately through two-way ANOVA tests (type III). This type of ANOVA is also referred to as the partial sum of squares and is appropriate for unbalanced data since it does not depend on the sampling structure or the particular order in the model (Shaw and Mitchell-Olds, 1993); hence, this approach adjusts best to our data set where treatments did not have equal number of values.

2.4 SOC decomposition models to predict Δ¹⁴C as a proxy of SOM persistence

The representation of the SOC dynamics has been commonly described through models that can be expressed as systems of linear differential equations (Manzoni et al., 2009; Sierra et al., 2012, 2017b; Sierra and Mueller, 2015) of the following form:

$$\begin{array}{}\text{(2)}& {\displaystyle \frac{\mathrm{d}\mathit{C}\left(t\right)}{\mathrm{d}t}}=\mathit{I}+\mathbf{A}\cdot \mathit{C}\left(t\right),\end{array}$$

where the vector C(t) is the rate of change of carbon over time in n different pools; the time dependent n-dimensional vector I represents the total input of carbon to each pool; and A represents the n×n-dimensional matrix with the rates of carbon processing for each pool in its main diagonal and the proportion of carbon transferred from one pool to another in the off-diagonals (Sierra et al., 2012, 2014, 2018b; Metzler and Sierra, 2018). This mass balance equation has a radiocarbon counterpart:

$$\begin{array}{}\text{(3)}& {\displaystyle \frac{{\mathrm{d}}^{\mathrm{14}}\mathit{C}\left(t\right)}{\mathrm{d}t}}={\mathit{I}}_{{}^{\mathrm{14}}\mathbf{C}}\left(t\right)+\mathbf{A}{\cdot }^{\mathrm{14}}\mathit{C}\left(t\right)-{\mathit{\lambda }}^{\mathrm{14}}\mathit{C}\left(t\right),\end{array}$$

where λ is the radiocarbon decay constant ($\mathrm{1}/\mathrm{8267}$ [yr⁻¹]).

We used a SOC decomposition model to calculate the Δ¹⁴C in bulk and CO₂ as well as their equivalent mean age and mean transit time of a theoretical soil for the year 2019 calibrated with the data set for the Northern Hemisphere Zone 3 (Hua et al., 2021). The simulated Δ¹⁴C results aimed at finding the same Δ¹⁴C ranges obtained from the incubated soils establishing different ranges of decomposition rates for peatland and grassland, with the objective to compare under which settings the models could describe best laboratory or field conditions and for this reason, model fitting was not necessary. For this purpose, we implemented a two-pool model (one slow and one fast pool) considering two different structures – parallel and series (Fig. 2) – and modified the parameters involved in the soil carbon processing (Manzoni et al., 2009; Falloon and Smith, 2000) (Table 1). In the parallel structure, carbon enters the soil and splits among the two pools according to γ and decomposition occurs independently in each pool according to their respective k. In the series structure, carbon enters only to the fast pool and it is either decomposed and emitted to the atmosphere or transferred to the slow pool according to the transfer coefficient α.

Figure 2 Parallel (a) and series (b) model structure and formulas implemented in the simulations. Carbon is fixed into the biomass from atmospheric CO₂ through photosynthesis. Subsequently, it is incorporated into the soil as litterfall. From this state on, it is considered input (I) and will split between the two pools according to the different structures (γ for parallel and α for series). The rate at which the carbon is decomposed (k) in each pool (C) will depend on microbial activity, environmental factors (temperature and soil moisture), and physical and chemical protection of the SOM (Blanco-Canqui and Lal, 2004; von Lützow et al., 2008; Manzoni et al., 2009). Boxes represent the soil pools, which account with an initial amount of carbon (C₁=200, C₂=5000), and arrows represent directions of input and output.

Initial parameters of the model such as the starting year of simulation and initial Δ¹⁴C values of each pool were considered separately for each type of ecosystem. Moreover, we defined I and C as constant since they can be adjusted depending on the specific soil characteristics. Additionally, we assumed that decomposition rates k reflect the effect of temperature and soil moisture on carbon cycling timescales (Manzoni et al., 2009) and therefore on Δ¹⁴C values. In our approach, only one parameter at time was modified for each simulation.

Table 1 Definition of parameters used to evaluate the variation of Δ¹⁴C values, mean age, and mean transit time and their ranges used in the simulation experiments.

We used the R package SoilR (Sierra et al., 2014) to simulate the temporal dynamics of Δ¹⁴C in the bulk soil and the respired CO₂ as well as the age and transit time distributions of carbon. Assuming steady state for the carbon stocks, the probability density function (PDF) of the age (Eq. 4) and the transit time (Eq. 5) (Metzler and Sierra, 2018) as well as their means can be calculated by the following expressions (Sierra et al., 2018b):

$$\begin{array}{}\text{(4)}& f\left(a\right)=-{\mathrm{1}}^{\text{T}}\cdot \mathbf{A}\cdot e^{a\cdot \mathit{A}}\cdot {\displaystyle \frac{{\mathit{C}}^{*}}{\sum {\mathit{C}}^{*}}},a\ge \mathrm{0},\end{array}$$

where a is the random variable age, 1^T is the transpose of the n-dimensional vector containing 1, e^a⋅A is the matrix exponential for each value of a, and $\sum {\mathbf{C}}^{*}$ is the sum of stocks of all pools at steady state.

The PDF of transit times can be obtained as

$$\begin{array}{}\text{(5)}& f\left(\mathit{\tau }\right)=-{\mathrm{1}}^{\text{T}}\cdot \mathbf{A}\cdot e^{\mathit{\tau }\cdot \mathit{A}}\cdot {\displaystyle \frac{\mathit{I}}{\sum \mathit{I}}},\mathit{\tau }\ge \mathrm{0},\end{array}$$

where τ represents the random variable transit time.

These PDFs measure the probability that a certain amount of carbon is above or below a specific age or transit time (Sierra et al., 2018b) (Eqs. 6 and 7, respectively). The mean of the age (i.e. the expected value of the PDF) can be computed by the following expression:

$$\begin{array}{}\text{(6)}& \mathbb{E}\left(a\right)=-{\mathrm{1}}^{\text{T}}\cdot {\mathbf{A}}^{-\mathrm{1}}\cdot {\displaystyle \frac{{\mathit{C}}^{*}}{\sum {\mathit{C}}^{*}}},\end{array}$$

The mean of the transit time can be computed by

$$\begin{array}{}\text{(7)}& \mathbb{E}\left(\mathit{\tau }\right)=-{\mathrm{1}}^{\text{T}}\cdot {\mathbf{A}}^{-\mathrm{1}}\cdot {\displaystyle \frac{\mathit{I}}{\sum \mathit{I}}}.\end{array}$$

Based on the age and transit time distribution for the different simulated cases, metrics such as the mean, median and quantiles can be used as proxies of SOM persistence. The relation between ages and transit times may present three cases: 𝔼(a)=𝔼(τ) (type I – well-mixed homogeneous system), 𝔼(a)<𝔼(τ) (type II – retention system) and 𝔼(a)>𝔼(τ) (type III – non-retention system) (Sierra et al., 2018b). In type I, the probability of mineralization and release as CO₂ is the same for every C atom. In type II, the carbon is retained for a relatively long time before it is released. In type III, most of the C atoms stay in the system for a short period of time, but some atoms remain for a long time (Bolin and Rodhe, 1973). Comparisons between these distributions will provide detailed information on mixing, store, recycling, transport and transformation processes of the SOM (Sierra et al., 2018b).

Figure 3 Boxplot with the variation of the Δ¹⁴C values in CO₂ (a) and bulk soil (b) for peatland and grassland incubated soils.

3 Results

3.1 Δ¹⁴C values of incubated soils under temperature and soil moisture variation

Values of Δ¹⁴C from the incubation experiment contrasted strongly between grassland and peatland soils (Figs. 3 and 4). Generally, Δ¹⁴C of bulk soil and heterotrophic CO₂ respiration from the peatland soil was more depleted than from the grassland soil. For example, in peatland soil the Δ¹⁴C values of bulk soil (−85.9 ‰ to −60.9 ‰, mean = −80.1, n=36, including outlier of −180) were clearly more depleted than those of the heterotrophic respired flux (−20.9 ‰ to 23.9 ‰, mean = 4, n=36), which indicated that the peatland behaved as a retention system (type II). In contrast, for the grassland soil, the Δ¹⁴C of bulk soil (21.1 ‰ to 73.9 ‰, mean = 43.3, n=33) fell similarly around the 1:1 line compared to the Δ¹⁴CO₂ (13.9 ‰ to 83.4 ‰, mean = 38.5, n=33, including outliers of −227.1 and −105.1), indicating that the samples behaved mostly as a well-mixed homogeneous system (type I). Such similar values can also reflect a type II system where Δ¹⁴C values in the bulk and the respired flux are equal for the year of sampling (Fig. A7).

The temperature treatments did not systematically affect the radiocarbon signature of the bulk or the heterotrophic respired CO₂ in peatland or grassland soils (Table 2). There was significant evidence that manipulations in WFPS resulted in changes in the Δ¹⁴C values of bulk soil and CO₂ (p values = 0.09 in grasslands and p value = 0.01 in peatlands from an ANOVA test, Table 2, Fig. 5a–b) except for the bulk soil in peatlands. When the interacting effects of WFPS and temperature were evaluated together, there was no evidence that their interplay affected the radiocarbon signature of both soils (p values > 0.05 for all the analysis). Outliers presented extremely depleted Δ¹⁴C values and occurred both in bulk (peatland) and CO₂ (grassland) at the combined treatment WFPS = 60 and temperature = 10 °C. Such wide variation in Δ¹⁴C both between ecosystems and treatments could be potentially explained by intrinsic processes affecting soil carbon dynamics, which will be explored in the following sections with a SOC decomposition model.

Figure 4 Relationship between Δ¹⁴C in bulk soil and Δ¹⁴C in heterotrophic respired CO₂ for incubated grassland (green box) and peatland (brown box) soils of the QTP discriminated by temperature and WFPS. Possible types of system according to the Δ¹⁴C relations between bulk and respired carbon.

Higher temperature and WFPS caused an increase of CO₂ fluxes from heterotrophic respiration in the treated incubated soils (Table A1, Fig. 5c–d). In both ecosystems, wetter conditions showed higher respiration rates and higher slopes as the temperature increased. The absolute amounts of CO₂ produced from peatland soil were on average 14 times higher than from grassland soils for every independent treatment.

Table 2 Summary of p values obtained from ANOVA tests for the Δ¹⁴C of bulk soil and the Δ¹⁴CO₂ of grassland and peatland soils after incubation. p values are given for the independent effect of temperature (T) and WFPS as well as the integrated effect of temperature and WFPS (T ⋅ WFPS).

Figure 5 Comparison between Δ¹⁴C of heterotrophic respiration from incubated grassland (a) and peatland (b) soil at different temperature levels under WFPS = 605 % and 95 %. Black points represent minimum and maximum values out of the range between quartile 1 and 3 (25 % to 75 % of the data). Quartile 50 (median) represented by the line inside the box indicates the midpoint value in the frequency distribution. The box for the treatment WFPS = 60 % and T=10 °C shows a large dispersion of the 50 % of the data, which is explained by the outliers observed in Fig. 4. Additional panels indicate the heterotrophic respiration rates (mg CO₂ g soil⁻¹ d⁻¹) for each treatment based on the duration of the flux accumulation (see Table A1). Response of mean heterotrophic respiration rates to temperature and WFPS treatments in grassland (c) and peatland (d) soils is also shown.

3.2 Effect of changes of decomposition rates (k) on the Δ¹⁴C values

Changes in the vertical and horizontal direction of the Δ¹⁴C bulk versus Δ¹⁴C–CO₂ space are more evident across ecosystem type, which at the same time implies specific environmental conditions for the stability of SOM. To understand possible drivers of Δ¹⁴C changes in these directions and how they can be interpreted in terms of ages and transit times, we used a SOC decomposition model. We evaluated how model structure, decomposition rates of carbon, and the partitioning (γ) and transfer (α) coefficients of a two-pool parallel and a series model affected Δ¹⁴C values and, as a consequence, the mean age and the mean transit time. Our simulations mimicked a fast cycling grassland and a slow cycling peatland by differentiating the ranges of decomposition rates (Table 3) and showed how the modelled conditions affected the type of system (I, II and III).

Table 3 Range of parameters used for simulations in a SOC decomposition model for fast and slow systems. These ranges explore the required conditions to describe the Δ¹⁴C variation found in the incubated grassland and peatland soils and their equivalent age and transit time relationship. The target Δ¹⁴C intervals are shown as boxes in Figs. 6 to 9.

3.2.1 SOC decomposition in fast cycling systems (grasslands)

Variation of k₁ yielded Δ¹⁴C curves similar to mean age and mean transit time. In other words, in these simulations the change in the parameter k₁ resulted in similar trends in Δ¹⁴C bulk versus Δ¹⁴C–CO₂ space as in the mean age versus mean transit time space (Fig. 6). In the parallel structure (Fig. 6a, b), high values of k₁ yielded more enriched Δ¹⁴C values in the bulk soil than in the heterotrophic respired flux. This was expressed as a fast respiration of SOC and reflected in a short transit time. Simultaneously, low values of k₂ resulted in a slow respiration, which was registered in ages longer than transit times, at the initial state of the carbon stay in the system. However, as k₁ decreased and approached k₂, mean age and mean transit time became similar and converged to the 1:1 line. In the case of the series structure (Fig. 6c, d), the total amount of input travelled first through the fast pool, which meant that while k₁ is high, most of the carbon decomposition occurs under the dynamics of the fast pool. Hence, the system showed equal mean age and mean transit time. As k₁ decreased, the transfer from the fast to the slow pool becomes more relevant, and then the mean transit time increased compared to the mean age.

The behaviour of Δ¹⁴C values contrasted significantly with that of mean age and mean transit time when k₂ tended to lower values (Fig. 7). Generally, when k₂ equals k₁, Δ¹⁴C and mean age and mean transit time showed similar values. However, as k₂ decreased, the carbon stayed in the slow pool for a longer time and the Δ¹⁴C enriched in the bulk soil until reaching a peak (112 ‰ and 89 ‰ for series and parallel structures, respectively), from where it subsequently depleted to the initial Δ¹⁴CO₂ value. For the parallel structure (Fig. 7a, b), transit time kept increasing slowly since most of the carbon stayed in the fast pool for a short time, so the carbon remaining in the slow pool contributed to increasing the mean age. As for the series structure, most of the carbon was transferred from the fast to the slow pool due to the high transfer rate α, which contributed to increasing the transit time.

Figure 6 Predictions of Δ¹⁴C in bulk soil vs. Δ¹⁴CO₂ with their equivalent simulation of mean age in bulk soil vs. mean transit time in CO₂ for parallel (a, b) and series model structure (c, d) for a fast cycling system. Variation of k₁ with α=0.8 and γ=0.8. Complementary prediction with α=0.1 and γ=0.1 can be seen in Fig. A2. The green box represents the range of measured Δ¹⁴C values obtained from incubated grassland soils for bulk (20 ‰ to 74 ‰) and CO₂ (10 ‰ to 85 ‰), excluding outliers.

3.2.2 SOC decomposition in slow cycling systems (peatlands)

Simulations of Δ¹⁴C values in slow cycling systems resulted in significantly more depleted values in both bulk soil and respired CO₂ (Fig. 8) than in fast cycling systems. Generally, Δ¹⁴C values were more depleted in the bulk soil than in the respiration flux. Also, the results from these simulations showed very different patterns in the Δ¹⁴C bulk versus Δ¹⁴C–CO₂ space than in the mean age versus mean transit time space: only with the exception of equal k₁ and k₂, which resulted in similar Δ¹⁴C values as well as similar mean age and mean transit time in the parallel structure.

Figure 7 Predictions of Δ¹⁴C in bulk soil vs. Δ¹⁴CO₂ with their equivalent simulation of mean age in bulk soil vs. mean transit time in CO₂ for parallel (a, b) and series model structure (c, d) for a fast cycling system. Variation of k₁ with α=0.8 and γ=0.8. The green box represents the range of measured Δ¹⁴C values obtained from incubated grassland soils for bulk (20 ‰ to 74 ‰) and CO₂ (10 ‰ to 85 ‰), excluding outliers.

Additionally, we looked at the variation of Δ¹⁴C keeping the same k values but reducing α and γ to 0.1. Simulations indicated that the behaviour of Δ¹⁴C values was opposite between series and parallel structures when input to each pool was inversely proportional (Fig. 7 versus Fig. A3 and Fig. 8 versus Fig. A4). For example, for the case of low values of k₁ and k₂ (Fig. 9) (γ=0.2 and α=0.8), both model structures showed similar patterns due to relatively similar amounts of input going to the slow pool. When γ=0.8, Δ¹⁴C showed values out of the target range (−20 to 23 for Δ¹⁴C–CO₂ and −86 to −60 for Δ¹⁴C bulk based on incubation results) since decomposition occurred as in a fast cycling system. As k₂ decreased, Δ¹⁴CO₂ enriched and Δ¹⁴C and mean age remained higher than mean transit time since the slow decomposition of the slow pool dominated the system response.

3.3 Variation in the proportion of input (γ and α) modulated mean age and mean transit time

Our simulations showed the importance that the partitioning and transference coefficients (γ and α) have in defining the amount of carbon entering to each pool and in consequence the mean age and mean transit time (Fig. A1). For a parallel structure, when the total amount of input entered only to the fast pool (γ=1), the system behaved as a one pool system and mean age and mean transit time were equal. In the extreme opposite (γ=0), all carbon entered only to the slow pool and mean age and mean transit time converged, although with a longer transit through the system. In contrast, for a series structure, the extreme value α=1 resulted in a longer mean transit time than mean age as the contribution to the respired flux comprehended the sum of the decomposition of the newly absorbed carbon and the initial stocks stored in the fast pool. In the case of slow cycling systems (Figs. A5 and A6), Δ¹⁴CO₂ moved in a wider range including very depleted values (−100 ‰ to 70 ‰) and approached a steady state of Δ¹⁴C in bulk soil from the right when k₁=0.01 or from the left when k₁=0.1.

Figure 8 Predictions of Δ¹⁴C in bulk soil vs. Δ¹⁴CO₂ with their equivalent simulation of mean age in bulk soil vs. mean transit time in CO₂ for parallel (a, b) and series model structure (c, d) for a slow cycling system. Variation of k₁ with α=0.8 and γ=0.8. The brown box represents the range of measured Δ¹⁴C values obtained from incubated peatland soils for bulk (−90 ‰ to −65 ‰) and CO₂ (−18 ‰ to 25 ‰), excluding outliers.

4 Discussion

4.1 Did changes in temperature and soil moisture result in the release of old carbon in the heterotrophic respiration flux?

We measured the Δ¹⁴C values of heterotrophic respired CO₂ from soils incubated at different temperatures and WFPS levels to assess whether changes in these variables would result in destabilization of stored carbon. In general, higher temperature and WFPS levels resulted in larger CO₂ fluxes although with a faster increment in peatland soils, which could be explained by the high amount of total organic carbon available for decomposition. Heterotrophic respiration rates, however, did not correlate with the Δ¹⁴C values. Our results showed that changes in temperature did not affect systematically the radiocarbon values of heterotrophic respired CO₂ in any of the incubated soils. Nevertheless, changes in WFPS had a significant effect on the Δ¹⁴CO₂ and Δ¹⁴C bulk of grassland soils and only on the Δ¹⁴CO₂ of peatland soils. Our experiments indicated that higher WFPS levels led to depleted Δ¹⁴CO₂ values in peatlands at all temperature treatments, except at 10 °C.

In contrast, higher WFPS resulted in more enriched Δ¹⁴CO₂ values for grasslands. This suggests that the direction of changing WFPS depends on ecosystem characteristics. It is interesting to observe that the Δ¹⁴C in the bulk soil was always more negative than the Δ¹⁴CO₂ in the peatlands, while the grasslands registered similar Δ¹⁴C values in both states. This indicates that peatlands are systems that stabilize organic matter over time and release it once the stable conditions change. In comparison, for grasslands two interpretations may arise: first, that there was not organic matter stabilization under the incubated conditions since the respired flux reassembles the bulk soil Δ¹⁴C signature or, second, that the incubation conditions were not strong or long enough to destabilize the existing old SOM. Our simulations indicate that similar Δ¹⁴C values in the bulk and in the respired CO₂ can also occur at a specific year due to the dynamics of the atmospheric bomb curve (Fig. A7) independently of the different SOM cycling times in the soil pools.

In tundra ecosystems, Kwon et al. (2019) suggested that drainage of shallow soil layers may have accelerated old carbon decomposition. In addition, Estop-Aragones et al. (2020) concluded that old carbon would increase in proportion from “cold” across “warm wet” to “warm and dry” for high arctic tundra as well as from “undisturbed” to “burnt active layer” for peatland plateau. These previous findings are in disagreement with our observations since after our experiments, lower WFPS resulted in the release of relatively enriched Δ¹⁴CO₂ instead of releasing old carbon. One potential explanation for this discrepancy is that the existing old SOM could have reacted too slowly to drier conditions (WFPS = 60) compared to the fast response of “younger” SOM, causing an evident prevalence of young carbon given the short duration of the incubations for peatlands. Consistent with this, it has been found that some SOM components or compartments may be more sensitive to modified conditions than others (Feng and Simpson, 2008) and, therefore, transit time would react strongly when such sensitivity is higher in the slower pool (Manzoni et al., 2009).

For this reason, our results are indicative, but not conclusive, of the influence of environmental factors on the Δ¹⁴C signature of heterotrophic respired carbon. It is possible that our Δ¹⁴C values are limited by the short time of the incubations since the dynamics of carbon transfers cannot be properly observed in short timescales (Crow and Sierra, 2018). For example, during a 1-year incubation experiment, most of the CO₂ was derived from labile SOM as the temperature increased (Leifeld and Fuhrer, 2005). Additionally, other factors such as yearly climate seasonality, daily freeze–thaw cycles and water table oscillation were not replicated in the laboratory due to equipment and time availability. An important aspect that has been observed in soil incubations is that CO₂ accumulation decreases or even stops after a certain period probably due to the CO₂ saturation of the limited headspace in the incubation flasks, which depends on the SOC content. From this we could deduce that old carbon, which usually needs longer times to be decomposed, can only be recovered if soil respiration is not space limited. For example, Azizi-Rad et al. (2022) found the respiration rate to decline after 14 d holding the soil at 10 °C. In this sense, incubations with high total organic carbon may run out of headspace soon and affect the interpretation of the Δ¹⁴C values. Finally, although temperature has been found to be the main modulator of SOC decomposition rates (Azizi-Rad et al., 2022), dominating over the effect of soil moisture and oxygen availability, it is important to consider that an increase in respiration rates, as observed in our incubations, does not necessarily involve the release of Δ¹⁴C of lower or higher signatures as shown in our model simulations.

4.2 Are there differences in the age of heterotrophic respired CO₂ between grassland and peatland?

We found that Δ¹⁴C values from peatland were more strongly depleted than those in grasslands, indicating the presence and respiration of older carbon in peatlands. Generally, low temperature, low soil microbial activity and anoxic conditions (Xiang et al., 2009; Ma et al., 2016) have posed favourable conditions for the stabilization of organic matter since the Early Holocene in the Zoigê region (Chen et al., 2014; Sun et al., 2017). Older bulk soil than CO₂ demonstrated that the peatland soil behaved as a retention system (type II) where carbon was stored for relatively long time before its respiration (Sierra et al., 2018b). The release of that old carbon may occur when the SOM destabilization dominates over the stabilization mechanisms, for example when peatlands are drained (increasing oxidation through water table reduction) and their temperature is raised (Dutta et al., 2006; Hicks Pries et al., 2013; Lupascu et al., 2014; Estop-Aragonés et al., 2018), both factors mimicked in our incubations.

In contrast, Δ¹⁴CO₂ from grasslands were mostly similar to Δ¹⁴C in bulk soil, indicating that the soil behaved predominantly as a well-mixed homogeneous system (type I) where most of the SOM has the same probability of being mineralized and released as CO₂ (Sierra et al., 2018b). However, Δ¹⁴C relation also indicated that the grasslands shift from retention to non-retention system and vice versa. In grasslands, SOM from the topsoil is permanently under changing conditions (daily temperature and soil moisture variation, grazing and mechanical alteration, among others) that promote its fast cycling instead of its stabilization (Han et al., 2017). These factors, along with climatic, vegetation and edaphic properties, caused a low SOM stability expressed in a large proportion of labile carbon (Hou et al., 2021). Nonetheless, even though SOM accumulation still occurs in grassland areas (Tian et al., 2009) registered in the extremely depleted outliers from the grassland soil respiration, it happens at a slower rate than in peatlands. As a result, most of the fixed atmospheric carbon is respired at short timescales, which are in turn registered as young (“post bomb”) Δ¹⁴C values. Finally, although extremely depleted Δ¹⁴C values are rare, they are probable since the measured carbon particles represent one value, which can fall on the tail of the system–age probability distribution. We interpret that outlier results belong to carbon particles that have remained for very long time in the system (out of the mean values) whose specific soil characteristics were not captured in our model structure or internal soil conditions. The reasons behind the occurrence of these outliers at the specific 10 °C remain to be investigated.

Figure 9 Predictions of Δ¹⁴C in bulk soil vs. Δ¹⁴CO₂ with their equivalent simulation of mean age in bulk soil vs. mean transit time in CO₂ for parallel (a, b) for a slow cycling system. Series model presented exactly the same behaviour. Variation of k₂ with α=0.8 and γ=0.2. The brown box represents the range of measured Δ¹⁴C values obtained from incubated peatland soils for bulk (−90 ‰ to −65 ‰) and CO₂ (−18 ‰ to 25 ‰), excluding outliers.

Taken together, SOM stability and ecosystem characteristics defined the relation between Δ¹⁴C in bulk soil and Δ¹⁴CO₂, which in turn explained the difference in the age of heterotrophic respiration between grassland and peatlands. The Δ¹⁴CO₂ values obtained in our incubated peatlands indicate that such carbon was captured through photosynthesis during a time when the ¹⁴C levels of the atmosphere were depleted and have remained stored in the soil for a long enough time for radioactive decay to become relevant. However, it is imprecise to attribute a specific year to the entire bulk soil or CO₂ fractions, since they are formed by SOM accumulated in multiple steps and of different qualities and their respective decomposition.

In that sense, a Δ¹⁴C value is not per se indicative of age or transit time of SOC. To disentangle the utility of radiocarbon as a tool for SOM persistence in soil and be able to shed light on the timescales of carbon cycling in different ecosystems, we used SOC decomposition models. Our models contribute to understanding the time component of SOM persistence through the use of mean transit time as a parameter that integrates processes of SOM dynamics (Manzoni et al., 2009).

4.3 How can Δ¹⁴C values be interpreted to understand the effect of decomposition rates on mean age and mean transit time?

We simulated the dynamics of Δ¹⁴C as decomposition rates change and how this behaviour differs from those of mean age and mean transit time. Our simulations were able to reassemble the Δ¹⁴C ranges obtained from incubations and showed that modelled Δ¹⁴C values differed widely between slow cycling systems and fast cycling systems. Generally, decomposition rate, transfer rates and partitioning coefficients of a given model structure modulated Δ¹⁴C values and, in consequence, mean age and mean transit time (Bruun et al., 2004; Manzoni et al., 2009). For example, fast cycling systems with high decomposition rates resulted in a fast respiration represented by short mean transit times (lower than 20 years). At the same time, low decomposition rates resulted in longer mean transit times, which in turn represented older respired CO₂ (lower than 140 years). In slow cycling systems, where the decomposition rates were lower, transit times were consequently longer (lower than 20 000 years). Independently of the decomposition rates values, the series model structure increased the transit time of carbon compared to parallel structure due to the fact that a proportion of the carbon atoms had to pass through the two pools before being respired.

Other factors such as observation time and starting year of the simulation modified the response of Δ¹⁴C. This is due to the radioactive nature of the ¹⁴C isotope and the influence of the “bomb ¹⁴C”. For instance, the Δ¹⁴C enrichment during the last 60 years was clearly observed in fast cycling systems only with low k₂. Furthermore, the inflexion points on the Δ¹⁴C curves occurred at different values for slow and fast cycling systems. Hence, we could assume that a simulation with starting year after 1962 would not show an inflexion point in the Δ¹⁴C due to the absence of the bomb peak.

Our modelling results suggest that the increase of decomposition rates contributes to the release of older carbon in the respired flux as we hypothesized, but this depends on the initial state of the system. This interpretation is not straightforward since the stability of SOM depends on the specific combination of temperature and moisture for different ecosystems. For peatlands, only drier conditions and consequent increase of oxygen might increase decomposition rates and therefore cause the destabilization of SOM. In contrast, grasslands would need an increase of temperature to facilitate SOM decomposition, provided that moisture and oxygen are available (Sierra et al., 2017a; Azizi-Rad et al., 2022). SOM decomposition rates are expected to increase at higher temperatures (Leifeld and Fuhrer, 2005) and towards the extremes of the moisture range (Sierra et al., 2015). However, the destabilization of SOM can occur in any direction of WFPS variation (depending on the oxygen content) and in turn affect the age and transit time of carbon in a non-linear pattern.

From our results we could observe how the relation Δ¹⁴C bulk versus Δ¹⁴C–CO₂ properly represented the relation mean age versus mean transit time. We found that there is a good correspondence between both relations in the fast cycling systems as long as the decomposition rates (k) remain high. Such correspondence did not occur as k became smaller (typical for slow cycling systems) since the appearance of the bomb peak may have introduced anomalies that modified the equivalence between the two relations. In addition, the outlier Δ¹⁴C values found in our experimental data, where Δ¹⁴C changed drastically in the vertical and the horizontal direction, may be related to this complex response of the ¹⁴C tracer.

To quantify cycling times of carbon, radiocarbon can be used as a tool to understand SOM destabilization and persistence through the use of the concepts of age and transit time and their mutual relation. Nonetheless, it is essential to couple Δ¹⁴C measured values with a model that involves the dynamics of soil carbon in different pools and their interaction with the environment. Therefore, the acquisition of empirical data from soils (number of pools, I, C, k, γ, α) along with the correct setting of model structure will improve predictions on terrestrial and atmospheric carbon interactions.

5 Conclusions

Based on the incubation results of soils from the QTP, we showed that the Δ¹⁴C values of the peatland are significantly more depleted than the ones of the grassland both in the bulk soil and the heterotrophic respired CO₂. Our results indicated that changes in temperature did not systematically affect the radiocarbon values of heterotrophic respired CO₂ in any of the incubated soils despite increases in heterotrophic respired C at higher temperatures. Nevertheless, changes in WFPS had a relatively small effect on the ¹⁴CO₂ in grasslands but a significant influence on the ¹⁴CO₂ of peatlands, where higher WFPS levels led to more depleted ¹⁴CO₂ values except at 10 °C. In peatlands, more depleted ¹⁴C bulk values than ¹⁴C–CO₂ indicated that SOM stabilizes over time, and it is released once the stable conditions change. In grasslands, similar Δ¹⁴C values in bulk and heterotrophic respired CO₂ indicated that the soil behaved as a well-mixed homogeneous system due to either an absence of SOM stabilization or that the manipulation treatments were not long or strong enough to destabilize the existing old SOM. In this sense, the short duration of our incubations might have been an obstacle to register the influence of long-term factors such as climate seasonality and water table oscillation on SOM dynamics.

From our modelling approach, we conclude that radiocarbon can be used as a tool to understand SOM persistence through the use of the concepts of mean age and mean transit time and their mutual relation. Our simulations were able to reassemble the Δ¹⁴C values obtained from incubations and showed that modelled Δ¹⁴C values differed widely between slow cycling systems and fast cycling systems. We found that low values of decomposition rates, more common in slow cycling systems, modified the behaviour of Δ¹⁴C patterns due to the incorporation of ¹⁴C bomb in the soil system. Hence, the correspondence between these mutual relations strongly depended on the internal dynamics of the soil and its interaction with the environment. For this reason, the acquisition of empirical data from soils (number of pools, I, C, k, γ and α) along with the correct setting of model structure will improve our understanding on the stability of carbon in the soils of a changing QTP. In this way, current changes in climate patterns and land cover alteration may have a larger impact on the Zoigê peatlands than on the grasslands, given the vulnerability of large carbon stocks to be destabilized by changes in temperature. Nevertheless, the interaction with moisture may dampen or amplify the temperature effect, adding uncertainty on the future trajectories of soil carbon in the Qinghai–Tibetan Plateau.

Appendix A

Table A1 Mean daily CO₂ heterotrophic respiration for incubated soils under temperature and WFPS variation.

Table A2 Range of parameters used for simulations in a SOC decomposition model for fast and slow systems (Tangarife-Escobar et al., 2024). These ranges explore the required conditions to describe the Δ¹⁴C variation found in the incubated grassland and peatland soils and their equivalent age and transit time relationship. The target Δ¹⁴C intervals are shown as boxes in Figs. A1 to A6.

Δ¹⁴C variation in fast cycling systems

Figure A1 Predictions of Δ¹⁴C in bulk soil vs. Δ¹⁴CO₂ with their equivalent simulation of mean age in bulk soil vs. mean transit time in CO₂ for parallel (a, b) and series model structure (c, d) for a fast cycling system. k₁=0.8, k₂=0.1, with α=0–1 and γ=0–1. The green box represents the range of measured Δ¹⁴C values obtained from incubated grassland soils for bulk (20 ‰ to 74 ‰) and CO₂ (10 ‰ to 85 ‰), excluding outliers.

Figure A2 Predictions of Δ¹⁴C in bulk soil vs. Δ¹⁴CO₂ with their equivalent simulation of mean age in bulk soil vs. mean transit time in CO₂ for parallel (a, b) and series model structure (c, d) for a fast cycling system. Variation of k₁ with α=0.1 and γ=0.1. The green box represents the range of measured Δ¹⁴C values obtained from incubated grassland soils for bulk (20 ‰ to 74 ‰) and CO₂ (10 ‰ to 85 ‰), excluding outliers.

Figure A3 Predictions of Δ¹⁴C in bulk soil vs. Δ¹⁴CO₂ with their equivalent simulation of mean age in bulk soil vs. mean transit time in CO₂ for parallel (a, b) and series model structure (c, d) for a fast cycling system. Variation of k₂ with α=0.1 and γ=0.1. The green box represents the range of measured Δ¹⁴C values obtained from incubated grassland soils for bulk (20 ‰ to 74 ‰) and CO₂ (10 ‰ to 85 ‰), excluding outliers.

Δ¹⁴C variation in slow cycling systems

Figure A4 Predictions of Δ¹⁴C in bulk soil vs. Δ¹⁴CO₂ with their equivalent simulation of mean age in bulk soil vs. mean transit time in CO₂ for parallel (a, b) and series model structure (c, d) for a slow cycling system. Variation of k₁ with α=0.1 and γ=0.2. The brown box represents the range of measured Δ¹⁴C values obtained from incubated peatland soils for bulk (−90 ‰ to −65 ‰) and CO₂ (−18 ‰ to 25 ‰), excluding outliers.

Figure A5 Predictions of Δ¹⁴C in bulk soil vs. Δ¹⁴CO₂ with their equivalent simulation of mean age in bulk soil vs. mean transit time in CO₂ for parallel (a, b) and series model structure (c, d) for a slow cycling system. k₁ = 0.01, k₂=0.0001, with α=0–1 and γ=0–1. The brown box represents the range of measured Δ¹⁴C values obtained from incubated peatland soils for bulk (−90 ‰ to −65 ‰) and CO₂ (−18 ‰ to 25 ‰), excluding outliers.

Figure A6 Predictions of Δ¹⁴C in bulk soil vs. Δ¹⁴CO₂ with their equivalent simulation of mean age in bulk soil vs. mean transit time in CO₂ for parallel (a, b) and series model structure (c, d) for a slow cycling system. k₁=0.1, k₂=0.0001, with α=0–1 and γ=0–1. The brown box represents the range of measured Δ¹⁴C values obtained from incubated peatland soils for bulk (−90 ‰ to −65 ‰) and CO₂ (−18 ‰ to 25 ‰), excluding outliers.

Δ¹⁴C curves for bulk soil, respired CO₂ and atmosphere for a two-pool system in parallel and series structure

Figure A7 Variation of Δ¹⁴C values in atmosphere, bulk soil and respired CO₂ for a two-pool soil during the period 1950–2019. The similarity between Δ¹⁴C values may indicate both a well-mixed homogeneous system (type I) or a retention system (type II) where the bulk soil and the respired flux record the same value. Nonetheless, the mean transit time for the bulk and the soil respiration might be different due to the contrasting decomposition rates of the fast (k₁=0.8) and the slow pool (k₂=0.1). Simulation conducted with α=0.8 and γ=0.8.