Controls on Soil Organic Carbon Partitioning and Stabilization in the California Sierra Nevada

Abstract

There is a critical need to quantify the role of soil mineral composition on organic carbon (C) stabilization in forest soils. Here, we address this need by studying a matrix of forest ecosystems and soil parent materials with the objective of quantifying controls on the physical partitioning and residence time of soil organic carbon. We sampled soil profiles across a climate gradient on the western slope of the California Sierra Nevada, focusing on three distinct forest ecosystems dominated by ponderosa pine, white fir, or red fir, on three igneous parent materials that included granite, andesite, and basalt. Results indicated that short-range order mineral phases were the dominant factors accounting for the variation in soil carbon content and residence time. The results further suggested an interaction between ecosystem fire regime and the degree of soil weathering on the partitioning, chemical composition, and residence time of C in density separated soil physical fractions. These results suggest a link between the degree of soil weathering and C storage capacity, with a greater divergence in storage capacity and residence time in the Inceptisols, Entisols, and Andisols of the white fir and red fir ecosystems relative to minimal variation in the highly weathered Ultisols and Alfisols of the ponderosa pine ecosystem.

1. Introduction

Understanding the fate and stabilization of organic carbon (C) in temperate forest soils is central to determining the role of terrestrial ecosystems in mitigating climate change, predicting effects of changing forest management and fire regimes, and quantifying local and regional terrestrial C budgets. Temperate forests contain roughly 118 Pg of C, or approximately 15% of the total C in forested ecosystems globally [1]. Over half of this C is stored belowground, serving as a long-term store (10²–10⁴ yrs) and an active sink of C [2,3]. Here, we quantify how variation in soil parent material and soil mineral assemblage control the physical partitioning and stabilization of soil C across a range of regionally important, highly productive conifer ecosystems in the Sierra Nevada of California. These ecosystems play a significant and disproportionately large role in regional soil C budgets relative to their land area [4], highlighting the critical need to understand the mechanisms controlling soil C stabilization in these systems.

The current state of understanding of soil C dynamics indicates that the storage and stabilization of soil C is controlled by complex interactions among climate, vegetation type, fire regime, topography, and lithology [5,6,7,8,9]. Climate, in the form of precipitation, temperature, and potential evapotranspiration, is considered a central control on primary production and belowground C accumulation, with the assumption that systems evolve toward a steady-state belowground C pool maintained by plant inputs and decomposition [10]. Once in the belowground mineral matrix, C is stabilized through a combination of mechanisms including (i) organic matter chemical composition, (ii) adsorption to mineral surfaces and complexation with metals, and (iii) occlusion of C within aggregates. Vegetation type controls the “quality” or chemical composition of inputs [11]; lithology controls aggregation, organo-mineral association, and microbial community structure and function [11,12,13]; and fire and topography the control production and storage of pyrolyzed C and erosive redistribution of soil C [14,15,16]. At pedon and catchment scales, all of these factors exhibit complex interactive control on soil C stabilization.

Biogeochemical models commonly rely on soil physical parameters, namely clay content, to partition C into respiration and soil C pools of varying size and mean residence time (MRT), with the assumption that the greater reactive surface area associated with clay minerals facilitates C protection via organo-mineral interactions and aggregation [17,18,19]. However, in terms of soil physical properties, substantial literature indicates that the composition of the soil mineral assemblage is more important than clay content in controlling soil C dynamics, particularly the presence of short-range order mineral phases and Al–humus complexes [17,18,19,20,21,22,23]. The mineral-specific preservation is directly related to the chemical and physical properties of the mineral matrix with organo-mineral interaction mechanisms including: (i) the direct sorption of organics to mineral surfaces through entropy-driven reactions, (ii) hydrophobic interactions, van der Waal’s forces, and H-bonding, (iii) and through multivalent cations (such as Ca²⁺, Al³⁺, and Fe³⁺) that bridge the negative charge of mineral and organic surfaces [24]. Variation in the mean residence time of mineral-associated C has been suggested to largely be a function of the mineral assemblage, with a strong correlation between soil C content, residence time, and pedogenic iron oxyhydroxides, short-range order phases such as allophane, imogolite, and ferrihydrite, and Al–humus complexes [25,26,27,28]. These relationships are particularly evident in Andisols, where the soil mineral assemblage is dominated by short-range order phases, and an abundance of Al–humus complexes.

The mineral assemblage also directly controls the abundance and strength of soil aggregates. The aggregation of mineral particles through mineral–mineral and organo-mineral interactions stabilizes belowground C by occlusion within aggregate structures [29,30] where access by the microbial community is limited by physical separation [31]. Additionally, strong oxidative gradients are possible within soil aggregates with limited oxygen availability in interior aggregate regions [32]. Diverse families of C compounds enhance aggregate formation, such as polysaccharides and compounds extractable by polar and non-polar organic solvents [7]. Emphasis has also been placed on the role of hydrophobic compounds in promoting aggregation by repelling water and reducing wettability [33,34,35]. A combination of physical protection and the preferential accumulation of hydrophobic compounds within aggregates may explain the relative stability of aggregate-occluded C.

The combination of density fractionation and physical dispersion techniques isolate C fractions that capture a number of these stabilization mechanisms, and can be used for modeling belowground C dynamics [36,37]. The density fractionation scheme isolates a free light fraction, a light fraction putatively occluded with aggregate structures or coated with minerals, and a mineral-associated C fraction. Chemical properties and the residence time of the C varies significantly across these fractions for forested and agroecosystems [38]. In dry conifer forests of the western United States (US), free light fractions typically consist of intact and partially degraded plant parts and pyrolyzed C with relatively fast turnover, occluded fractions include a mix of organic materials, degraded plant parts, and pyrolyzed C with a relatively long turnover time, and mineral-associated fractions comprise highly altered organic materials and microbial metabolites of intermediate turnover time [12,39]. The relative stability of occluded C in dry western US conifer forests contrasts to results observed in wet temperate and tropical ecosystems [40,41,42], where occluded C generally exhibits residence times similar to that of free light C, and mineral C fractions exhibit the longest residence times.

Our previous work in forested ecosystems in the Sierra Nevada [43,44,45,46] documented significant correlations between short-range order phases and C mineralization, and suggested that mineral assemblage may also control soil C stocks and soil C residence times [22]. The objective of this research was to build from this framework and quantify the physical partitioning and residence time of soil organic carbon and their controlling factors across a range of conifer ecosystems and parent materials, using the western slope of the California Sierra Nevada as a model set of ecosystems. We placed particular focus on the role parent material and soil mineral assemblage exert on controlling C partitioning and stabilization with the hypothesis that soils enriched in short-range order phases exhibit both greater C storage and longer residence time.

2. Materials and Methods

2.1. Study Sites and Field Sampling

The western slope of the Sierra Nevada provides an ideal setting to address parent material and soil mineral assemblage controls on C partitioning and residence time. The western slope encompasses a large environmental gradient that varies predictably with elevation (Figure 1), where temperature decreases with elevation at the same time that precipitation increases and changes from predominantly winter rain to snow at high elevations (

Table 1. Dominant vegetation, climate parameters, and soil taxonomy for sampled ecosystems ¹.

Ecosystem	Dominant Vegetation	Elevation (m a.s.l.)	MAP (mm yr−1)	MAT (°C)	Parent Material	Soil Taxonomy
PP	Pinus ponderosa	920–1400	80–130	10–13	GR	fine-loamy, mixed, semiactive, mesic Ultic Haploxeralf
	Pinus lambertiana		Mostly rain		BS	fine, kaolinitic, mesic Xeric Haplohumult
	Quercus kelloggii				AN	fine, parasesquic, mesic Andic Palehumult
WF	Abies concolor	1500–1800	80–130	8–10	GR	coarse-loamy, mixed, superactive, mesic Humic Dystroxerept
	Pinus ponderosa		Rain/snow		BS	loamy-skeletal, mixed, superactive, mesic Typic Haploxerept
	Pinus lambertiana				AN	medial-skeletal, amorphic, mesic Humic Haploxerand
RF	Abies magnifica	2200–2400	100–130	5–6	GR	mixed, superactive, frigid Dystric Xeropsamment
	Pinus jeffreyi		snow		BS	sandy-skeletal, mixed, superactive, frigid Typic Xerorthent
					AN	medial-skeletal, amorphic, frigid Humic Vitrixerand

¹ MAP—mean annual precipitation; MAT—Ecosystem abbreviations: PP—ponderosa pine; WF—white fir; RF—red fir. Dominant vegetation is listed in order of over-story dominance. Parent material abbreviations: GR—granite; BS—basalt; AN—andesite.

). Vegetation communities change concurrent with climate grading from oak woodlands to alpine and subalpine communities. The absolute elevation boundaries of ecosystem transitions change with latitude, i.e., ecosystem elevation limits decrease with increasing latitude, following macroclimate variation. Historical fire regimes also vary with elevation, grading from a regime typified by frequent, low to mixed intensity fires at lower elevations to one typified by low frequency, high intensity fire at higher elevations [47,48].

We sampled three parent material transects (granite (GR), andesite (AN), basalt (BS)) across elevation transects representing the three primary conifer ecosystems in the region with ponderosa pine (PP; Pinus ponderosa Laws.) at lower elevations with precipitation as rainfall, white fir (WF; Abies concolor (Gord. and Glend.) Lindl.) at midslope elevations with precipitation mixed as rain and snow, and red fir (RF; Abies magnifica A. Murr.) at the highest elevations with precipitation falling as snow. Ecosystems were named for the dominant overstory species and referred to by a combination of parent material and biome name, e.g., andesite and ponderosa pine, ANpp, or basalt and red fir, BSrf, etc. The sampled range of parent materials also translates into strong physicochemical gradients of soil mineral assemblage within a given ecosystem, particularly in terms of pH, selective dissolution extracts, particle size distributions, and soil taxonomy (Table 1 and

Table 2. Summary of depth weighted average (±1 standard error) soil properties for sampled locations ¹.

	Depth (m)	pH 1:1 H2O	Clay (g kg−1)	Amorphous 2 (g kg−1)	Alp (g kg−1)	Alo (g kg−1)	Feo (g kg−1)	Fed (g kg−1)	Feo/Fd
Ecosystem
PP	0.91 ± 0.04 A	5.8 ± 0.2 A	335 ± 61 A	225 ± 94 A	1.9 ± 0.7 A	5.6 ± 2.6 B	2.5 ± 1.1 A	38 ± 13 A	0.1 ± 0.0 B
WF	0.96 ± 0.07 A	6.1 ± 0.1 A	71 ± 6 B	337 ± 130 A	3.6 ± 1.4 A	21.5 ± 9.4 A	3.9 ± 1.3 A	11 ± 3 A	0.4 ± 0.1 A
RF	0.94 ± 0.09 A	5.7 ± 0.3 A	56 ± 9 B	215 ± 84 A	3.4 ± 0.9 A	16.3 ± 7.8 A	3.6 ± 1.2 A	7 ± 2 A	0.5 ± 0.1 A
Parent Material
AN	0.87 ± 0.05 A	6.0 ± 0.1 A	169 ± 103 A	371 ± 64 A	4.7 ± 1.1 A	25.6 ± 7.9 A	5.6 ± 0.6 A	25 ± 12 A	0.4 ± 0.1 A
BS	0.90 ± 0.04 A	6.1 ± 0.1 A	184 ± 115 A	339 ± 73 A	2.7 ± 0.1 AB	14.7 ± 5.5 A	2.7 ± 0.4 AB	26 ± 14 A	0.2 ± 0.1 A
GR	1.04 ± 0.06 A	5.5 ± 0.2 A	108 ± 54 A	67 ± 7 B	1.5 ± 0.5 B	1.8 ± 1.1 B	1.7 ± 0.5 B	6 ± 3 A	0.4 ± 0.2 A

¹ Values followed by a different letter are significantly different (α = 0.05) using a two one-way ANOVAs with ecosystem or parent material as the main effects followed by Tukey’s honestly significant difference post hoc test. Values were formed to meet assumptions of normality and homogeneity of variance; non-transformed means and standard errors are reported in the table. ² Amorphous content determined by the quantitative X-ray diffraction of bulk soil samples.

).

Three pedons at each site were excavated to a depth of 1 m with a spade, and samples were collected by genetic horizon. Sampling sites were constrained to similar landform (summit locations), slope (<10%), aspect (W and SW facing slopes), and canopy position (outside of tree canopy) to minimize landscape and microclimate variability among sites. Soil morphological properties were described for each profile following standard protocols [49], and bulk samples were collected from each genetic horizon. Additionally, samples for bulk density were collected from each horizon using a modified compliant cavity method [50].

2.2. Soil Characterization

All of the soil samples were air dried and sieved to isolate the <2 mm fine earth fraction, and all of the soil characterization conducted on this fraction followed standard laboratory protocols unless otherwise stated [50]. Soil pH was measured on all of the samples at weight to volume ratios of 1:1 (soil:water) solution. Total organic C and N was measured for each fraction using high temperature dry combustion (Carlo-Erba Elemental Analyzer, CE Elantech, Inc., Lakewood, NJ, USA).

Particle size was determined by laser diffraction using a Beckman Coulter LS 13 320 Laser Diffraction Particle Size Analyzer at the University of Arizona, Center for Environmental Physics and Mineralogy. Following pretreatment to remove organics using NaOCl adjusted to pH 9.5 and carbonates using Na-acetate adjusted to pH 5.0, roughly 0.2 g and 0.1 g of sample were weighed into tubes and mixed for 24 h with 5 mL of deionized water using a Thermo Scientific Labquake^® shaker/rotator, followed by the addition of 5 mL of 5% sodium hexametaphosphate solution for an additional 24 h to ensure the dispersion of soil particles prior to particle size analysis. Laser particle size data were reported equivalent to US Department of Agriculture particle size classification.

The mineral composition of bulk soils were identified by quantitative X-ray diffraction (QXRD). Soil samples were pretreated to remove organic matter [51] with a 100-mL solution of 6% NaOCl adjusted pH to 9.5 with HCl, rinsed with deionized water, centrifuged, dried, and mixed. A known amount of internal standard (corundum) was added to each sample to allow the quantitative interpretation of XRD peaks. Samples were ground using a McCrone Micronizing Mill and prepped for XRD following [52]. All of the sample preparation steps were intended to maximize random orientation and increase the exposed surface area of the included minerals. Samples were run as random powder mounts, measured from 5° 2-theta to 65° 2-theta, with a step size of 0.02° 2-theta, a 3-s dwell time per step, spinning at one revolution per second, a 1° divergent anti-scatter slit, a 10-mm divergent mask, a 1° incident anti-scatter slit, a 0.6-mm fixed receiving slit using a PANalytical X’Pert PRO-MPD X-ray diffraction system (PANalytical, Almelo, The Netherlands) generating Cu–Kα radiation at an accelerating potential of 45 kV and a current of 40 mA. The resulting diffractograms were analyzed using Rietveld analyses to identify crystalline mineral phases and abundance [53].

The additional quantitative mineralogical analysis of bulk soils by selective dissolution procedures from previous studies at these locations [44,54,55,56] were combined with the data collected here to complete the soil characterization database. These data included Fe, Al, and Si extracted with Na-pyrophosphate, acid ammonium oxalate, and citrate dithionate following standard protocol [50]. Extractions were non-sequential. Oxalate extracts Al, Fe, and Si (Al_o, Fe_o, Si_o) from organic complexes and short-range order (SRO) Fe-oxyhydroxides (ferrihydrite) and aluminosilicates (allophane and imogolite). Pyrophosphate extracts Al (Al_p) bound in organo-metal complexes, sodium dithionite extracts Fe (Fe_d) from organic complexes and secondary forms of Fe-oxyhydroxides, both crystalline and non-crystalline [57,58] (Table S1).

2.3. Density Fractionation

The physical distribution of soil C was determined for a select set of surface and subsurface horizons for parent material by ecosystem combination following the methods detailed in Rasmussen, Torn and Southard [22]. Depths of horizons selected for density separation were 0–20 cm and 30–60 cm to facilitate the direct comparison of radiocarbon values among the different sites without the confounding effects of radiocarbon variation with depth.

Briefly, 30 g of soil material was mixed with 150 mL of sodium polytungstate (SPT) at a density of 1.8 g cm⁻³, isolating the “free” light fraction. The heavy fraction remaining was re-suspended in SPT and treated with ultrasonic energy of 1500 J g⁻¹ soil (equivalent to 300 J mL⁻¹) using a Branson Sonifier 450 (Branson Ultrasonics, Danbury, CT, USA) with the ultrasonic probe tip inserted 5 cm below the liquid surface during disruption. The rate of energy output was calibrated by measuring the change in temperature of 100-mL of DI H2O after treatment with ultrasonic energy [59]. The sonication treatment releases a combination of light organic matter that was either “occluded” within aggregate structures and/or had substantial mineral coating such that is was dense enough to not float off with the initial density separation. The remaining dense or heavy “mineral” fraction represents soil C strongly associated with mineral surfaces. Each fraction was thoroughly washed to remove excess SPT. Total organic C and N were measured for each fraction using high temperature dry combustion (Carlo-Erba Elemental Analyzer, CE Elantech, Inc., Lakewood, NJ, USA). The naming for the convention used herein for the three fractions resulting from this separation includes “free light fraction” (fLF), “occluded light fraction” (oLF), and “mineral fraction” (mF).

2.4. Radiocarbon Analyses

Following density separation, ¹⁴C abundances were measured on composited bulk soils and the three density fractions. Samples were composited by horizon for each ecosystem by parent material combination. Radiocarbon analysis was facilitated by the Radiocarbon Collaborative, which was a program supported by Michigan Technological University, the USDA Forest Service, and Lawrence Livermore National Laboratory (LLNL). Samples were graphitized at the Carbon, Water, and Soils Laboratory in Houghton, Minnesota (MI). Briefly, each ground sample (~2 mg for fLF and oLF, ~45 mg for mF) was combusted in a sealed glass test tube containing CuO and Ag to generate CO₂ (g) that was subsequently reduced to graphite by heating with H₂ (g) and a Fe catalyst as described by Vogel, et al. [60]. Radiocarbon abundances were measured at the Center for Accelerator Mass Spectrometry at LLNL in Livermore, California (CA) [61] (Table S2). Resulting measurements were corrected for radioactive decay through normalization to the absolute activity of oxalic acid i, which is an international radiocarbon standard. Quality control was accomplished by analyzing a set of additional well-established radiocarbon standards (Oxalic Acid II, Australian National University sucrose, and TIRI wood) with the sample unknowns. Backgrounds were subtracted based on the activity of a coal standard that was processed along with the samples. Reported fraction modern (Fm) values were converted to Δ¹⁴C notation: $Fm=\left[\left({\Delta }^{14}\mathrm{C}/1000\right)+1\right]/\left[exp\left(1950-yr\right)/8267\right]$, where yr is the year of radiocarbon measurement. The Δ¹⁴C values range from positive to negative, with positive values indicating “modern” C, and increasingly negative values corresponding to longer residence times (all of the radiocarbon data reported in Table S2).

2.5. Data and Statistical Analyses

General soil characterization data were compiled and compared by parent material and the ecosystem using profile depth weighted average values for select soil properties. Data selected for summary were chosen to describe variation in physicochemical and mineralogical properties across sites and included profile depth, pH in H₂O, percent clay, amorphous content as determined from QXRD analysis, pyrophosphate extractable Al, oxalate extractable Al and Fe, dithionite extractable Fe, and the ratio of Fe_o/Fe_d. Characterization data were summarized for one profile per site because selective dissolution data were compiled from previous studies, giving n = 9 profiles with 42 horizons. The data were compared using two separate one-way ANOVAs with parent material or ecosystem as the main effects, significance determined at α of 0.05, and means comparison using Tukey’s honestly significant difference (HSD) post hoc test (Table 2). Data were transformed as needed to achieve normal distributions prior to ANOVA; values reported in tables are means of non-transformed data for ease of interpretation.

Soil organic C data were collected for all of the sampled profiles and horizons, with 27 profiles (three per parent material by ecosystem combination) and a total of 124 horizons. Organic C weight percent values for each horizon were converted to C concentration with units of kg m⁻³ as:

$$C_{conc}=\mathrm{C}·{\rho }_b·\left(1-rf\right)$$

(1)

where C_conc is the C concentration [kg m⁻³], C is the weight fraction organic C [kg kg⁻¹], ρ_b is the bulk density [kg m⁻³], and rf is the volume rock fraction described in the field; C_conc values allow for direct comparison of organic C content across horizons with varying bulk density and rock fragment content. Pedon scale C stocks, with units of kg m⁻², were calculated for each sampled pedon as:

$$C_{stock}={\displaystyle \sum }_{i=1}^k{C}_{conc,i}·h_i$$

(2)

where C_stock is pedon C stock [kg m⁻²], k is the total number of horizons i per pedon, and h is the horizon depth [kg].

Carbon data were compared using a two-way ANOVA with parent material and ecosystem as the main effects, and a parent material × ecosystem interaction term. ANOVAs were calculated separately using C_stock, depth weighted average C [g kg⁻¹], and C:N as the dependent terms, with significance determined at α of 0.05 and means compared using Tukey HSD post hoc tests. Data were transformed as necessary to meet assumptions of normality and homogeneity of variance prior to ANOVA with the non-transformed means reported in

Table 3. Summary of depth weighted average (±1 standard error) organic C and radiocarbon content for all of the sampled locations ¹.

	C (g kg−1)	C:N	C Stock (kg m−2)	Δ14C 2 (‰)
Ecosystem
PP	19.1 ± 1.3 A	22.6 ± 0.7 B	16.2 ± 1.6 A	−48.5 ± 28.8 A
WF	24.9 ± 5.4 A	22.8 ± 0.8 B	15.4 ± 2.7 A	−66.0 ± 21.6 A
RF	21.0 ± 4.0 A	26.4 ± 1.5 A	13.2 ± 2.1 A	−56.6 ± 5.0 A
Parent Material
AN	33.3 ± 4.0 A	23.8 ± 0.8 A	22.2 ± 1.4 A	−86.2 ± 34.1 A
BS	20.0 ± 1.0 B	23.3 ± 0.9 A	11.6 ± 0.9 B	−44.3 ± 16.4 A
GR	11.6 ± 1.4 C	24.7 ± 1.7 A	11.0 ± 1.4 B	−40.6 ± 28.8 A
Parent Material × Ecosystem
ANpp	18.9 ± 1.8 C	22.2 ± 1.1 BC	19.9 ± 2.8 AB	−104.0
BSpp	22.4 ± 0.7 BC	23.7 ± 0.7 BC	14.0 ± 0.9 ABC	−34.0
GRpp	16.0 ± 2.5 CD	22.0 ± 1.8 BC	14.9 ± 3.5 ABC	−7.5
ANwf	45.4 ± 3.2 A	22.3 ± 0.5 BC	25.6 ± 2.1 A	−107.7
BSwf	19.0 ± 1.2 C	24.9 ± 1.9 BC	10.4 ± 2.0 C	−35.7
GRwf	10.3 ± 1.7 DE	21.1 ± 0.8 C	10.1 ± 0.8 C	−54.7
ANrf	35.6 ± 1.1 AB	26.9 ± 0.1 AB	21.0 ± 1.8 A	−46.9
BSrf	18.7 ± 2.3 C	21.2 ± 1.4 C	10.5 ± 0.5 BC	−63.2
GRrf	8.6 ± 0.3 E	31.0 ± 0.5 A	8.0 ± 0.5 C	−59.7

¹ Values followed by a different letter are significantly different (α = 0.05) using two one-way ANOVAs with ecosystem or parent material as the main effects, followed by Tukey HSD post hoc test. Values were formed to meet the assumptions of normality and homogeneity of variance; non-transformed means and the standard errors are reported in the table. ² Radiocarbon content was only measured for one profile in each parent material × ecosystem, so the values presented here are the depth weighted average values for that one profile.

.

The density fractionation was performed on one surface and one subsurface horizon from each of the three replicate profiles sampled at each parent material by ecosystem combination. The fraction data summarized in Table S3 include the C content in each fraction, fraction C:N, and the relative proportion of total C in each fraction. Separate two-way ANOVAs were performed on surface and subsurface samples with parent material and ecosystem as the main effects with a parent material × ecosystem interaction term. Radiocarbon was measured on one sample composited from the three replicates for each fraction, and hence were not included in the full two ANOVA.

We also performed regression analysis to determine the factors that explained the most variance in soil C content [g kg⁻¹] and Δ¹⁴C [‰]. Soil C and Δ¹⁴C were regressed against a suite of independent soil physicochemical and mineralogical variables, including horizon midpoint depth, percent clay, Fe_d, Fe_o, Al_o, Al_p, Fe_o/Fe_d, percent amorphous content as determined by QXRD, and pH in H₂O, using mixed-model regression including both fixed and random effects. Regression analyses were performed using average soil parameter values for each parent material by ecosystem factor with a total of 42 separate horizons. Horizon depth was included as a fixed effect in all of the models, and the parent material by ecosystem factor was set as the random effect to account for the lack of independence among multiple horizons from one profile. All of the data were transformed to normal distributions as needed, and standardized to a mean of zero and standard deviation of one prior to regression analysis.

We initially performed simple exploratory linear regression analyses between horizon depth and soil C and Δ¹⁴C, and then compared the residuals of these regressions to the various soil physicochemical/mineralogical properties using Spearman’s correlation to guide the selection of fixed effects to include in the mixed models (Table S4). The correlation analyses indicated all of the oxalate extractable parameters, Al_p, and the percent of amorphous material by QXRD were positively correlated with C content, and also highly collinear, with r > 0.7. Fe_o exhibited the strongest correlation to the soil C by depth residuals, and was thus chosen as the best parameter for describing the amorphous content relationship to soil C. Based on these analyses, we performed a mixed regression model using depth and Fe_o as the fixed effects, and parent material by ecosystem as the random effect. A similar analysis for Δ¹⁴C also indicated that Fe_o was the overall best predictor, with a final mixed model where depth and Fe_o were the fixed effects, and parent material by ecosystem was the random effect (

Table 4. Fixed-effect parameters from the linear mixed models.

	Model Parameter	Regression Coefficient	95% Lower	95% Upper	F-Ratio	Prob > F
Bulk C (g kg−1)	Depth	−0.74	−0.89	−0.6	111.1	<0.0001
	Feo	0.34	0.11	0.58	11.4	0.0103
	Intercept	0.003	−0.23	0.24	-	-
Bulk Δ14C (‰)	Depth	−0.91	−1.05	−0.76	161.1	<0.0001
	Feo	−0.21	−0.37	−0.06	10.4	0.0109
	Intercept	0.001	−0.15	0.15	-	-

).

All of the statistical analyses were performed in JMP Pro v13.0 (SAS Institute Inc., Cary, NC, USA), graphs were generated using SigmaPlot v11.0 (Systat Software Inc., San Jose, CA, USA), and maps were produced in ArcGIS v10.3.1 (ESRI, Redlands, CA, USA).

3. Results

3.1. General Soil Properties

General soil properties and taxonomic variation have been detailed in previous studies, with a brief summary here by ecosystem type and each parent material to provide context (Table 1 and Table 2). Generalizing by parent material and ecosystem indicated relatively few statistical differences because of the large variability in properties when grouping across the large biogeochemical gradients that these sites encompass. The degree of weathering and mineral alteration was greatest in the relatively warm, wet PP systems, as evidenced by the significantly greater clay content, a greater proportion of pedogenic Fe oxyhydroxides dominated by crystalline forms, and a relative dominance of the secondary phyllosilicates by kaolin and gibbsite. The PP soils were all classified as Ultisols or Alfisols, reflecting their advanced weathering state. The WF systems presented intermediate levels of weathering and mineral alteration, with soils classified as Andisols and Inceptisols. The RF soils expressed the least amount of weathering with the lowest clay percentages and pedogenic Fe-oxyhydroxides, and were classified as Entisols or Andisols. Generalizing by parent material indicated significantly greater short-range order components in AN soils, generally greater clay content in BS derived soils, and more acidic pH values in GR soils.

3.2. Bulk Carbon Concentration and Stocks

Soil C content and stocks varied significantly across parent materials, ecosystems, and the parent material by ecosystem interaction term (Table 3). There was little variation among ecosystems, but significant variation among parent materials. AN soils had an average soil C content of 33.3 g kg⁻¹, which was 1.5 times that of the BS soils with an average of 20.0 g kg⁻¹, and three times greater than that of the GR soils with an average of 11.6 g kg⁻¹. Soil C stocks were also significantly greater in AN soils with an average of 22.2 kg m⁻², and twice that of stocks in the BS and GR soils. Parsing by the parent material by ecosystem interaction term indicated that ANwf soils contained the greatest soil C content, with an average of 45.4 g kg⁻¹, and C stocks, with an average of 25.6 kg m⁻². The least amount of soil C was found in the GRrf soils, with an average soil C content of 8.6 g kg⁻¹ and C stock of only 8.0 kg m⁻². Of note is that C:N tended to be greater in RF ecosystems, with the highest average value of 31.1 in the GRrf soils.

Carbon concentrations (kg m⁻³) expressed clear variation with depth and among ecosystems (Figure 2). All of the profiles indicated that C concentration decreased with depth. In the WF and RF ecosystems, the AN soils showed greater C concentrations with depth than both the GR and BS soils, which followed similar depth patterns to one another.

Regression analysis indicated the dominant physicochemical control on soil C content (g kg⁻¹) was Fe_o (Table 4). All of the oxalate extractable components, Al_p, and amorphous content by QXRD were positively correlated with soil C content (Table S4), but only Fe_o was included in the regression model because of issues with parameter collinearity.

3.3. Bulk Carbon Radiocarbon Content

Bulk Δ¹⁴C decreased with depth in all of the profiles, with the surface values all presenting a modern signature with positive values ranging from 0‰ to +125‰, and declining to values between −75‰ to −250‰ in the deepest horizons (Figure 3). The AN soil Δ¹⁴C values were generally more negative than the BS and GR soils at a given depth across all of the ecosystems, but particularly in the WF and PP ecosystems. This pattern was also present in the depth weighted average parent material by ecosystem Δ¹⁴C data, where AN soils exhibit substantial radiocarbon depletion relative to the BS and GR soils (Table 3). The RF soils exhibited minimal variation among averages or by depth. Averaged by the main effects of ecosystem or parent material indicated no significant variation due to the large variability in Δ¹⁴C values. Regression analyses indicated that along with depth, the best predictor of bulk Δ¹⁴C was Fe_o (Table 3).

3.4. Physical Carbon Partitioning and Radiocarbon Content

The fraction data, when pooled by fraction across all of the samples, indicated significant variation in fraction C content, C:N, Δ¹⁴C, and the relative partitioning of C to each fraction (Figure 4). The oLF exhibited greater C enrichment with a median C content of 400 g kg⁻¹, the highest C:N values with a median of 44, and the oldest Δ¹⁴C signature with a median of −59 ‰; however, the least amount of total C was partitioned to the oLF. In contrast, the mF exhibited the lowest C content and C:N, with Δ¹⁴C comparable to the oLF, while accounting for a median of 34% of the total soil C.

Similar patterns were observed when data were split by parent material (Figure 5) and ecosystem (Figure 6) with some key observations. AN soil C was predominantly partitioned to the mF fractions, and also exhibited the highest C content in the mF. When parsed by ecosystem, there was a trend of greater fLF with increasing elevation, grading from PP to WF to RF ecosystems. Of note is that the PP ecosystem exhibited the greatest relative partitioning of C into the oLF and also showed relatively high C content and C:N in both the fLF and oLF.

There was generally little variation in fraction C content by parent material and ecosystem in the fLF and oLF, with fLF ranging from 20% to 32% C and oLF ranging from 29% to 45% C (Table S3). The mF fraction did exhibit significant variation, with AN soils containing the greatest mF C content across all of the ecosystems, and surface and subsurface horizons. The greatest difference was observed in the ANwf surface and subsurface horizons, where AN soils had on the order of three times the mF C relative to the GR and BS mineral fractions.

Fraction C:N generally did not vary within a fraction across parent materials and ecosystems in the surface horizons, but did exhibit some significant differences among parent materials and ecosystems in the subsurface. In particular, the C:N of the fLf and oLf in the PP ecosystems were relatively high, with oLF values ranging from ~59 to 69. Similarly, PP mF C:N values were significantly lower across all of the parent materials relative to the WF and RF ecosystems.

4. Discussion

Carbon concentration, total C stocks, and residence times were generally greatest in the andesite-derived soils across all of the ecosystems. This pattern was most pronounced in the white fir ecosystem where the andesite-derived soils classified as Humic Haploxerands with an amorphic mineralogy class. These findings agree well with previous work comparing Andisols and andesite-derived soils, indicating that these soils preferentially store and stabilize organic carbon relative to non-andesitic soils [21,22,62,63]. Relative differences in C stocks and residence time among parent materials were greatest in the white fir and red fir ecosystems, particularly in the red fir ecosystem where andesite soils contained two to three times the total C stock relative to the basalt and granite-derived soils, respectively. The red fir ecosystem is the coldest and wettest ecosystem, with a relative lack of mineral weathering and subsurface structural development. In contrast, soil C stocks exhibited minimal variation among parent materials in the warm, wet ponderosa pine ecosystems with highly developed soils in terms of clay content, degree of mineral alteration, and subsurface structural development. The relative variation in differences among the ecosystems variation suggests a link between degree of soil development and variation in C storage capacity among parent materials, with greater divergence in storage capacity in the Inceptisols, Entisols, and Andisols of the white fir and red fir ecosystems relative to the highly weathered Ultisols and Alfisols of the ponderosa pine ecosystem. These patterns of variation in C storage and residence with degree of soil development correspond with similar patterns observed across soil chronosequences on volcanic and mixed alluvial parent materials demonstrating a peak in C storage and residence time in moderately developed Andisols and Inceptisols [23,25,64].

Regression analyses for both C and radiocarbon content indicated that short-range order materials were the soil physicochemical parameters that accounted for the greatest variance in both, suggesting that these phases represent the dominant control on C stabilization across these ecosystems. In this study, oxalate extractable Fe presented the strongest relationship, and thereby served as a representative proxy for the other short-range order phases, but all of the measures of short-range materials and amorphous phases were correlated with both C and radiocarbon content. In contrast, clay content exhibited a minimal relation to either C content or residence time. These findings are similar to those in other large data syntheses [28,65,66,67] demonstrating that short-range order phases and mineral type are more important controls on soil C content than clay content, particularly in cool, wet, forested Inceptisols and Andisols. Mechanistically, the short-range order phases provide a high reactive surface area that promotes strong organo-mineral interactions and aggregation, leading to C storage and stabilization [68,69,70].

The physical partitioning data indicated the majority of C is located in a combination of the free light and mineral fractions, the relative distribution of which varies significantly by ecosystem, parent material, and depth. The contribution of the free light fraction increases substantially with elevation from ponderosa pine to red fir-dominated ecosystems, with a less pronounced trend of decreased partitioning to the mineral fraction with elevation. The variation in partitioning with elevation is most likely related to the greater abundance of clay-sized materials and potentially reactive surface area in the ponderosa pine Ultisols and Alfisols relative to the clay-poor Entisols and Andisols in the red fir ecosystems. It is important to note that the relative partitioning of C to the mineral fraction was greatest in andesite soils despite not having significantly greater clay size materials, which is a trend attributable to the greater abundance of short-range order materials.

The free light fraction generally exhibits relatively high C content and C:N ratios that are typical of plant materials, and indeed visual observation of the free light fraction indicates the presence of identifiable plant parts, similar to Wagai, Mayer and Kitayama [38], who noted that free light fractions consist of partially decayed organic materials in a synthesis of density fractionation studies. This fraction also generally exhibits the greatest radiocarbon enrichment, with positive Δ¹⁴C values indicative of modern carbon, suggesting a dominance of this fraction by recent inputs with rapid turnover. In contrast, the mineral fractions exhibit low C content and C:N values <20, and radiocarbon values that are more negative than the corresponding free light fractions. The low C:N of the mineral fraction indicates that it consists predominantly of microbially processed organic matter, which is likely some combination of both live and dead microbial biomass and organo-mineral complexes, with a combination of mineral protection and the recycling of older C within the microbial community, contributing to the older radiocarbon values.

Carbon content and C:N are generally greatest in the occluded fraction, and this fraction generally exhibits the most depleted, and therefore oldest, radiocarbon signatures. The high C content and C:N in the occluded fractions suggest a substantial contribution of pyrolyzed carbon. Visual examination of the occluded fractions indicates that the organic matter is very dark to black in color and fine grained, which is suggestive of charred materials. Furthermore, the greatest proportion of occluded C is found in the ponderosa pine ecosystems, particularly in surface horizons. The ponderosa pine ecosystem has a fire regime typified by frequent, low-intensity fires that combust surface litter and understory plant growth, and foster frequent inputs of charred material to surface horizons. Heckman et al. [12] observed a similar trend in ponderosa pine-dominated forests of Arizona where the occluded fraction contained from 12% to 63% pyrolyzed C based on ¹³C-NMR analysis, and was enriched two to five times in pyrolyzed C relative to bulk soil C. The ponderosa pine ecosystem free light fractions also exhibit generally higher C content and C:N ratios that suggest a contribution of pyrolyzed C to this fraction as well. These findings agree well with other studies in dry western conifer ecosystems [22,39] that indicate an important contribution of pyrolyzed C to the occluded fraction.

The occluded fractions in the white fir and red fir ecosystems contribute to a smaller proportion of total C, but still exhibit characteristics of pyrolyzed C in terms of their visual appearance, C content, and C:N ratio. This pattern follows the change in the fire regime with elevation whereby the frequency of fire decreases, but the intensity increases. Therefore, we speculate that across all of these systems, the occluded organic matter is dominated by pyrolyzed materials, and that the longer residence times of C in these fractions is a function of both C chemistry and physical occlusion by mineral phases. The change in fire regime from low intensity to high intensity may also affect the chemistry and chemical stability of the pyrolyzed C preserved in the occluded and free light fractions. The role of fire in local soil C cycling will play a growing role as climate change, land use, and other management factors alter the intensity and return intervals, erosion, and C flux to the atmosphere [71,72,73,74].

5. Summary

The effect of parent material mineral composition and its degree of weathering were important factors controlling soil C content and age. The key findings on the factors influencing soil C and its age in these systems included: (i) significant variation in the partitioning and residence time of soil C was directly related to variation in soil secondary mineral assemblage, particularly short-range order materials, and the degree of soil weathering and subsurface physicochemical alteration; (ii) oxalate extractable Fe was highly correlated to soil C content, age, and partitioning and should be considered a proxy to describe soil C stabilization potential; (iii) clay was not a significant predictor of storage or residence time, highlighting the need to include parameters other than clay in soil biogeochemical models, particularly in forested ecosystems; and (iv) the data suggested a significant role of fire in the generation of occluded materials and mineral assemblage in partitioning soil C to mineral and free light fractions.

The results from this study indicating a limited to no role for clay content and a central role for short-range order phases in the control of soil C content and residence time demonstrate that we need to “rethink” the basic assumption of biogeochemical models to include mineralogical and transition metal data, particularly Fe, in describing the fate of soil C in temperate forest soils. The results also indicated that the degree of variation in C content, partitioning, and stabilization varied among soil orders, with the highly weathered Ultisols exhibiting minimal variation, whereas the weakly weathered Entisols exhibited substantial variation. This was driven by the relative convergence of soil physicochemical properties in the highly weathered soils. Taken in a broader context, these results suggest there may be consistent trends in soil C stabilization mechanism type and variability among soil taxonomic units that may be used to better inform biogeochemical models. Finally, the results from this study, while indicating an important role for fire regime in the relative partitioning of soil C, the interaction among mineral and physicochemical controls with ecosystem fire regime and the amount and type of pyrolyzed C produced, and how this impacts soil C dynamics remains unclear. Disentangling these interactions is needed in order to improve understanding of soil C dynamics in forested ecosystems.