Deep in the Sierra Nevada critical zone: saprock represents a large terrestrial organic carbon stock

Large uncertainty remains in the spatial distribution of deep soil organic carbon (OC) storage and how climate controls belowground OC. This research aims to quantify OC stocks, characterize soil OC age and chemical composition, and evaluate climatic impacts on OC storage from the soil surface through the deep critical zone to bedrock. These objectives were carried out at four sites along a bio-climosequence in the Sierra Nevada, California. On average, 74% of OC was stored below the A horizon, and up to 30% of OC was stored in saprock (friable weakly weathered bedrock). Radiocarbon, spectroscopic, and isotopic analyses revealed the coexistence of very old organic matter (OM) (mean radiocarbon age = 20 300 years) with relatively recent OM (mean radiocarbon age = 4800 years) and highly decomposed organic compounds with relatively less decomposed material in deep soil and saprock. This co-mingling of OM suggests OC is prone to both active cycling and long-term protection from degradation. In addition to having direct effects on OC cycling, climate indirectly controls deep OC storage through its impact on the degree of regolith weathering (e.g. thickening). Although deep OC concentrations are low relative to soil, thick saprock represents a large, previously unrealized OC pool.


Introduction
Soil is the largest terrestrial organic carbon (OC) reservoir on earth, storing over 3000 Pg of C (Le Quéré et al 2014). Most soil OC inventories focus on soil horizons near the surface, and fewer studies report OC below soil in materials such as saprock, even though 60% of OC is found below 1 m (Richter and Markewitz 1995, Jobbágy and Jackson 2000, Jackson et al 2002, Harper and Tibbett 2013. Although key works have expressed the importance of deep soils in understanding the OC cycle, 60% of papers published in four major soil science journals did not specify the depth sampled and the remaining studies had an average depth of only 24 cm (Yost and Hartemink 2020). Investigations of deep OC dynamics (>1 m) are rare in part because the thickness of overlying soil imposes unique sampling challenges related to cost, labor, and accessibility of equipment. Evidence suggests constraining OC stock estimates to the upper layers of soil has led to an underestimation of global soil OC stocks, limiting our understanding of how the entire regolith OC responds to changes in climate.
Compared to subsoil, OC concentration in nearsurface horizons is typically higher with shorter mean residence times (Jobbágy andJackson 2000, Chabbi et al 2009). The rate of organic matter (OM) decomposition is slower in subsoil, largely due to lower concentrations of OM and lower microbial abundance, as well as increased persistence of OC due to burial, aggregation, and organomineral associations (Chabbi et al 2009, Rumpel and Kögel-Knabner 2011, Schrumpf et al 2013. Despite its persistence, OC in subsoils can be vulnerable to mineralization due to changes in climate (Jobbágy and Jackson 2000, Melillo et al 2002, Davidson and Janssens 2006, Berhe et al 2008, Chabbi et al 2009, including warming that can stimulate biogeochemical reaction rates (Hicks Pries et al 2017, Soong et al 2021). While many studies have demonstrated the relevance of subsoil OC dynamics, only a handful of studies have documented OC below the solum (A, E, and B horizons) into the lower reaches of soil commonly referred to as parent material.
We investigated the distribution, stability, and chemical composition of OM in regolith to the depth of hard bedrock, which at our study sites includes the mineral soil and saprock. Saprock is slightly weathered bedrock, where the original rock fabric is maintained, but has become more porous and friable (Graham et al 2010). We sampled to depths spanning 1.0-10.7 m along a 2300 m bioclimosequence located on the western slope of the Sierra Nevada in California (location of the Southern Sierra Nevada Critical Zone Observatory; figure 1and table S2 (available online at stacks.iop.org/ERL/16/ 124059/mmedia)). The objectives of this study were to: (a) quantify OC stocks in soil and saprock across a bioclimatic gradient; (b) characterize the age and relative degree of decomposition of OM in soil and saprock across the bioclimatic gradient; and (c) identify regional trends in deep OC storage as influenced by climate, vegetation, and degree of weathering.

Site descriptions
The four study sites were located along an elevation transect (i.e. bio-climatic gradient), from 405 m to 2700 m, on the western slope of the southern Sierra Nevada, California (figure 1). As elevation increases, mean annual air temperature (MAT) decreases, while mean annual precipitation (MAP) and deep-water percolation (DWP, calculated as the difference between precipitation and evapotranspiration (ET) (Goulden and Bales 2014)) increases. Soils include Entisols, Inceptisols and Alfisols, representing a range of soil properties. Soils reflecting the highest degree of development (as indicated by clay content and pedogenic iron) are found at mid elevations (Tian et al 2019). All sites are underlain by granodiorite bedrock, except the lowest elevation site that is tonalite, a slightly more felsic rock (Dahlgren et al 1997). Because we were unable to constrain the exact age of soil-landscapes, we assumed erosion and weathering rates have reached steady state with local climate conditions, as was the case with other studies located at and near our sites (Dahlgren et al 1997, Dixon et al 2009. Vegetation along the bio-climosequence varies systematically with climate, consisting of oak savannah at the lowest elevation site (405 m asl elevation), pine-oak forest (1160 m), mixed-conifer forest (2015 m), and subalpine forest (2700 m). Vegetation changes with elevation in response to changes in climate (figure 1). Gross primary productivity (GPP) and ET reflect vegetation patterns, where ET and GPP are low at oak savannah (350 mm yr −1 , 480 g C m −2 yr −1 ), and high at pine-oak forest (600 mm yr −1 , 1900 g C m −2 yr −1 ) and mixedconifer forests (550 mm yr −1 , 1700 g C m −2 yr −1 ), decreasing at subalpine forest (450 mm yr −1 , 700 g C m −2 yr −1 (Goulden et al 2012)) (figure 1).

Regolith sampling
A total of 36 regolith profiles (soil + saprock) were analyzed (table S2) from soil pits, Geoprobe cores, and hand augers. Hand augering (with a depth limitation of 7.6 m) was used at a subset of sites that could not be accessed by Geoprobe due to terrain. Geoprobe and hand auger samples at adjacent locations were compared to ensure no sampling effect on OC distributions.
Soil pits were sampled by genetic horizon to a depth of hard bedrock or the maximum excavation depth of 2.7 m. Geoprobe cores were collected in 1 m sections to the depth of refusal and were subsampled by 20-50 cm intervals based on color change. Hand-auger samples were collected to a depth of hard bedrock or 7.5 m, whichever was shallower, and subsampled in 20-50 cm intervals based on color. Regolith was partitioned into soil and saprock based on friability and brittleness (Graham et al 2010).

Carbon and nitrogen analyses
Soil and saprock materials were air-dried, gently crushed with mallet, and sieved (2 mm). The coarse fraction (i.e. >2 mm dia) was assumed to have no OC. Total C and N concentrations (%) of the fineearth fraction (i.e. ⩽2 mm dia) were determined on oven-dried (70 • C for 48 h) samples ground to pass a 180 µm sieve and analyzed by dry combustion (Costech Analytical ECS 4010 instrument, Costech Analytical Technologies, Inc., Valencia, CA). The average signal detected from samples of low concentration (354 ± 162 (sd), N = 10) was significantly larger than the blank (15 ± 14 (sd), N = 9). Most samples had pH values below six; however, a subset of samples was pretreated with 1 M hydrochloric acid to verify that the samples had no inorganic C. In all cases, the acid pretreatment did not result in a reduction in total C so we assumed no inorganic C was present. Hence, total C was equivalent to OC in our samples. Bulk density of soil and saprock was measured using the core method (Blake 1965). Multiple depth profiles (see table S2 for data) of OC were aggregated over 25 cm depth intervals to create an aggregated OC depth profile (figure 2) using the R package Algorithms for Quantitative Pedology (Beaudette et al 2013). The algorithm summarizes OC along depth slices from a collection of regolith cores (figure S2). Variability of OC by depth was also calculated as the 1st and 3rd quartiles.
Equations (1)-(3) describe calculations of OC stock and OC density for soil, saprock, and regolith. Thickness (d c (cm)) was calculated using the actual thickness of horizons (d (cm)) and correcting for the volumetric percent of the coarse fraction (r (%); equation (1)) OC stock (C s (kg C m −2 )) was calculated as the sum of the product of OC concentration (C (g C kg −1 )), layer thickness (d c (cm)), and bulk density (ρ b (Mg m −3 ); equation (2)) OC density (C d (kg C m −3 )) of both soil and saprock was calculated using the respective OC stock (C s (kg C m −2 )) divided by its associated (i.e. soil or saprock) thickness (d (cm); equation (3)) Polynomial and linear regression models were performed to explore correlations between climatic parameters (temperature, precipitation, and deep percolation) and C s and C d for each core. The square of the correlation coefficient (adjusted R 2 ) was calculated for soil (N = 55) and saprock (N = 22), and we evaluated statistical significance by p-value < 0.001.

Diffuse reflectance infrared Fourier transform (DRIFT) spectroscopic analyses
Soil and saprock were analyzed for OM chemical composition via DRIFT. Mid-IR spectra (400-4000 cm −1 ) were recorded using an IR spectrophotometer (Bruker IFS 66 v s −1 , Ettlingen, Germany) and a Praying Mantis DRIFT sampler (Harrick Scientific Corporation, Ossining, NY). Before analysis, all samples were ground to a powder using a ball mill (SPEX Sample Prep Mixer Mill 8000C, Metuchen, NJ, USA) and dried at 60 • C for 24 h. Potassium bromide (Aldrich Chemical Company, Saint Louis, MO, USA, FT-IR grade) was used as a background reference. We used 500 background scans and 500 sample scans with a 4 cm −1 resolution. All spectra were tangentially baseline corrected. Peak areas were calculated that corresponded to local baselines to

Radiocarbon analysis
Radiocarbon analyses were conducted after sealedtube combustion of OC to CO 2 (with CuO and Ag) that was then reduced onto Fe powder in the presence of H 2 (Vogel et al 1984). Radiocarbon values were measured in 2018 on the model FN van de Graff tandem accelerator mass spectrometer at the Center for Accelerator Mass Spectrometry (Lawrence Livermore National Laboratory, Livermore, CA, USA). OM δ 13 C values were determined using a DELTA V Plus isotope ratio mass spectrometer (Thermo Fisher Scientific, Inc., Waltham, MA, USA, Stable Isotope Ecosystem Laboratory, UC Merced). Radiocarbon isotopic values were corrected for mass-dependent fractionation with measured δ 13 C values, and were reported in ∆notation corrected for 14 C decay since 1950 (Stuiver and Polach 1977).

DRIFT-FTIR mineral content testing
DRIFT-FTIR was analyzed on regolith samples before and after OM removal (ashed at 430 • C for 8 h (Ellerbrock and Gerke 2004)) to ensure there was no mineral interference. Bulk samples and post-ashing samples were almost identical, verifying that even though regolith was low in OM, peaks of interest in the FTIR spectra represent OC in these samples.

Depth trends in OC
Depth trends in OC concentration were similar across the bio-climosequence, decreasing exponentially from soil surface horizons to uniformly low values in saprock (figure 2). The boundary between soil and saprock was between 1.7 and 2 m at all sites, except subalpine forest, where saprock did not exist (i.e. soil thickness = regolith thickness) (figure 2). Considering all sites, average OC concentration decreased from 26 g kg −1 in topsoil (A horizon) to 0.2 g kg −1 in deep saprock. The OC concentration in saprock was lowest at oak savannah. Pine-oak forest and mixed-conifer forests had similar saprock OC concentrations (figure 2).

Regional partitioning of OC stocks in soil and saprock
Soil organic C stock increased with elevation from 4.1 ± 0.6 kg C m −2 in the oak savannah site (400 m asl elevation) to a maximum of 20.2 ± 12.8 kg C m −2 in the mixed-conifer forest (2015 m; figure 3). At the subalpine forest site (2700 m), soil stocks were lower (8.1 ± 4.1 kg C m −2 ) and similar to levels at the pineoak forest site (9.6 ± 3.2 kg C m −2 , 1160 m; figure 3 and table S4). While soils exhibited the largest OC stock at every site, saprock accounted for over 25% of OC stock at mid-elevations (figure 3).
Small but detectable OC concentrations in saprock became significant contributions to overall regolith stocks (soil + saprock) when integrated over large thicknesses (figure 3). OC stock in saprock was greatest at mid-elevations (5-10 kg C m −2 ) where mean thickness ranged from 2.25 to over 9 m (figure 1). The relative proportion of saprock OC stock in regolith increased with elevation, from oak savannah (8.9%), to pine-oak forest (22.5%), and mixed-conifer forest (29.2%).

Characterization of OM in soil and saprock
Soil ∆ 14 C generally decreased with depth, but fluctuated in saprock indicating the co-existence of older and newer C ( figure 4(a)). Across all sites, ∆ 14 C values of topsoil (A horizons) were relatively less depleted (∆ 14 C −70‰-120‰), indicative of soil OC dominantly composed of newly fixed OC. The most depleted ∆ 14 C value (−921‰ or 20 300 years BP average radiocarbon age) was found in saprock at 6 m depth at mixed-conifer forest, while younger OC was found below at 9.5 m depth (−456‰ or 4800 years BP average radiocarbon age). While deep saprock ∆ 14 C values suggest, on average, the presence of relatively old OC, fluctuations between more and less depleted ∆ 14 C values with depth imply a mixture of older and younger OC. Similar fluctuations with depth in saprock were observed at pine-oak forest and oak savannah sites.
C:N ratios tended to decrease with depth in soil, but fluctuated more dramatically throughout saprock ( figure 4(b)). Fluctuations in saprock C:N generally coincided with that of ∆ 14 C values and FTIR indices of decomposition (figures 4(c)-(f)).
Ratios of peak areas derived from DRIFT-FTIR spectroscopy were used as biochemical proxies for decomposition-related processing of OM (Kaiser et al 2012, Hall et al 2018. Interestingly, at all sites, certain depths within saprock show similar ratios to soil, whereas the signature from other depths depart drastically from that of soil (figure 4). This suggests that OM in certain locations in the saprock profile, which tended to be older than soil OM (more negative ∆ 14 C), is of varying degrees of decomposition.
Although FTIR indices of OM decomposition fluctuated with depth at all sites, certain changes in the depth signature coincided with changes in ∆ 14 C. For example, mixed-conifer forest had a large peak in the ratio of aromatic C=C:carboxylate COO− (figure 3(c)) at the soil-saprock transition (∼1.75 m). This same site had a decrease in all FTIR indices at depths between 4.5 and 5.75 m (figures 4(c)-(f)), which coincided with a zone where the most depleated ∆ 14 C values were measured (figure 3(a)). At pine-oak forest, three of the four FTIR proxies showed an increase with depth (figures 4(d)-(f)) between 200 and 400 cm. This increase coincided with a decrease in ∆ 14 C values from −100‰ to −300‰. Similar coincident changes between FTIR and ∆ 14 C values occurred at the base of the regolith. At the base of regolith (4.75 m) at pine-oak forest, ∆ 14 C values were more depleated relative to the overlying layer (from −300‰ to 500‰), and FTIR indices decreased (figures 4(e)-(f)). Soil OM signatures in oak savannah were distinct from other sites, likely indicating the difference between grassland and forest vegetation on OM chemical composition, distribution, and transformations in soil.

Climate-plant-regolith interactions
Bioclimatic parameters MAT, MAP, and DWP showed varying degrees of correlation with OC stock in soil and saprock (figures 5(a)-(c)). Trends show peak storage at mid elevations (pine-oak and mixedconifer forests) where water and temperature are optimal for productivity. In saprock, correlations between OC stocks and bio-climatic parameters were stronger, showing positive linear correlations with DWP and MAP and negative correlation with MAT. Trends in OC density mirrored that of soil stocks, but not saprock (figures 5(d)-(f)). With a slope near zero, OC density in saprock is insensitive to bioclimatic conditions of this region and reflects the uniformly low concentration of OC throughout its thickness and across sites where saprock is present. It indirectly demonstrates that thickness drives OC stock trends in saprock. These bio-climatic factors, however, drive regolith thickening along the gradient.

Regional trends and controls on deep OM
Few studies have documented deep OC in the critical zone from the topsoil to bedrock. The range in OC concentration in saprock (0.15-0.94 g kg −1 ) was consistent with the few studies previously published (0.3 and 0.4 g kg −1 ) in Australia (depths of 5-38 m) and the eastern USA (depth of 6 m) (Richter andMarkewitz 1995, Harper andTibbett 2013). No studies to our knowledge have documented systematic trends in deep OC stocks in response to state factors (i.e. climate and biota).
Climate indirectly controls OC in regolith along the western slope of the southern Sierra Nevada through its impact on vegetation and weathering (figures 1 and 5). Climate exerts strong controls on saprock weathering primarily through thickening and subtle transformations of primary minerals (Dahlgren et al 1997, Tian et al 2019. Trends in regolith stocks are also influenced by GPP, which was highest at mid-elevations, contributing to high C inputs, and perhaps, to regolith thickening through biophysical weathering processes (figure 1) (Hinsinger et al 2001, Goulden et al 2012, Kelly and Goulden 2016, O'Geen et al 2018. Plant productivity and regolith thickening are limited by low precipitation in the oak savannah (lowest elevation), and low temperature in the subalpine forest (highest elevation) (figure 1). Historic glaciation has scoured away saprock at subalpine forest and, combined with the shorter growing season, soils remain shallow overlying hard bedrock (Giger andSchmitt 1993, Dahlgren et al 1997). This analysis highlights the biotic and climatic interactions that appear to control C storage for this region (figure 5). For example, precipitation and temperature directly influence OC inputs by governing GPP (figure 1) (Goulden and Bales 2014). The extent of DWP allows for the deep penetration of roots and translocation of DOM (Bales et al 2011). All three parameters influence regolith thickening, which ultimately controls saprock C stocks through preservation at depth and increased storage capacity.

Factors influencing depth distributions of deep OC
In most studies and especially in global soil C modeling, radiocarbon age increases with depth to 1 m (Shi et al 2020). The greater depths explored in our study displays a fluctuating trend in radiocarbon age within saprock, indicating that residence time of OC varies with depth. This illustrates that current soil carbon models that predict an exponential increase in radiocarbon age with depth and a corresponding decrease in C% to zero may be inaccurate. Although fluctuations in FTIR indices of OM decomposition coincided with major shifts in OC age at specific depths, correlations between age and degree of decomposition were not consistent (figure 4). This irregularity with depth is likely a result of varying degrees of preservation of old and more modern OM and variability in the vertical distribution of C inputs. Moreover, soils have a mixture of carbon compounds that are heterogeneous with respect to residence time (Sierra et al 2017). Hence, it is important to note that bulk soil radiocarbon values reflect the mean residence time across this diversity of soil carbon compounds.
The OM may be derived from roots and mycorrhizae that are known to explore deep regolith in the Sierras (Stone and Kalisz 1991, Hubbert et al 2001, Allen 2007. Moreover, sharp fluctuations in OM, have been found to indicate the presence of plant-derived OM, such as roots (Ryals et al 2014, Hall et al 2018. Changes in C:N, like those observed here in deep saprock, can indicate different types of organic materials (various residues, dissolved organic matter (DOM), plant types), whereas low C:N can indicate zones where OM is more microbially processed (Rumpel andKögel-Knabner 2011, Batjes 2016). In fact, a sagebrush steppe existed at midelevations in the Sierra Nevada during the last glacial maximum, and a transition from steppe to forest occurred about 10 000 years ago (Anderson 1990, 1996, Woolfenden 1996. These shifts in vegetation coincide with radiocarbon ages estimated to be over 20 000 years at 6 m depth at the mid-elevations ( figure 4(a)). Considering changes in vegetation history and the deep exploration of contemporary forest and shrub vegetation, it is plausible that variation in deep OC age is a result of many generations of root exploration.
Recent OM can also come from the continuous downward flux of DOM, via deep percolation, facilitated by coarse-textured soils and episodic flushing of snowmelt (Kalbitz et al 2000, Bales et al 2011, Kaiser and Kalbitz 2012, Gross and Harrison 2019. Soils in this region contain low concentrations of reactive and short-range order minerals that sorb DOM (Rasmussen et al 2005), and low sorption favors deep DOM transport (Hemingway et al 2019). Organomineral stabilization of ancient DOM is less likely given the low degree of mineral transformation in saprock (Tian et al 2019), but possible if mineral transformations are localized around roots exploiting fractures (Frazier and Graham 2000). Although, DOM can be transformed and degraded with depth (Roth et al 2019), little is known about its fate in saprock.
The existence of both old and new OC in various stages of decomposition (figure 4) is likely a result of dampened microbial activity or microbial isolation Zech 2000, Shen et al 2015). There is a large difference in OM concentrations between saprock and topsoil, and energy starvation of microbes and microbial response to changes in resource availability precludes their ability to produce enzymes that degrade OM (Allison et al 2007, Dove et al 2020. In addition, microbial activity at depth may be slowed due to the temperature effect on biogeochemical reactions being dampened by cool, stable temperatures relative to soils, which experience dynamic swings in temperature and moisture (Tian et al 2019).
In some lithologies, presence of C in deep regolith may suggest contribution of geogenic C to saprock (Zethof et al 2019). However, multiple lines of evidence suggest that this is not the case at our sites. Mineralogical analysis of saprock, through quantitative x-ray diffraction, showed no evidence of calcite (Tian et al 2019), as did acid pretreatments before and after C analysis of a subset of saprock samples. In addition, it is highly unlikely that igneous bedrock would contain geogenic OC by virtue of the nature of the type of rock formation.

Implications for global warming
Given the current trajectory of contemporary climate change, atmospheric temperatures are expected to increase by up to 4 • C by the end of the century, warming the land and underlying soil (Hartmann et al 2013). There is considerable uncertainty on how much and how fast deep regolith will warm relative to the atmosphere. One study reported that soils down to 3 m have warmed by as much as 2 • C since 1961 in China (Zhang et al 2016). At our study sites, measurements and modeling at the pine-oak and mixedconifer forest sites showed that, with depth, there is a dampening of diurnal and seasonal variability in temperature (Tian et al 2019). Soil temperature at 10 m was uniform and relatively low compared to 0.5 m, which is seasonally variable (Tian et al 2019 (Min et al 2020). One study estimated that with a 4 • C increase in temperature, deep soils have the potential to release C to the atmosphere at a rate equivalent to 30% of current fossil fuel emissions (Hicks Pries et al 2017). Regression analysis shows that MAT explained more variability in OC stock than other climatic variables tested (figure 5), suggesting that OC storage in soil may be vulnerable to warming. However, the relationship between saprock OC density and MAT suggests deep OC may be less sensitive to warming. Even though temperature warming in the next decades may affect GPP and OC decomposition in soils, there is a need to understand the vulnerability of OC in deep soils and saprock. Our findings of large OC stocks, comprised throughout the depth column as a mixture old and seemingly young, and of varying degrees of decomposition, indicates that C dynamics operate differently in deep saprock than in soil. Deep OM characteristics suggests a significant degree of C protection in saprock despite its minimal degree of mineral transformations (Tian et al 2019).

Conclusions
We discovered large OC stocks in saprock representing a previously unexplored pool of OC. Regolith thickness directly influences OC stock, which is governed in part by climate and vegetation. Overall, we observed significant amounts of deep OC storage at mid-elevations in the Sierra Nevada, where the climate is most conducive to weathering (saprock thickening) and GPP (OM inputs to regolith). The magnitude of this stock relative to soil suggests high uncertainty in the global terrestrial C inventory, especially in areas with climate that is conducive to thick regolith.
Over the next century, subsoil is projected to warm at roughly the same amount as topsoil (IPCC 2014). The effect of a temperature increase on C cycling is less certain in deep regolith than in surficial soils, where temperature fluctuations are uniformly low and OC is isolated from high rates of microbial activity. Most importantly, our findings suggest that shallow soil sampling provides a biased estimate of OC storage and leads to misleading conclusions about deep OC dynamics in ecosystems that form thick saprock layers. Our study provides compelling evidence to justify further investigations of OC in different ecosystems with deep regolith to fully understand the how one of the largest pools of OC will respond to changes in climate.

Data availability statement
All data are available in the manuscript, the supplementary materials; datasets have been archived on Figshare  The data that support the findings of this study are openly available at the following URL/DOI: https:// doi.org/10.6084/m9.figshare.13106102.