Dynamics of deep soil carbon – insights from ¹⁴C time series across a climatic gradient

Quantitative constraints on soil organic matter (SOM) dynamics are essential for comprehensive understanding of the terrestrial carbon cycle. Deep soil carbon is of particular interest as it represents large stocks and its turnover times remain highly uncertain. In this study, SOM dynamics in both the top and deep soil across a climatic (average temperature ∼ 1–9 ^∘C) gradient are determined using time-series (∼20 years) ¹⁴C data from bulk soil and water-extractable organic carbon (WEOC). Analytical measurements reveal enrichment of bomb-derived radiocarbon in the deep soil layers on the bulk level during the last 2 decades. The WEOC pool is strongly enriched in bomb-derived carbon, indicating that it is a dynamic pool. Turnover time estimates of both the bulk and WEOC pool show that the latter cycles up to a magnitude faster than the former. The presence of bomb-derived carbon in the deep soil, as well as the rapidly turning WEOC pool across the climatic gradient, implies that there likely is a dynamic component of carbon in the deep soil. Precipitation and bedrock type appear to exert a stronger influence on soil C turnover time and stocks as compared to temperature.

1 Introduction

Within the broad societal challenges accompanying climate and land use change, a better understanding of the drivers of carbon turnover time in the largest terrestrial reservoir of organic carbon, as constituted by soil organic matter (SOM), is essential (Batjes, 1996; Davidson and Janssens, 2006; Doetterl et al., 2015; Prietzel et al., 2016). Terrestrial carbon turnover time remains one of the largest uncertainties in climate model predictions (Carvalhais et al., 2014; He et al., 2016). At present, there is no consensus on the net effect that climate and land use change will have on SOM stocks (Crowther et al., 2016; Gosheva et al., 2017; Melillo et al., 2002; Schimel et al., 2001; Trumbore and Czimczik, 2008). Deep soil carbon is of particular interest because of its large stocks (Jobbagy and Jackson, 2000; Balesdent et al., 2018; Rumpel and Kögel-Knabner, 2011) and perceived stability. The stability is indicated by low ¹⁴C content (Rethemeyer et al., 2005; Schrumpf et al., 2013; van der Voort et al., 2016) and low microbial activity (Fierer et al., 2003). Despite its importance, deep soil carbon has not been frequently studied and remains poorly understood (Angst et al., 2016; Mathieu et al., 2015; Rumpel and Kögel-Knabner, 2011). The inherent complexity of SOM and the multitude of drivers controlling its stability further impede the understanding of this globally significant carbon pool (Schmidt et al., 2011). In this framework, there is a particular interest in the portion of soil carbon that could be most vulnerable to change, especially in colder climates (Crowther et al., 2016). Water-exactable organic carbon (WEOC) is seen as a dynamic and potentially vulnerable carbon pool in the soil (Hagedorn et al., 2004; Lechleitner et al., 2016). Radiocarbon (¹⁴C) can be a powerful tool to determine the dynamics of carbon turnover time over decadal to millennial timescales because of the incorporation of bomb-derived ¹⁴C introduced in the atmosphere in the 1950s as well as the radioactive decay of ¹⁴C naturally present in the atmosphere (Torn et al., 2009). Furthermore, ¹⁴C can also be employed to identify petrogenic (or geogenic) carbon in the soil profile. Understanding the potential mobilization of stabilized petrogenic carbon is key because it could constitute an additional CO₂ source to the atmosphere (Hemingway et al., 2018). Time-series ¹⁴C data are particularly insightful because they enable the tracking of recent decadal carbon. Furthermore, single-time-point ¹⁴C data can yield two estimates for turnover time, whilst time-series data yield a single turnover time estimate (Torn et al., 2009). Given that the so-called “bomb radiocarbon spike” will continue to diminish in the coming decades, time-series measurements are increasingly a matter of urgency in order to take full advantage of this intrinsic tracer (Graven, 2015). Several case studies have collected time-series ¹⁴C soil datasets and demonstrated the value of this approach (Baisden and Parfitt, 2007; Prior et al., 2007; Fröberg et al., 2010; Mills et al., 2013; Schrumpf and Kaiser, 2015). However, these studies are rare, based on single specific sites, and have been rarely linked to abiotic and biotic parameters. Much more is yet to be learned about carbon cycling through time-series observations in top- and subsoils along environmental gradients. Furthermore, to our knowledge, there are no studies with pool-specific ¹⁴C soil time series focusing on labile carbon.

This study assesses two-pool soil carbon dynamics as determined by time-series (∼20 years) radiocarbon across a climatic gradient. The time-series data are analyzed by a numerically optimized model with robust error reduction to yield carbon turnover time estimates for the bulk and dynamic WEOC pools. Model output is linked to potential drivers such as climate, forest productivity, and physicochemical soil properties. The overall objective of this study is to improve our understanding of shallow and deep soil carbon dynamics in a wide range of ecosystems.

2 Materials and methods

2.1 Study sites, sampling strategy, and WEOC extraction

The five sites investigated in this study are located in Switzerland between 46–47^∘ N and 6–10^∘ E and encompass large climatic (mean annual temperature (MAT) 1.3–9.2 ^∘C, mean annual precipitation (MAP) 864–2126 mm m⁻² yr⁻¹), and geological gradients (Table 1). The sites are part of the Long-term Forest Ecosystem Research (LWF) program at the Swiss Federal Institute for Forest, Snow and Landscape Research, WSL (Schaub et al., 2011; Etzold et al., 2014). The soils of these sites were sampled between 1995 and 1998 (Walthert et al., 2002, 2003) and were resampled following the same sampling strategy in 2014 with the aim to minimize noise caused by small-scale soil heterogeneity. In both instances 16 samples were taken on a regular grid on the identical 43 m by 43 m (∼1600 m²) plot (Fig. 1; see Van der Voort et al., 2016, for further details). For the archived samples taken between 1995 and 1998, mineral soil samples down to 40 cm depth (intervals of 0–5, 5–10, 10–20 and 20–40 cm) were taken on an area of 0.5 m by 0.5 m (0.25 m²). For samples >40 cm (intervals of 40–60, 60–80, and 80–100 cm), corers were used to acquire samples (n=5 in every pit, area $\sim \mathrm{2.8}\times {\mathrm{10}}^{-\mathrm{3}}$ m²). The organic layer was sampled by use of a metal frame (30 cm×30 cm). The samples were dried at 35–40 ^∘C, sieved to remove coarse material (2 mm), and stored in hard plastic containers under controlled climate conditions in the “Pedothek” at WSL (Walthert et al., 2002). For the samples acquired in 2014 the same sampling strategy was followed and samples were taken on the exact same plot proximal (∼10 m) to the legacy samples. For the sampling, a SHK Martin Burch AG HUMAX soil corer ($\sim \mathrm{2}\times {\mathrm{10}}^{-\mathrm{3}}$ m²) was used for all depths (0–100 cm). For the organic layer, a metal frame of 20 cm×20 cm was used to sample. Samples were sieved (2 mm), frozen, and freeze-dried using an oil-free vacuum-pump-powered freeze dryer (Christ, Alpha 1–4 LO plus). For the time-series radiocarbon measurements, all samples covering ∼1600 m² were pooled to one composite sample per soil depth using the bulk density. In order to determine bulk density of the fine earth of the 2014 samples, stones >2 mm were assumed to have a density of 2.65 g cm⁻³. For the Alptal site, 16 cores were taken on a slightly smaller area (∼1500 m²), which encompasses the control plot of a nitrogen addition experiment (NITREX project) (Schleppi et al., 1998). For this site, no archived samples are available and thus only the 2014 samples were analyzed. Soil carbon stocks were estimated by multiplying soil organic carbon (SOC) concentrations with the mass of soil calculated from measured bulk densities and stone contents for each depth interval (Gosheva et al., 2017). For the Nationalpark site, the soil carbon stocks from 80 to 100 cm were estimated using data from a separately dug soil profile (Walthert et al., 2003) because the HUMAX corer could not penetrate the rock-dense soil below 80 cm depth. In order to understand very deep soil carbon dynamics (i.e., >100 cm), this study also includes single-time-point ¹⁴C analyses of soil profiles that were dug down to the bedrock between 1995 and 1998 as part of the LWF program on the same sites (Walthert et al., 2002). The sampling of the profiles has not yet been repeated.

Figure 1 Sample locations, all of which are part of the Long-term ecosystem research program (LWF) of the Swiss Federal Institute for Forest, Snow and Landscape Research, WSL: (1) Othmarsingen, (2) Lausanne, (3) Alptal, (4) Beatenberg, and (5) Nationalpark. Image made using 2016 swisstopo (JD100042).

Table 1 Overview of sampling locations and climatic and ecological parameters.

^a Walthert et al. (2003). ^b Etzold et al. (2014). ^c Von Arx et al. (2013). ^d Krause et al. (2013) for Alptal data.

2.2 Climate and soil data

Temperature and precipitation data are derived from weather stations close to the study sites that have been measuring for over 2 decades, yielding representative estimates of both variables and over the time period concerned in this study (Etzold et al., 2014). The pH values for all sites and concerned depth intervals were acquired during the initial sampling campaign (Walthert et al., 2002). At Alptal, pH values were determined as described in Xu et al. (2009), values of 10–15 cm were extrapolated to the deeper horizons because of the uniform nature of the gley horizon. The Beatenberg Podzol is marked by strong eluviation (∼ 4–35 cm) and illuviation (∼ 35–60 cm) (Walthert et al., 2003). Exchangeable cations were extracted (in triplicate) from the 2 mm sieved soil in an unbuffered solution of 1 M NH₄Cl for 1 h on an end-over-end shaker using a soil-to-extract ratio of 1:10. The element concentrations in the extracts were determined by inductively coupled plasma atomic emission spectroscopy (ICP-AES) (Optima 3000, PerkinElmer). Contents of exchangeable protons were calculated as the difference between the total and the Al-induced exchangeable acidity as determined (in duplicate) by the KCl method (Thomas, 1982). This method was applied only to soil samples with a pH (CaCl₂) < 6.5. In samples with a higher pH, we assumed the quantities of exchangeable protons were negligible. The effective cation-exchange capacity (CEC) was calculated by summing up the charge equivalents of exchangeable Na, K, Mg, Ca, Mn, Al, Fe, and H. The base saturation (BS) was defined as the percental fraction of exchangeable Na, K, Mg, and Ca of the CEC (Walthert et al., 2002, 2013). Net primary production (NPP) was determined by Etzold et al. (2014) as the sum of carbon fluxes by woody tree growth, foliage, fruit production, and fine-root production. Soil texture (sand, silt, and clay content) on plot-averaged samples taken in 2014 have been determined using grain size classes for sand, silt, and clay of 0.05–2, 0.002–0.05, and <0.002 mm, respectively, according to Klute (1986). The continuous distribution of grain sizes was also determined after removal of organic matter (350 ^∘C for 12 h) using the Mastersizer 2000 (Malvern Instruments Ltd.). Soil water potential (SWP) was measured on the same sites as described in Von Arx et al. (2013). In accordance with Mathieu et al. (2015), topsoil refers to the mineral soil down to 20 cm depth, and deep soil refers to mineral soil below 20 cm. Out of the five sites, two are hydromorphic (Gleysol and Podzol in Alptal and Beatenberg, respectively), whilst the others are nonhydromorphic (Luvisol, Cambisol, and Fluvisol in Othmarsingen, Lausanne, and Nationalpark, respectively).

2.3 Isotopic (¹⁴C, ¹³C) and compositional (C, N) analyses

Prior to the isotopic analyses, inorganic carbon in all samples was removed by vapor acidification for 72 h (12 M HCl) in desiccators at 60 ^∘C (Komada et al., 2008). After fumigation, the acid was neutralized by substituting NaOH pellets for another 48 h. All glassware used during sample preparation was cleaned and combusted at 450 ^∘C for 6 h prior to use. Water-extractable organic carbon (WEOC) was procured by extracting dried soil with 0.5 wt % precombusted NaCl in ultrapure Milli-Q (MQ) water in a 1:4 soil : water mass ratio (adapted from Hagedorn et al., 2004, details in Lechleitner et al., 2016).

In order to determine absolute organic carbon and nitrogen contents, as well as ¹³C values, an elemental analyzer and isotope ratio mass spectrometer system was used (EA-IRMS, Elementar, vario MICRO cube – Isoprime, visION). Atropine (Säntis) and an in-house standard peptone (Sigma) were used for the calibration of the EA-IRMS for carbon concentration, nitrogen concentration and C : N ratios, and ¹³C, respectively. High ¹³C values were used to flag if all inorganic carbon had been removed by acidification.

The ¹⁴C measurements of the bulk soil samples were first graphitized using an EA-AGE (elemental analyzer and automated graphitization equipment, Ionplus AG) system at the Laboratory of Ion Beam Physics at ETH Zürich (Wacker et al., 2009). Graphite samples were measured on a MICADAS (MIniturised radioCArbon DAting System, Ionplus AG) also at the Laboratory of Ion Beam Physics, ETH Zürich (Wacker et al., 2010). For three samples (Alptal depth intervals 40–60, 60–80, and 80–100 cm) the ¹⁴C signature was directly measured as CO₂ gas using the recently developed online elemental analyzer (EA) – stable isotope ratio mass spectrometers (IRMS)–AMS system at ETH Zürich (McIntyre et al., 2016). Oxalic acid (NIST SRM 4990C) was used as the normalizing standard. Phthalic anhydride and in-house anthracite coal were used as blank. Two in-house soil standards (Alptal soil 0–5 cm, Othmarsingen soil 0–5 cm) were used as secondary standards. For the WEOC, samples were converted to CO₂ by wet chemical oxidation (WCO) (Lang et al., 2016) and run on the AMS using a gas ion source (GIS) interface (Ionplus). To correct for contamination, a range of modern standards (sucrose, Sigma, ${\mathit{\delta }}^{\mathrm{13}}\mathrm{C}=-\mathrm{12.4}$ ‰ VPDB, $F{}^{\mathrm{14}}\mathrm{C}=\mathrm{1.053}\pm \mathrm{0.003}$) and fossil standards (phthalic acid, Sigma, ${\mathit{\delta }}^{\mathrm{13}}\mathrm{C}=-\mathrm{33.6}$ ‰ VPDB, F¹⁴C<0.0025) were used (Lechleitner et al., 2016).

2.4 Numerical optimization to find carbon turnover time and size of the dynamic pool

2.4.1 Turnover time based on a single ¹⁴C measurement

The ¹⁴C signature of a sample can be used to estimate turnover time of a carbon pool (Torn et al., 2009).

$$\begin{array}{}\text{(1)}& {\displaystyle }{R}_{\mathrm{sample},t}=k\times R_{\mathrm{atm},t}+\left(\mathrm{1}-k-\mathit{\lambda }\right)\times R_{\mathrm{sample}\left(t-\mathrm{1}\right)}\text{(2)}& {\displaystyle }{R}_{\mathrm{sample},t}={\displaystyle \frac{\mathrm{\Delta }{}^{\mathrm{14}}{\mathrm{C}}_{\mathrm{sample}}}{\mathrm{1000}}}+\mathrm{1}\end{array}$$

In Eqs. (1)–(2), the constant for radioactive decay of ¹⁴C is indicated as λ; the decomposition rate k (inverse of turnover time) is the only unknown in this equation and is hence the variable for which the optimal value that fits the data is sought using the model. The R value of the sample is inferred from Δ¹⁴C, hence accounting for the sampling year, as shown in Eq. (2) (Herold et al., 2014; Solly et al., 2013). In order to avoid ambiguity, the term turnover time – and not, for example, mean residence time – is used solely in this article (Sierra et al., 2016).

For the turnover time estimation, we assumed the system to be in steady state over the modeled period ($\sim \mathrm{1}\times {\mathrm{10}}^{\mathrm{4}}$ years, indicating soil formation since the last glacial retreat; Ivy-Ochs et al., 2009), hence accounting both for radioactive decay and incorporation of the bomb-testing-derived material produced in the 1950s and 1960s (Eq. 1) (Herold et al., 2014; Torn et al., 2009). We assumed an initial fraction of the modern (F_m) ¹⁴C value of 1 at 10 000 BCE. For the period after 1900, atmospheric fraction modern (F_m) values of the Northern Hemisphere were used (Hua et al., 2013). This equation could be solved in Microsoft Excel with manual iterations (e.g., Herold et al., 2014), or alternatively a numerical optimization can be used to find the best fit automatically. In this paper, we used a numerical optimization constructed in MATLAB version 2015a (MathWorks, Inc., Natick, Massachusetts, US) to find the best fit. The numerical optimization is exhaustive, meaning that every single turnover time value from 1 to 10 000 years with a 0.1-year interval is tested. The error is defined as the difference between the fitted value of R and the measured value (Eq. 3). The turnover time value with the lowest error is then automatically selected.

$$\begin{array}{}\text{(3)}& {\displaystyle }{\mathrm{Error}}_{\mathrm{single}\mathrm{timepoint}}=\mathrm{|}{R}_{\mathrm{calculated}}-R_{\mathrm{measured}}\mathrm{|}\end{array}$$

The residual errors of each fit are provided in Table S3 in the Supplement. Turnover times determined with the numerical optimization match the manually optimized turnover time modeling published previously (Herold et al., 2014; Solly et al., 2013).

2.4.2 Turnover time based on two ¹⁴C measurements

A single ¹⁴C value could yield many possible turnover time values (Torn et al., 2009; Graven, 2015). If there is a time-series ¹⁴C dataset, this problem can be eliminated. In this paper, we have time-series data of both the bulk soil and the vulnerable fraction (WEOC). For all samples a time-series dataset is available, both data points are employed to give the best estimate of turnover time. The same numerical optimization (Eqs. 1 and 2) as we did for a single time point, except that we try to find the best fit for both time points whilst reducing the compounding residual mean square error (RMSE, Eq. 4). As can be seen in Fig. 2a, single time points can yield two likely turnover times. But when two datapoints are available, a single value can be found. The input data for Fig. 2 can be found in Table S1. The results of the time-series turnover time modeling for both the bulk and WEOC pools of the subalpine site Beatenberg are shown in Fig. 3.

$$\begin{array}{}\text{(4)}& \begin{array}{rl}& {\mathrm{Error}}_{\mathrm{two}\mathrm{time}\mathrm{points}}=\\ & \sqrt[\mathrm{2}]{\mathrm{|}{R}_{\mathrm{calculated}}-R_{\mathrm{measured}}{\mathrm{|}}_{\mathrm{time}\mathrm{point}\mathrm{1}}^{\mathrm{2}}}\\ & \stackrel{\mathrm{‾}}{+\mathrm{|}{R}_{\mathrm{calculated}}-R_{\mathrm{measured}}{\mathrm{|}}_{\mathrm{time}\mathrm{point}\mathrm{2}}^{\mathrm{2}}}\end{array}\end{array}$$

Figure 2 Numerical optimization of least mean-square error reduction, showing the reduction of error spread for two soil depths. For the Beatenberg Podzol organic layer (a) the individual ¹⁴C time points for both 1997 and 2014 both yield two solutions and are almost equally likely (i.e., the error nears zero). The combined optimization using both time points reveals the likeliest option. For the (b) 0–5 cm layer, the single time points only have a single likely solution.

Figure 3 (a) Time series of soil carbon turnover time in years as determined by numerical modeling for (b) subalpine site Podzol Beatenberg. The bulk turnover time in the organic layer is rapid (14 years), followed by the turnover time of the water-extractable organic carbon (WEOC) (191 years) and the bulk turnover time of the soil (410 years) at 0–5 cm depth. Photo soil profile courtesy of Marco Walser, WSL.

2.4.3 Vegetation-induced lag

In order to account for vegetation lag, two scenarios were run: firstly (1) with no assumed lag between the fixation of carbon from the atmosphere and input into to the soil and (2) model run with a lag of fixation of the atmospheric carbon as inferred from the dominant vegetation (Von Arx et al., 2013; Etzold et al., 2014). In the case of full deciduous trees coverage a lag of 2 years was assumed, and for the case of 100 % conifer-dominated coverage a lag of 8 years was incorporated (Table 1).

2.4.4 Turnover time and size of vulnerable pool based on the two-pool model

As SOM is complex and composed of a continuum of pools with various ages (Schrumpf and Kaiser, 2015) and there are data available from two SOM pools, the ¹⁴C time-series data can be leveraged to create a two-pool model. The following assumptions were made: First, both pools (slow and fast) make up the total carbon pool (Eq. 5). Secondly, the total turnover time of the bulk soil is made up of the “dynamic” fraction turnover time multiplied by dynamic fraction pool size and the “slow” pool turnover time multiplied by slow pool size (Eq. 6). Furthermore, we assume that the signature of the sample (the time-series bulk data) is determined by the rate of incorporation of the material (atmospheric signal) and the loss of carbon the two pools (Eq. 7). Lastly, we assume that the radiocarbon signal of the WEOC pool is representative for a dynamic pool, as it could be representative for a larger component of rapidly turning-over carbon, even in the deep soil (Baisden and Parfitt, 2007; Koarashi et al., 2012). The turnover time of the slow pool was set between 100 and 10 000 years, with a time step of 10 years. The size of the dynamic pool was set to be between 0 and 0.5, with a size step of 0.01.

$$\begin{array}{}\text{(5)}& {\displaystyle }\mathrm{1}=F_{\mathrm{1}}+F_{\mathrm{2}}\text{(6)}& {\displaystyle }{k}_{\mathrm{total}}=(F_{\mathrm{1}}/k_{\mathrm{1}}+F_{\mathrm{2}}/k_{\mathrm{2}}{)}^{-\mathrm{1}}\text{(7)}& {\displaystyle }\begin{array}{rl}{R}_{\mathrm{sample},t}& =k_{\mathrm{total}}\times R_{\mathrm{atm},t}+F_{\mathrm{1}}[\left(\mathrm{1}-k_{\mathrm{1}}-\mathit{\lambda }\right)\\ & \times R_{\mathrm{sample}\left(t-\mathrm{1}\right)}]+F_{\mathrm{2}}\left(\mathrm{1}-k_{\mathrm{2}}-\mathit{\lambda }\right)\\ & \times R_{\mathrm{sample}\left(t-\mathrm{1}\right)},\end{array}\end{array}$$

where F₁ is the relative size of the dynamic pool and F₂ is the relative size of the (more) stable pool. The k₁ is the inverse of the turnover time of the dynamic or WEOC as determined using the numerical optimization of Eqs. (1)–(4). The k₂ is the inverse of the turnover time of the slow pool. The calculation of the error term becomes complex because it needs to be recalculated for each unique combination of pool-size distribution (Eq. 5) and turnover time (inverse of k, Eq. 6). Therefore, the error space changes from column vector to a two-dimensional matrix of length of the step size increments (F₁) and width of the inverse of the turnover time of the slow pool (k₂).

$$\begin{array}{}\text{(8)}& {\displaystyle }\begin{array}{rl}& {\mathrm{Error}}_{k_{\mathrm{2}},F_{\mathrm{1}}}=\\ & \sqrt[\mathrm{2}]{\mathrm{|}{R}_{\mathrm{calculated}}-R_{\mathrm{measured}}{\mathrm{|}}_{\mathrm{time}\mathrm{point}\mathrm{1}}^{\mathrm{2}}}\\ & \stackrel{\mathrm{‾}}{+\mathrm{|}{R}_{\mathrm{calculated}}-R_{\mathrm{measured}}{\mathrm{|}}_{\mathrm{time}\mathrm{point}\mathrm{2}}^{\mathrm{2}}}\end{array}\text{(9)}& {\displaystyle }\mathrm{Error}=\mathrm{Min}\left({\mathrm{Error}}_{k_{\mathrm{2}},F_{\mathrm{1}}}\right)\end{array}$$

The numerical optimization finds the likeliest solution for the given dataset. This model constitutes a best fit and more data would better constrain the results. Additional details can be found in the Supplement and Fig. S1. All MATLAB-based numerical optimization codes can be found in the Supplement. For correlations (packages HMISC and corrgram, method is Pearson), statistical software R version 1.0.153 was used.

3 Results

3.1 Changes in radiocarbon signatures over time

Overall, there is a pronounced decrease in radiocarbon signature with soil depth at all sites (Fig. 4). The time-series results show clear changes in radiocarbon signature over time from the initial sampling period (1995–1998) as compared to 2014, with the magnitude of change depending on site and soil depth. In the uppermost 5 cm of soils, the overarching trend in the bulk soil is a decrease in the ¹⁴C bomb-spike signature in the warmer climates (Othmarsingen, Lausanne), whilst at higher-elevation (colder) sites (Beatenberg, Nationalpark) the bomb-derived carbon appears to enter the topsoil between 1995–1998 and 2014.

Figure 4 (a–e) Changes in radiocarbon signature of both bulk soil and WEOC over 2 decades at four sites on a climatic gradient. For Alptal (Gleysol) (c) only the 2014 time point was available. For the warmer locations Othmarsingen and Lausanne (Luvisol, Cambisol MAT 9.2–7.6 ^∘C), depletion in bomb-derived radiocarbon occurs in the first 5 cm of soil in 2014 as compared to 1995–1998. The colder Beatenberg site (Podzol, MAT 4.7 ^∘C) is marked by a clear enrichment of ¹⁴C in the mineral soil in 2014 w.r.t. 1997. At the coldest site Nationalpark (Fluvisol, MAT 1.3 ^∘C) almost all samples, taken 2 decades after the initial sampling, show an enrichment in radiocarbon signature. WEOC contains bomb-derived carbon in the topsoil in 2014 at all sites.

Water-extractable organic carbon (WEOC) has an atmospheric ¹⁴C signature in the topsoil at all sites in 2014. The deep soil in the 1990s still has a negative Δ¹⁴C signature of WEOC at multiple sites. There are two distinguishable types of depth trends for WEOC in the 2014 dataset: (1) WEOC has the same approximate ¹⁴C signature throughout depth (Othmarsingen, Beatenberg) and (2) WEOC becomes increasingly ¹⁴C depleted with depth (Alptal, Nationalpark), or an intermediate form where WEO¹⁴C is modern throughout the topsoil but becomes more ¹⁴C depleted in the deep soil (Lausanne) (Fig. 4). The isotopic trends of WEOC covary with grain size as inherited from the bedrock type (Walthert et al., 2003). Soils with a relatively modern WEO¹⁴C signature in 2014 (down to 40 cm) are underlain by bedrock with large-grained (Fig. S2, Table S3) components (the moraines and sandstone at Othmarsingen, Lausanne, and Beatenberg). Soils where WEO¹⁴C signature decreases with depth are underlain by bedrock containing fine-grained components, e.g. the flysch in Alptal (Schleppi et al., 1998) and intercalating layers of silt and coarse-grained alluvial fan in Nationalpark (Walthert et al., 2003).

3.2 Carbon turnover time patterns

Incorporation of a vegetation-induced time lag (Tables 2 and S2) has an effect on modeled carbon dynamics in the organic layer, but this effect is strongly attenuated in the 0–5 cm layer in the mineral soil and virtually absent for the deeper soil layers. The residual errors associated with the carbon turnover time estimates converge to a single point (Fig. 2) and are low (i.e., <0.06R, Tables S3 and S4). Turnover times show two modes of behavior for well-drained soils and hydromorphic soils. The nonhydromorphic soils have relatively similar values with decadal turnover times for the 0–5 cm layer, increasing to an order of centuries down to 20 cm depth, and to millennia in deeper soil layers (∼980 to ∼3940 years at 0.6 to 1 m depth) (Fig. 5). In contrast, the hydromorphic soils are marked by turnover times that are up to an order of magnitude larger, from centennial in topsoil to (multi)millennial in deeper soils. At the Beatenberg Podzol, turnover time in the shallow layer which overlaps with the eluvial horizon is slower (20–40 cm, ∼1900 years) than the deepest layer (40–60 cm, ∼1300 years) which overlaps with the illuvial horizon (Fig. 5, Table S5).

Figure 5 Carbon (a) stocks in the mineral soil (kg C m⁻²), (b) turnover time bulk soil in years, and (c) turnover time water-extractable organic carbon soil in years. Locations are ordered from the warmest to coldest sites, i.e., (1) Othmarsingen (Luvisol), (2) Lausanne (Cambisol), (3) Alptal (Gleysol), (4) Beatenberg (Podzol), and (5) Nationalpark (Fluvisol). Gray boxes indicate the absence of material; black boxes indicate the occurrence of the C horizon (poorly consolidated bedrock-derived stony material or bedrock itself).

Table 2 Vegetation and soil data of the study sites. Soil water potentials (hPa) are for 15 cm depth.

^a Walthert et al. (2003). ^b Von Arx et al. (2013). ^c Krause et al. (2013).

Carbon stocks also show distinct differences between drained and hydromorphic soils with greater stock in the hydromorphic soils (∼15 kg C m⁻² at Beatenberg and Alptal vs. ∼6–∼7 kg C m⁻² at Othmarsingen, Lausanne, and Nationalpark; Fig. 5, Table 3).

Table 3 Soil properties as well as carbon stocks and fluxes in 0–20, 20–60, and 60–100 cm depths of the study sites for the bulk and water-extractable organic carbon (WEOC).

^a Walthert et al. (2002, 2003). Fe and Al content (mmol_c kg⁻¹) determined by NH₄Cl extraction. For the 0.2–0.6 depth interval, the CEC determined for 0.2–0.4 m was taken, and similarly for the
depth interval 0.6–1.0 m the values for 0.6–0.8 m were taken in the case of Othmarsingen, Lausanne, Beatenberg, and Nationalpark. ^b Krause et al. (2013). ^c Diserens et al. (1992). CEC determined
(mmol_c kg⁻¹) hydrogen, lead, and zinc ions were not included; aluminum content determined by the Lakanen method. CEC values for 0.2–0.4 m were extrapolated to 1 m. ^d Xu et al. (2009).
^e Depth to 0.8 m. ^f Depth to 0.4 m.

The turnover times of the WEOC mimic the trends in the bulk soil but are up to an order of magnitude faster. Considering WEOC turnover time in the nonhydromorphic soils only, there is a slight increase in WEOC turnover time with decreasing site temperature, but the trend is not significant (Table S4). The modeled estimate for dynamic fraction is variable at the surface but decreases towards the lower topsoil (from ∼0.2 at 0–5 cm to ∼0.01 at 10–20 cm in Othmarsingen). In the deep soil, the model indicates there could also be a non-negligible proportion of dynamic carbon (e.g., 0.10–0.23 at 20–40 cm). The residual errors associated with the error reduction of the two-pool model are also low (i.e., <0.06R), but do not converge as strongly as the single-pool model (Fig. S1).

3.3 Preglacial carbon in deep soil profiles

The turnover times of deep soil carbon exceed 10 000 years in several profiles, indicating the presence of carbon that predates the glacial retreat (Fig. 6). These profiles are located on carbon-containing bedrock and concern the deeper soil (80–100 cm) of the Gleysol (Alptal), as well as >100 cm in the Cambisol (Lausanne) (Fig. 6, Table S6).

Figure 6 Modeled turnover times (yrs) of single profiles sampled down to the bedrock between 1995 and 1998; Δ¹⁴C data published in Van der Voort et al. (2016). Results indicate the presence of petrogenic (bedrock-derived) carbon as modeled turnover time exceeds soil formation since the end of the last ice age (10 000 years) in Lausanne (>100 cm, Cambisol) and Alptal (80–100 cm, Gleysol). For Beatenberg (Podzol) and Nationalpark (Fluvisol), no petrogenic carbon was found.

3.4 Environmental drivers of carbon dynamics

Pearson correlation was used to assess potential relationships between carbon stocks and turnover time and their potential controlling factors (climate, NPP, soil texture, soil moisture, and physicochemical properties; Tables 4, S7, and S8). For the averaged topsoil (0–20 cm, n=5), carbon stocks were significantly positively correlated to mean annual precipitation (MAP). Turnover time in the bulk topsoil negatively correlated with silt content and positively with average grain size. Turnover time in the WEOC of the topsoil did not correlate significantly with any parameter except a weak positive correlation with grain size. Deeper soil bulk stock and turnover time positively correlated with MAP and iron content.

Table 4 Pearson correlations for averaged depth intervals for the topsoil (0–20 cm, n=5) and deep soil (20–60 cm, n=5). Significance is denoted with ^•, ^*, or ${}^{**}$ for p values smaller than 0.1 (marginally significant) 0.05, and 0.005 (significant), respectively. Nonsignificant correlations are indicated by the superscript “ns”. SWP or soil water potential used are the median values at 15 cm for each of these five sites (Von Arx et al., 2013). Water-extractable carbon is abbreviated to WEOC. Results indicate that no single climatic or textural factor consistently covaries with carbon stocks, or turnover time.

4 Discussion

4.1 Dynamic deep soil carbon

4.1.1 Rapid shifts in ¹⁴C abundance reflect dynamic deep carbon

The propagation of bomb-derived carbon into supposedly stable deep soil on the bulk level across the climatic gradient implies that SOM in deep soil contains a dynamic pool and could be less stable and potentially more vulnerable to change than previously thought. This possibility is further supported by the WEO¹⁴C, which is consistently more enriched in bomb-derived carbon than the bulk soil. Near-atmospheric WEO¹⁴C signature pervades down to 40 or even 60 cm depth. Hagedorn et al. (2004) also found WEOC to be a highly dynamic pool using ¹³C tracer experiments in forest soils.

We consider our ¹⁴C comparison over time to be robust because the grid-based sampling and averaging were repeated on the same plots, which excludes the effect of plot-scale variability (Van der Voort et al., 2016). Our ¹⁴C time-series data in the deep soil corroborate pronounced changes in ¹⁴C (hence substantial SOM changes) observed in subsoils of an area with pine afforestation (Richter and Markewitz, 2001). The findings are also in agreement with results from an incubation study by Fontaine et al. (2007) which showed that the deep soil can have a significant dynamic component. Baisden and Parfitt (2007) also found indications of a deep dynamic pool using modeling of ¹⁴C time series on the bulk level for a New Zealand soil under stable pastoral management.

4.1.2 Carbon dynamics reflect soil-specific characteristics at depth

Bulk carbon turnover time for the top and deeper soil fall in the range of prior observations and models, although the data for the latter category are rare (Scharpenseel and Becker-Heidelmann, 1989; Paul et al., 1997; Schmidt et al., 2011; Mills et al., 2013; Braakhekke et al., 2014). The carbon turnover time is related to soil-specific characteristics. The higher turnover times of hydromorphic as compared to nonhydromorphic soils are likely due to increased waterlogging and limited aerobicity (Hagedorn et al., 2001b), which is conducive to slow decomposition rates and enhanced carbon accumulation. The WEOC turns over up to an order of magnitude faster than the bulk and mirrors these trends, indicating that it indeed is a more dynamic pool (Hagedorn et al., 2004; Lechleitner et al., 2016). Results also reflect known horizon-specific dynamics for certain soil types, particularly in the deep soil. The hydromorphic Podzol at Beatenberg shows specific pedogenetic features such as an illuviation layer with an enrichment in humus and iron in the deeper soil (Walthert et al., 2003) where the turnover times of bulk and WEOC are more rapid and stocks are higher as compared to the eluvial layer above (Fig. 5). This is likely due to the input of younger carbon via leaching of dissolved organic carbon. The nonhydromorphic Luvisols are marked by an enrichment of clay in the deeper soil, which can enhance carbon stabilization (Lutzow et al., 2006). This is also reflected in the turnover time of the 60–80 cm layer in the Othmarsingen Luvisol; in this clay-enriched depth interval (Walthert et al., 2003), turnover time is relatively high as compared to the other (colder) nonhydromorphic soils (Fig. 5). These patterns are consistent with findings by Mathieu et al. (2015) that the important role of soil pedology on deep soil carbon dynamics.

4.1.3 Dynamic carbon at depth and implications for carbon transport

The analytical ¹⁴C data as well as turnover time estimates indicate that there is likely a dynamic portion of carbon in the deep soil. The estimated size of the dynamic pool can be large, even at greater depth than it was observed by other ¹⁴C time series (Richter and Markewitz, 2001; Baisden and Parfitt, 2007; Koarashi et al., 2012). The two-pool modeling indicates that the size of the dynamic pool in the deep soil can be upwards of ∼10 %. A deep dynamic pool is consistent with findings of a ¹³C tracer experiment by Hagedorn et al. (2001) that shows that relatively young (<4 years) carbon can be rapidly incorporated into the topsoil (20 % new C at 0–20 cm depth) but also in the deep soil (50 cm), and findings by Balesdent et al. (2018) which estimate that up to 21 % of the carbon between 30 and 100 cm is younger than 50 years. Rumpel and Kögel-Knabner (2011) have highlighted the importance of the poorly understood deep soil carbon stocks and a significant dynamic pool in the deep soil could imply that carbon is more vulnerable than initially suspected. One major input pathway of younger C into deeper soils is the leaching of DOC (Kaiser and Kalbitz, 2012; Sanderman and Amundson, 2009). Here, we have measured WEOC – likely primarily composed of microbial metabolites (Hagedorn et al., 2004) – carrying a younger ¹⁴C signature than bulk SOM, and thus representing a translocation of fresh carbon to the deep soil. The WEOC turnover time is in the order of decades, implying that it is not directly derived from decaying vegetation but rather composed of microbial material feeding on the labile portion of the bulk soil. In addition to WEOC, roots and associated mycorrhizal communities may also provide a substantial input of new C into soils in deeper soils (Rasse et al., 2005). Additional modeling such as in CENTURY and RotC could provide additional insights into the soil carbon dynamics and fluxes (Manzoni et al., 2009).

4.2 Contribution of petrogenic carbon

Our results on deep soil carbon suggest the presence of pre-aged, or ¹⁴C-dead (fossil), pre-interglacial carbon in the Alptal (Gleysol) and Lausanne (Cambisol) profiles, implying that a component of soil carbon is not necessarily linked to recent (< millennial) terrestrial productivity and instead constitutes part of the long-term (geological) carbon cycle (> millions of years). In the case of the Gleysol in Alptal, the ¹⁴C-depleted material could be derived from the poorly consolidated sedimentary rocks (flysch) in the region (Hagedorn et al., 2001a; Schleppi et al., 1998; Smith et al., 2013), whereas carbon present in glacial deposits and molasse may contribute to deeper soils at the Lausanne (Cambisol) site. The potential contribution of fossil carbon was estimated using a mixing model using the signature of a soil without fossil carbon, the signature of fossil carbon and the measured values (Table S4). Fossil carbon contribution in the Alptal profile between 80 and 100 cm (Fig. 6, Table S4) is estimated at ∼40 %. Below 1 m at the Lausanne site, the petrogenic percentages range from ∼20 % at 145 cm up to ∼80 % at 310 cm depth (Fig. 6, Table S4).

Other studies analyzing soils have observed the significant presence of petrogenic soil (geogenic in soil science terminology) in loess-based soils (Helfrich et al., 2007; Paul et al., 2001). Our results suggest that preglacial carbon may comprise a dominant component of deep soil organic matter in several cases, resulting in an apparent increase in the average turnover time of carbon in these soils. Hemingway (2018) highlighted that fossil carbon oxidized in soils can lead to significant additional CO₂ emissions. Therefore, the potential of soils to activate fossil petrogenic carbon should be considered when evaluating the soil carbon sequestration potential.

4.3 Controls on carbon dynamics and cycling

In order to examine the effects of potential controls on soil C turnover time and stocks, we explore correlations between a number of available factors which have previously been proposed, such as texture, geology, precipitation, temperature, and soil moisture (Doetterl et al., 2015; McFarlane et al., 2013; Nussbaum et al., 2014; Seneviratne et al., 2010; van der Voort et al., 2016).

From examination of data for all samples, it emerges that C turnover time does not exhibit a consistent correlation with any specific climatological or physicochemical factor. This implies that no single mechanism predominates and/or that there is a combined impact of geology and precipitation as these soil-forming factors affect grain size distribution, water regime, and mass transport in soils. Exploring potential relationships in greater detail, we see that carbon stocks in the topsoil and deep soil as well as turnover time are positively related to MAP, which could be linked to waterlogging and anaerobic conditions even in upland soils leading to a lower decomposition and thus to a higher build-up of organic material (Keiluweit et al., 2015). Our results are supported by the findings, based on >1000 forest sites, that precipitation exerts a strong effect on soil C stocks across Switzerland (Gosheva et al., 2017; Nussbaum et al., 2014). Furthermore, Balesdent et al. (2018) also highlighted the role of precipitation and evapotranspiration on deep soil organic carbon stabilization. Nonetheless, it has to be noted that for these sites, the precipitation range does not include very dry soils (MAP 864–2126 mm yr⁻¹). Topsoil carbon turnover time in bulk and WEOC positively correlates with grain size, which is likely caused by the large grain size of the waterlogged Podzol Beatenberg. Overall, geology seems to impact the carbon cycling in three key ways. Firstly, when petrogenic carbon is present in the bedrock from shale or reworked shale (Schleppi et al., 1998; Walthert et al., 2003), fossil carbon contributes to soil carbon. Secondly, porosity of underlying bedrock either prevents or induces waterlogging, which in turn affects turnover time. Thirdly, the initial components of the bedrock influence the final grain size distribution and mineralogy (Fig. S2, Table 3), which is also reflected in the carbon stock and bulk and pool-specific turnover times. Within the limited geographic and temporal scope of this paper, we hypothesize that for soil carbon stocks and their turnover times, temperature is not the dominant driver, which has been concluded by some (Giardina and Ryan, 2000) but refuted by others (Davidson et al., 2000; Feng et al., 2008). The only climate-related driver which appears to be significant for the deep soil is precipitation.

4.4 Modular robust numerical optimization

The numerical approach used here builds on previous work concerning turnover modeling of bomb-radiocarbon-dominated samples (Herold et al., 2014; Solly et al., 2013; Torn et al., 2009) and the approach used in numerous time-series analyses with box modeling using Excel (Schrumpf and Kaiser, 2015) or Excel Solver (Baisden et al., 2013; Prior et al., 2007). However, certain modifications were made in order to (i) provide objective repeatable estimates, (ii) incorporate longer time-series data, and (iii) identify samples impacted by petrogenic (also called geogenic) carbon. Identifying petrogenic carbon in the deep soil is important considering the large carbon stocks in deep soils (Rumpel and Kögel-Knabner, 2011) and the wider relevance of petrogenically derived carbon in the global carbon cycle (Galy et al., 2008). This approach is modular and could be adapted in the future to identify the correct turnover time for time-series ¹⁴C data, which is becoming increasingly important with the falling bomb-peak signature (Graven, 2015). For the single and time-series data, the results from the numerical solution were benchmarked to the Excel-based model and it was found that the results agree.

Other studies (e.g., Baisden and Canessa, 2012; Prior et al., 2007) also use time-series data to estimate the value for two unknowns simultaneously (size of the pool and turnover time). The error does not always converge to a single low point, but can have multiple minima (Fig. S1). This potential issue should be considered when interpreting the data. More time-series data are required to eliminate this problem.

5 Conclusion

Time-series radiocarbon (¹⁴C) analyses of soil carbon across a climatic range reveal recent bomb-derived radiocarbon in both upper and deeper bulk soil, implying the presence of a rapidly turning-over pool at depth. Pool-specific time-series measurements of the WEOC indicate this is a more dynamic pool which is consistently more enriched in radiocarbon than the bulk. Furthermore, the estimated modeled size of the dynamic fraction is non-negligible even in the deep soil (∼ 0.1–0.2). This could imply that a component of the deep soil carbon could be more dynamic than previously thought.

The interaction between precipitation and geology appears to be the main control on carbon dynamics rather than site temperature. Carbon turnover time in nonhydromorphic soils is relatively similar (decades to centuries) despite dissimilar climatological conditions. Hydromorphic soils have turnover times which are up to an order of magnitude slower. These trends are mirrored in the dynamic WEOC pool, suggesting that in nonwaterlogged (aerobic) soils the transport of relatively modern (bomb-derived) carbon into the deep soil and/or the microbial processing is enhanced as compared to waterlogged (anaerobic) soils.

Model results indicate certain soils contain significant quantities of preglacial or petrogenic (bedrock-derived) carbon in the deeper part of their profiles. This implies that soils not only sequester modern carbon but can also mobilize and potentially metabolize fossil or geogenic carbon.

Overall, these time-series ¹⁴C bulk data and pool-specific data provide novel constraints on soil carbon dynamics in surface and deeper soils for a range of ecosystems.