Does elevated atmospheric CO₂affect soil carbon burial and soil weathering in a forest ecosystem?

Sampling and sample preparation

Soil samples were collected from the ORNL FACE site ten years into the experiment in September 2008. We used a sharpened steel pipe (4.8 cm diameter) driven into the soil with a nylon-face mallet (Sánchez-de León et al., 2018) to obtain four cores from each of the ambient (control) plots and four from each of the elevated CO₂ plots. Intact soil cores were stored frozen and sectioned with a thin ice-core saw (while frozen and the blades cleaned between cuts) as follows: the top 8 cm of each core was sectioned into 1-cm depth increments, and from 8 to ∼20 cm depth the core was sectioned into 2-cm depth increments. Additional soil core samples from 20 to 25 and 25 to 30 cm were collected adjacent to each sampling spot to help constrain the maximum depth of measurable ¹³⁷Cs activity. No samples were collected below 30 cm depth for these experiments. We compared our samples with soils samples collected in 1997 (prior to the initiation of FACE experiment) and archived. Pre-treatment core samples were for depth ranges of 0–5, 0–15, 15–30, and 30–45 cm from both the ambient and elevated-CO₂ plots.

After sectioning the soil cores, rocks and roots were manually removed from each section. Samples were dried at 80 °C, gently crushed and sieved to pass through a 2-mm sieve. Dry bulk densities were calculated from separate samples by comparing the 2 mm-sieved soil dry weight with the core section volume after correction for the occasional small rocks being removed.

Soil organic carbon concentration and stable C isotope ratios

Soil samples were ground to a fine powder for analysis of organic C concentration and stable C isotope ratios. Carbonates were removed before analyses as explained elsewhere (O’Brien et al., 2015). Analyses were performed at the Ecology Stable Isotope Laboratory (UIC) using a Costech ECS 4010 elemental analyzer with a zero-blank autosampler interfaced with a ThermoFinnigan Delta-Plus XL isotope ratio mass spectrometer in continuous flow. Soil organic C concentrations are reported in %C (dry weight basis). The ¹³C/¹²C isotope ratios are reported in the conventional delta notation, in units of per mil relative to the standard reference material VPDB (Coplen, 1996), according to: (1)${\displaystyle {\delta }^{13}\mathrm{C},‰=\left[\left(R_{\mathrm{sample}}∕R_{\mathrm{V\; PDB}}\right)-1\right]\times 1,000}$ where R is the atom ratio ¹³C/¹²C. Reproducibility of δ¹³C values is better than ±0.1‰  when compared to international standards.

Gamma spectrometry

Gamma spectrometry was performed at the Environmental Isotope Geochemistry Laboratory (UIC) by using a Canberra model GR3020 reverse-electrode intrinsic Ge detector system interfaced with a DSA-2000 digital spectrum analyzer. Dry homogenized sediment samples (5–10 g) were weighed into aluminum counting cans and these were sealed with Al foil. Gamma activities for ⁴⁰K, ¹³⁷Cs, ²¹⁰Pb, and ²²⁶Ra were measured at 1,460.5, 661.6, 46.5, and 186.2 keV, respectively, with cans centered on top of the detector. Detector efficiency was calibrated versus sample weight in the same geometry using the certified standards CANMET DL-1a (U-Th ore diluted in quartz sand) and NIST SRM-4357 (Ocean Sediment). Relative uncertainties of measured activities were less than ±10% for activities >4 Bq kg⁻¹, as calculated from counting statistics incorporating background subtraction and propagated errors. Activities were measured per sample dry mass and normalized to dry bulk density measurements for reporting in units of Bq cm⁻³ or Bq kg⁻¹.

Bioturbation model based on ¹³⁷Cs profiles

Mathematical models combining advection and diffusion have been developed to explain the downward movement of ¹³⁷Cs and the diffusion-like broadening of its profile in soils and sediments (Guinasso & Schink, 1975; Olsen et al., 1981; Robbins, 1986). The steady-state bioturbation model (Eq. (2) in Robbins, 1986) describes the total concentration [C(x, t)] of particle-bound radionuclides in the soil as a function of the vertical distance (x), and time (t). Biological agents and advection explain the downward transport of ¹³⁷Cs (see Eq. (2)). The biodiffusion coefficient (D_b) describes diffusive mixing of bulk soil by biological agents. Transport of ¹³⁷Cs can also be caused by advective processes involving motion of particles and pore fluid (ν). The net loss or gain of the ¹³⁷Cs within the soil profile is accounted for by the radioactive decay constant (λ) and the first-order feeding rate constant that describes net transport of ¹³⁷Cs by moving organisms (γ): (2)${\displaystyle \frac{\partial C}{\partial t}=\frac{\partial }{\partial x}\left(D_b\frac{\partial C}{\partial x}\right)-\frac{\partial }{\partial x}\left(\nu C\right)-\left(\lambda +\gamma \right)C.}$

Best-fits for the average ambient and elevated-CO₂ ¹³⁷Cs profiles used the following fixed values for the model (Eq. (2)): 0.75 Bq cm⁻² for the initial activity of ¹³⁷Cs (C₀) (Hardy et al., 1968; CDC-NCI, 2005); 45 years for time elapsed (t) between ¹³⁷Cs tracer deposition in 1963 and sample collection in 2008 (or 34 years for 1997); and, 0.023 yr⁻¹ for the ¹³⁷Cs decay constant (λ). For a pulse-like input of the tracer, the model represented in Eq. (2) has a well-known solution, which has been widely used to describe ¹³⁷Cs profiles in soils (Van Genuchten & Cleary, 1979; Ivanov et al., 1997; Bossew & Kirchner, 2004; Schuller et al., 2004). (3)${\displaystyle C\left(x,t\right)=C_0{e}^{-\left(\lambda +\gamma \right)t}\left\{\frac{1}{\sqrt{\pi D_{b}t}}{e}^{-{\left(x-\nu t\right)}^2∕\left(4D_{b}t\right)}-\frac{\nu }{2D_b}{e}^{\nu x∕D_b}\text{erfc}\left(\frac{x+\nu t}{2\sqrt{D_{b}t}}\right)\right\}.}$

We calculated error-weighted least-squares best fits of Eq. (3) to the data for each of our measured ¹³⁷Cs profiles by using MATLAB.

Bioturbation model based on unsupported ²¹⁰Pb profiles

Mathematical models combining advection and diffusion terms to account for downward transport and dispersion of unsupported ²¹⁰Pb differ from those for ¹³⁷Cs, because ¹³⁷Cs is deposited in a pulse-like manner whereas ²¹⁰Pb is deposited continuously as it is produced from decay of atmospheric ²²²Rn (Robbins, 1978). We used the following steady-state equation to describe advective-diffusive transport of ²¹⁰Pb (Kaste et al., 2011), where A_Z is the initial activity of unsupported ²¹⁰Pb at the surface (Bq/kg); A_Z is the activity of unsupported ²¹⁰Pb at depth z (cm); v is the advection rate (4)${\displaystyle A\left(z\right)=A_0\exp \left[\frac{\nu -\sqrt{{\nu }^2+4\lambda D}}{2D}\left(z\right)\right].}$

(cm yr⁻¹); D is the diffusion constant (cm yr⁻¹); and λ is the decay constant of ²¹⁰Pb (0.031 yr⁻¹).

Soil organic carbon and δ¹³C profiles

Soil organic carbon content was highest at the surface and decreased with depth (Fig. 1A). The top 6 cm of the elevated-CO₂ profiles, on average, have δ¹³C values significantly lower than in the ambient profiles (Fig. 1B), as indication of SOC inputs since the initiation of the experiment in 1998. As a result of new inputs, the top 2 cm of the elevated-CO₂ profiles, on average, have significantly higher SOC content than in the ambient profiles (p < 0.05), but below 2–4 cm depth the average profiles are not significantly different (Fig. 1A). When the δ¹³C values are compared with the inverse SOC content, it is apparent that the average elevated-CO₂ profile is depleted in ¹³C. The SOC being deposited at the soil surface has a δ¹³C value of about −38‰  in the elevated-CO₂ plot compared with −28‰  in the ambient plot (Fig. 2).

Soil bulk density was lower at the top 5 cm of the soil profile than at the rest of the soil depths (Table 1). Bulk density increased from values of about 0.5 g cm⁻³ at shallow depths to values greater than 1 g cm⁻³ at 5–6 cm depth and deeper. This may reflect degradation and mineralization of the litter layer which occurs during the first decades following deposition and produces denser residual material (Kaste et al., 2011).

¹³⁷Cs profiles

Detectable ¹³⁷Cs was measured from the surface to a depth of at least 20–30 cm in all soil profiles (Fig. 3). Total ¹³⁷Cs inventories of the soil profiles are less than or equal to that expected if the assumed initial activity of ¹³⁷Cs (0.75 Bq cm⁻³) remained in place and decayed for a period of 45 years from deposition in 1963 to sampling in 2008. Maximum measured activity for ¹³⁷Cs was 27.3 ± 0.6 mBq cm⁻³. Profiles of ¹³⁷Cs activity generally increase with depth from activities of about 2–6 mBq cm⁻³ at the surface to maximum activities at around 8–14 cm depth, followed by a general decrease to values of 0.5–2 mBq cm⁻³ at 30 cm depth.

Pre-treatment core samples for depth ranges of 0–5, 0–15, 15–30, and 30–45 cm collected in 1997 from both the pre-treatment ambient and elevated-CO₂ plots (prior to initiation of the FACE experiment) had cumulative ¹³⁷Cs activities equal to those of the post-treatment samples collected in 2008 (Fig. 4). There was no measurable activity of ¹³⁷Cs beyond 30 cm depth before and during the FACE experiment (Figs. 3 and 4).

For the 2008 samples, more ¹³⁷Cs activity was found at greater depth in elevated CO₂ plots when compared to ambient control plots (Fig. 3; p < 0.05). There was measureable ¹³⁷Cs activity at the 25–30 cm in soil samples collected in the elevated CO₂ plots, whereas the average depth of the deepest measurable ¹³⁷Cs activity in the ambient plots was 20.6 ± 2.5 cm. Similar values of ¹³⁷Cs activity were found for the 15–30 cm soil samples collected in 1997.

Bioturbation model

Bioturbation derived mixing rates were not significantly different between the ambient and elevated CO₂ plots (Table 2; Fig. 3). We solved Eq. (3) for the biogenic diffusivity (D_b), particle advection velocity (ν), and feeding rate constant (γ) values. Advection velocities are indistinguishable for both treatments with an average value of 0.18–0.19 cm yr⁻¹ (Table 2). The feeding rate constants were 0.008 ± 0.003 yr⁻¹ for ambient and 0.005 ± 0.001 yr⁻¹ for elevated-CO₂ plots. The D_b values were at 0.53 ± 0.20 cm² yr⁻¹ at ambient and 0.63 ± 0.29 cm² yr⁻¹ at elevated-CO₂. These biogenic diffusion coefficients (D_b) were used for calculating mixing time constants (τ) for a soil layer thickness L = 20 cm ($\tau =L^2{D}_b^{-1}$) (Kaste, Heimsath & Bostick, 2007). The top 20 cm layer of soil at the ORNL FACE site has estimated average mixing times ranging from about 640 to 750 years (Table 2) but with wide spatial variability.

²¹⁰Pb profiles

The activity ratio (²¹⁰Pb/²²⁶Ra) is a good indicator of excess (or deficient) ²¹⁰Pb relative to that expected from secular equilibrium with ²²⁶Ra (at secular equilibrium, (²¹⁰Pb/²²⁶Ra) = 1). The average (²¹⁰Pb/²²⁶Ra) activity ratio profiles in the ambient and elevated-CO₂ plots are similar, showing excess ²¹⁰Pb in the upper 5-to-10 cm and a deficit of ²¹⁰Pb below 10 cm depth (Fig. 5). Best-fit solutions of Eq. (4) to the (²¹⁰Pb/²²⁶Ra) profiles all yielded lower values of diffusion constant (near 0) and advection rate (∼0.9 cm yr⁻¹) than did the ¹³⁷Cs models. We show the best-fit steady-state advection-decay model in comparison with the mean (²¹⁰Pb/²²⁶Ra) profiles in Fig. 5. The parameters in this model were a constant initial (²¹⁰Pb/²²⁶Ra) value of 2.3 and a steady-state (²¹⁰Pb/²²⁶Ra) ratio of 0.75 at depths below 20-to-24 cm, where excess ²¹⁰Pb has decayed to <2% of its initial amount (Fig. 5). The steady-state value of 0.75, representing a 25% loss of in situ ²²²Rn production, is based on a survey of ²²²Rn loss in 119 soil cores from undisturbed landscapes in North America. As with the ¹³⁷Cs profiles, no significant difference in mean (²¹⁰Pb/²²⁶Ra) profiles is evident between the ambient and elevated CO₂ plots.

Soil mixing dynamics

There is a rapid decrease in (²¹⁰Pb/²²⁶Ra) ratios in the ambient and elevated-CO₂ plots at the FACE site, relative to that predicted by the simple ²¹⁰Pb advection-decay model (Fig. 5). The constant deposition of ²¹⁰Pb from the atmosphere to the soil surface creates a condition of radioactive disequilibrium where ²¹⁰Pb in the shallow parts of soil profiles is in excess of that produced in situ by decay of ²²⁶Ra and intermediate daughters. The profile of excess ²¹⁰Pb can be modeled in terms of soil or sediment accumulation and erosion rates and mixing parameters (Olsen et al., 1981; Robbins, 1986; Kaste, Heimsath & Bostick, 2007; Kaste et al., 2011; Matisoff & Whiting, 2012). The ²¹⁰Pb/²²⁶Ra profiles depicted in Fig. 5 decay too rapidly with depth to be consistent with simple downward advection at 0.18 cm yr⁻¹ and radioactive decay of ²¹⁰Pb. Using Eq. (4), we obtained a best-fit value for advective transport (v) of about 0.9 cm yr⁻¹. The relatively low and constant value of the ²¹⁰Pb/²²⁶Ra activity ratio at depth indicates diffusive escape of ²²²Rn (likely via transpiration stream and soil porosity) that was produced in situ from ²²⁶Ra decay. This ²²²Rn escape is possibly enhanced by bioturbation and transpiration occurring in the shallow root zone (top 10 cm) where the bulk density is the lowest (Table 1). The ¹³⁷Cs profiles and the rapid decrease in excess ²¹⁰Pb suggest a distinct soil boundary at about 4 cm deep, below which most of the ¹³⁷Cs activity resides (Fig. 5). This 0–4 cm depth soil layer is also evident in isotope profiles shown in Figs. 1 and 6. This soil multi-isotope boundary at 4 cm depth is consistent with the enhanced SOC accumulation at elevated CO₂ conditions when compared to ambient seen at the site between 0 and 5 cm (Fig. 1 of Jastrow et al., 2005). Further, these results are also consistent with the presence at the site of endogeic earthworms (Sánchez-de León et al., 2018), which avoid the soil surface likely preventing predation or competition with litter layer fauna.

A second soil isotope boundary is detected at about 16 cm deep, where the ²¹⁰Pb/²²⁶Ra activity ratio approaches the typical disequilibrium value of 0.75 seen in deeper soils (Graustein & Turekian, 1990) (Fig. 5). This 4–16 cm soil section may represent a soil mixing compartment influenced by root proliferation and earthworm activity, potentially redistributing and homogenizing SOC concentration within this depth range (Fig. 1A). This may partly prevent detection of soil C accrual in response to elevated CO₂ levels at these depths (Fig. 1A), despite the isotopic evidence for input of new SOC (Fig. 1B). More research using multiple radioisotope tracers to detect soil profile mixing sections may allow better determinations of C dynamics than are possible by using the traditional arbitrary depth comparisons.