Carbon Dioxide Production in Bedrock Beneath Soils Substantially Contributes to Forest Carbon Cycling

Soils are widely considered the primary terrestrial organic matter pool mediating carbon transactions with the atmosphere and groundwater. Because soils are both a host and a product of rhizosphere activity, they are thought to mark the location where photosynthetic fixation of carbon dioxide (CO2) is balanced by the oxidation of organic matter. However, in many terrestrial environments, the rhizosphere extends below soils and into fractured bedrock, and it is unknown if the resulting biological and hydrologic dynamics in bedrock have a significant impact on carbon cycling. Here we show substantial production of CO2 in weathered bedrock at 4–8 m below the thin soils (<0.5 m thick) of a Northern California forest using innovative monitoring technology for sampling gases and water in fractured rock. The deep CO2 production supports a persistent upward flux of CO2(g) year‐round from bedrock to soil, constituting between 2% and 29% of the average daily CO2 efflux from the land surface. When water is rapidly transported across the fractured bedrock vadose zone, nearly 50% of the CO2 produced in the bedrock dissolves into water, promoting water‐rock interaction and export of dissolved inorganic carbon (DIC) from the unsaturated zone to groundwater, constituting as much as 80% of the DIC exiting the hillslope. Such CO2 production in weathered bedrock is subject to unique moisture, temperature, biological, and mineralogical conditions which are currently missing from terrestrial carbon models.


Introduction
The conversion of photosynthetically derived organic carbon to CO 2 makes soils, and the suite of autotrophic and heterotrophic respiratory pathways that they host, a principal component of the terrestrial carbon cycle (Wenzel et al., 2016). Future climate change could increase the flux of CO 2 from soils to the atmosphere in some locations (Melillo et al., 2017). A principal source of uncertainty in current climate forecasts is how the global balance between net primary productivity and respiration in soils, and therefore the functioning of the terrestrial ecosystem as a source or sink of CO 2 , will shift (Bradford et al., 2016;Melillo et al., 2017;Wenzel et al., 2016). New carbon modeling paradigms are emerging in response to increasing recognition of CO 2 sources and transport pathways which operate separately from or in tandem with the basic balance between photosynthesis and soil respiration (Blankinship et al., 2018), questioning our understanding of subsurface carbon storage (Bradford et al., 2016).
Significant stores of organic carbon beneath depths typically considered in Earth system models (Harper & Tibbett, 2013) are increasingly recognized as potentially labile (Gross & Harrison, 2019;Keller & Bacon, 1998). This includes petrogenic organic carbon, or organic carbon hosted in the bedrock, that is present in shales exposed over a quarter of the terrestrial surface (Amiotte Suchet et al., 2003), and which can be a source of carbon for CO 2 production in deep vadose zones (Keller & Bacon, 1998;Longbottom & Hockaday, 2019;Wood & Petraitis, 1984). In addition to the degradation of petrogenic organic carbon, root networks commonly penetrate the bedrock beneath soils (Schwinning, 2010). Deep roots are known to alter water dynamics at depth (e.g., Fan et al., 2017;Rempe & Dietrich, 2018), and a recent study of a humid hardwood forest demonstrated that the respiration rate of roots extracted from bedrock fractures (upper 1.5 m) was comparable to roots extracted from the soil (Hasenmueller et al., 2017). The recycling of root tissue and the support of the microbiome surrounding root structures are all a product of deep rooting and contribute to a supply of organic carbon that can be respired at depths which are typically excluded from field samples and large-scale models (Gross & Harrison, 2019;Richter & Billings, 2015). While it may seem intuitive that these deep roots would drive a dynamic carbon cycle at depth, the relative importance of root carbon cycling at depth has not been quantified.
In addition to the respiration of deep, labile carbon pools, deep vadose zones are also thought to play an important role in carbon cycling through weathering (e.g., Winnick et al., 2017). The solubilization and aqueous transport of respired CO 2 within weathered bedrock (Angert et al., 2015;Sánchez-Cañete et al., 2018) can facilitate CO 2 -driven mineral weathering reactions such as carbonate weathering (Bergel et al., 2017) and silicate weathering (Brantley et al., 2014). The enhancement of air-filled pore space generated through CO 2 -driven weathering may also promote oxygen delivery to depth. These advance sulfide weathering fronts which decrease alkalinity, may promote the dissolution of carbonates and result in a source of CO 2 (Winnick et al., 2017). It is unclear if the participation of respired CO 2 in these weathering reactions within the unsaturated zone results in a net source or sink of CO 2 from the ecosystem scale to the global scale (Carroll et al., 2018;Torres et al., 2017;Winnick et al., 2017). However, it is clear that the role of "deep" carbon cycles is now emerging as a potential missing feature of the carbon budget (Gross & Harrison, 2019;Harper & Tibbett, 2013;Rumpel & Kögel-Knabner, 2011;Wan et al., 2018).
Such new perspectives extend conceptual frameworks deeper into the subsurface and these observations imply that deep carbon cycling must occur within weathered bedrock, but as yet the fluxes and rates of these pathways have not been directly constrained. The objective of this study is to establish the potential for deep carbon cycling to occur within weathered bedrock and directly measure and calculate fluxes of inorganic carbon from the unsaturated zone in weathered bedrock. We use discrete gas and water monitoring within a 16 m thick weathering profile underlying a steep, forested hillslope to document spatiotemporal patterns of CO 2 production and the fate of this CO 2.

Site Description
The field site is a steep (average slope of 30 o ) intensively monitored, north facing hillslope within the old-growth mixed broadleaf needleleaf evergreen forests of the Northern California Coast Ranges. The instrumented hillslope (Figure 1) is associated with the Eel River Critical Zone Observatory (ERCZO). The bedrock is mapped as argillaceous marine turbidites of the Coastal Belt of the Franciscan Complex that, in some places, exhibit 90 o dip (McLaughlin et al., 2000). The region experiences active uplift and erosion rates of about 0.2-0.4 mm/yr (Fuller et al., 2009). The dense argillaceous parent bedrock is exposed in the actively incising channel at the base of the hillslope and is poorly fractured and dark gray in color where unweathered. Parent bedrock organic matter content is 0.07-0.88% (Gu et al., 2020) and the mineralogy is a combination of primary minerals and secondary minerals including quartz (24.8%), albite (21.9%), chlorite (18.3%), smectite (12.6%), illite (74%), kaolinite (3.5%), and carbonates (2.2%) (Gu et al., 2020).
Boreholes distributed across the study hillslope reveal that the hillslope is primarily argillite with minor sandstone interbeds. Soils, defined as the unconsolidated, geomorphically mobile layer, are typically between 30 and 70 cm thick. Underlying the mobile soil, the bedrock maintains its relic structure but is visibly weathered (Figure 2, inset). The weathered bedrock is characterized by blocks of fine-grained dark gray argillaceous bedrock bounded by fractures of variable frequency and aperture width that are oxidized along fracture surfaces. The transition between weathered and unweathered bedrock is marked by a shift in fracture density, mechanical strength, iron oxidation along fracture surfaces, and loss of pyrite (Gu et al., 2020;Rempe & Dietrich, 2018). The depth to unweathered bedrock increases from~4 m near the channel to about 30 m at the ridgetop (Rempe & Dietrich, 2018). Sampling ports and sensors traverse the bedrock weathering profile via a vadose-zone monitoring system (VMS). The upper 12 m remains unsaturated year-round and is the focus of gas and pore water sampling. Sampling of soil pore water and groundwater is accomplished via lysimeter and groundwater well, respectively. Photo (right) illustrates the thin soils and fractured bedrock that are characteristic of the upper 1 m. below this depth, drilling, roadcuts, and exposures in the stream reveal variably weathered and fractured argillaceous marine turbidites of the Franciscan complex. Hydrologic monitoring at the site via groundwater monitoring, moisture sensors, and neutron probe logging reveal that the water table seasonally fluctuates at the base of the weathering profile (Rempe & Dietrich, 2018;Salve et al., 2012). Above the water table, the soils and weathered bedrock remain unsaturated year-round. In the Mediterranean climate of the site, most rainfall occurs between January and April leading to seasonal wetting and drying of the unsaturated zone. Runoff to the stream channel is recharged by water transiting the meters-thick bedrock vadose zone. Recharge occurs only after the water content throughout the unsaturated soil and weathered bedrock reaches a field capacity-like state, which takes~500 mm of precipitation. During the summer dry season, deeply rooted trees withdraw water from the weathered bedrock zone to depths of up to 12 m (Rempe & Dietrich, 2018).
At the vadose zone monitoring system (VSM) location ( Figure 1) where our monitoring is concentrated, soils are 70 cm thick and the transition from weathered to unweathered bedrock is~13.5 m. The water table fluctuates between~8 and 14 m and the upper 5 m experiences significant swings in water content associated with plant water uptake (Rempe & Dietrich, 2018). There is virtually no understory at the site and the trees are a mix of Douglas-fir (Pseudotsuga menziesii), Pacific madrone (Arbutus menziesii), tanoak (Notholithocarpus densiora), and live oaks (Quercus wislizenii and chrysolepis) (Oshun et al., 2015).

The VMS
We use a VMS to discretely sample water and gas within the bedrock weathering profile at~1.5 m intervals ( Figure 2). The VMS is composed of sensors and samplers installed within two inclined boreholes that traverse the length of the unsaturated zone, to~16 m ( Figure 2). Flexible sleeves that house sampling cells are placed~1.5 m apart (Rimon et al., 2007) and were installed within the inclined boreholes. One sleeve contains 0.25 m long ports designed to capture water held under tension within the unsaturated zone. This sleeve also houses flexible time-domain-transmission (TDT) sensors which measure temperature and water content. A second sleeve contains 1.5 m long sampling ports designed to collect freely draining water traversing the unsaturated zone. Perforated tubing is embedded within these ports for gas sampling. Table S1 in the supporting information reports the depths of gas and water samples and temperature and water content measurements in the VMS. Ports 1-9 are within weathered bedrock and Port 10 is located within unweathered bedrock. Ports 1-4 are unsaturated year-round and Port 10 is saturated year-round, while Ports 5-9 are located within the water table fluctuation zone and experience seasonal saturation. Sample tubing and sensor cables are all accessible from a centralized control panel. Technical details regarding the operation of the VMS are described in a number of studies including Dahan et al. (2003) and Rimon et al. (2007). In addition to the 10 VMS ports, a static flux chamber was installed in 2018 to measure gaseous concentrations at the base of the soil and top of the weathered bedrock at 0.8 m.

CO 2 and O 2 Concentration
CO 2 and O 2 concentrations in the gas phase (CO 2(g) or P CO2 , O 2(g) or P O2 ) were measured in the field using a Quantek Instruments Oxygen and Carbon Dioxide Analyzer (Model 902P) by pumping gas at~300 cc/min through VMS gas sampling ports. The static flux chamber at 0.8 m was used as an additional sampling port to increase sampling resolution at the soil-weathered bedrock boundary. The CO 2 sensor was calibrated periodically throughout the study by Quantek instruments and the O 2 sensor was calibrated before and during each sampling campaign. Each sensor has an error of 0.2% of the measured value. Between 18 June 2018 and 9 February 2019, an in-line sampling port was used to collect gas samples in 30 ml and 60 ml serum vials which were analyzed by gas chromatography (GC) at the UT Jackson School of Geosciences. A calibration curve between field (Quantek) and lab (GC) CO 2 was created and applied to the field data where it was noted in the data spreadsheet to account for drift in the sensor.

CO 2(g) Flux
Concentrations of CO 2 were used to calculate CO 2 flux according to the flux-gradient method outlined in Jury et al. (1991) using an adaptation of Fick's law appropriate for use in porous media: is the flux of CO 2 , dCO 2(g) /dz is the change in concentration over vertical distance (in this study 1 to 4 m) and D s is the diffusivity of CO 2(g) in porous media. The diffusion coefficient of CO 2 was constrained using Penman's equation (Penman, 1940), D s = D a 0.66θ a where θ a is the air-filled porosity and D a is the diffusivity of CO 2. Following Massman (1998), we vary through D a for each sampling location and time to account for temperature differences. Although we do not know the exact θ α for each sample time, we use a range of 5-10% θ a to calculate a range of D s for each sample date to encapsulate both the wet and dry seasons. θ a is lowest when water contents are high (January to June) and highest when water contents are low (July to December). θ a is equal to the year-round air-filled porosity, which we estimate as 5%, plus the changing water-filled porosity, that is, water content, which varies by~5% in weathered bedrock throughout the site (Rempe & Dietrich, 2018). Gu et al. (2020) performed a porosity analysis on core from Well 10 at the ERCZO which is adjacent to the VMS ( Figure 1) and calculated total weathered bedrock porosities in the upper 5 m approaching 15% with~5% water-accessible porosity, consistent with the observations by Rempe and Dietrich (2018). We estimate that 5% remains air-filled year round due to the extensive fracturing and microfracturing at the site, in part evidenced by the~5% epoxy-accessible porosity measured in rock chips alone by Gu et al. (2020).
A Licor 6,400 with soil chamber attachment was used to measure flux from the surface of the soil in collars installed near the VMS in 2018. A minimum of three measurements were made at each collar ( Figure 1).

Apparent Respiratory Quotient (ARQ) Analysis
Concentration gradients of CO 2(g) and O 2(g) from 1 to 4 m were used to calculate a flux of CO 2 based on an apparent respiratory quotient (F ARQ ) (Angert et al., 2015). This value reflects the expected CO 2(g) flux if the observed O 2(g) gradient is attributed to oxidation of organic carbon and corrected for the difference in diffusivities of O 2 and CO 2 . F ARQ is therefore an estimate of respiration, or R CO2 . Departures from a range of ARQ values that represent variation in organic carbon source (~0.7-0.9, Hicks Pries et al., 2020) are caused by processes that consume more CO 2 or O 2 including CO 2 dissolution, weathering of carbonate, weathering of iron and sulfur bearing minerals, and differences in ARQ between species (Angert et al., 2015). In subsurface environments, CO 2 dissolves in pore water, leading to CO 2 removal and lowering CO 2(g) concentrations. Because F b is an estimate of CO 2 production in the subsurface that is calculated using CO 2(g) concentration gradients, CO 2(g) dissolution leads to an underestimation of CO 2 production. F ARQ accounts for this discrepancy by comparing the CO 2(g) gradient to the O 2(g) gradient.
At the study hillslope, pyrite has been completely removed throughout the depths of interest to this study (Gu et al., 2020), and so we assume O 2 depletion is not due to weathering of iron and sulfur bearing minerals. To calculate F ARQ for each sampling date, we first calculated the ARQ according to Angert et al. (2015): ARQ = −0.76 (ΔCO 2(g) /ΔO 2(g) ). We make the simplifying assumption of a constant reference ARQ, ARQ ref (CO 2 respired to O 2 consumed during respiration), of 0.9, similar to the methods of Sánchez-Cañete et al. (2018). F ARQ is calculated by correcting F b using the O 2(g) gradient:

Dissolved CO 2
Aqueous samples were collected from the VMS ports (both freely draining and held under tension), filtered with a 0.2 μm filter, stored in either 20 or 40 ml glass VOA vials, refrigerated, and analyzed at the UT Jackson School of Geosciences Aqueous Geochemistry Lab. Dissolved inorganic carbon (DIC) was measured using a Teledyne Tekmar Apollo 9,000 carbon analyzer by acidification and quantification of the evolved CO 2 by nondispersive IR. The pH was immediately measured in the field with an Accumet pH meter and combustion electrode using unfiltered samples. Cation and anion concentrations were measured at the University of Illinois Urbana Champaign on a Thermo Fisher Scientific iCAP-Q Inductively Coupled Plasma Mass Spectrometer (ICP-MS) and Dionex ion chromatograph, respectively.
We model the aqueous system of each sample using Phreeqc ver. 3.5.0 to account for the possibility of exsolution of CO 2 during sampling from the VMS (e.g., Suarez, 1987). CO 2 degassing, which increases pH and decreases DIC, is a well-known problem for pH and DIC measurements (Suarez, 1987) especially from lysimeters. In our simulations, we make three necessary assumptions: (1) The aqueous phase is in equilibrium with the gas phase. This is routinely assumed in the unsaturated zone (Keller & Bacon, 1998;Stinchcomb et al., 2018). (2) The solution is undersaturated with respect to calcite, prohibiting calcite precipitation in the lysimeter during sample collection. This is consistent with the observation that the upper eight ports (above 13 m) of the VMS are always significantly undersaturated with respect to calcite and that carbonate minerals are present only in small amounts (inorganic carbon varies between 0.03% and 0.07%) at these 10.1029/2020JG005795

Journal of Geophysical Research: Biogeosciences
depths.
(3) The system is stable over the span of 1 day, leading to an invariant CO 2(g) concentration in this short time interval. We have found that the CO 2(g) in weathered bedrock does not change over the length of 1 day within the span of instrument error.
The equilibrium modeling approach is derived from the method of Elberling and Jakobsen (2000). In our simulation, the CO 2(g) is specified using our field measurements and Henry's Law is used to partition CO 2 to the aqueous phase. A set of equilibrium relationships are iteratively solved until pH and DIC convergence, maintaining charge and proton balance. Cation and anion concentrations, temperature, and CO 2(g) from the samples were kept constant while field pH and measured DIC were specified as initial conditions (Table 1) but were allowed to change in order to bring the solution into equilibrium with the CO 2(g) . We applied this method to samples that were within ±5% charge balance and with an SI calcite < −0.5, to fulfill the assumption that carbonate mineral dissolution is negligible. The use of DIC instead of alkalinity is a departure from the method outlined by Elberling and Jakobsen (2000) and was necessary because of sample size limitations. To check this approach, we estimate an alkalinity from the remaining anion charge needed to reach charge balance in the aqueous sample (Table 1). We find that the charge balanced, alkalinity approach produces a pH which agrees within 1.1% of our first method. Similarly, the DIC agrees within 5.5% for all dates.
We estimate soil [DIC] at 0.3 m depth (Table 1) using the following additional assumptions: Due to the unavailability of gas samples from the soil lysimeter we use a CO 2 of 1% which falls within the concentration gradient profile (Figure 3). We use charge balance with respect to the cation and anion data, attributing imbalance in our measured values to bicarbonate. We then calculate the soil DIC as CO 2(aq) + [HCO 3 − ]. We assume that [HCO 3 − ] makes up the charge balance deficit. Figure 3 illustrates the seasonally dynamic CO 2(g) and O 2(g) concentrations within the weathering bedrock for the period 2017-2019. Figure 4 shows CO 2 concentrations, O 2 concentrations, temperature, and relative rock moisture in weathered bedrock. Relative rock moisture is the amount of rock moisture relative to the maximum observed water content. It varies from 0 to 1 where 0 represents the lowest observed water content and 1 represents the highest observed water content. During the wet season, water contents are high relative to the dry season (Figure 4c). When water contents are high in the wet season (blue, Figure 4b), O 2(g) is low relative to the dry season (orange, Figure 4b). CO 2 is less seasonally dynamic than O 2 . Nonetheless, wet season CO 2 is generally higher than dry season CO 2 in the upper 5 m of the profile (Figures 4a and S1). Below 5 m, the opposite is true so that CO 2(g) is lower in the wet season relative to the dry season. CO 2 concentrations remain high in the summer even as water contents Note. The sample on 3/8 at 0.3 m is a soil lysimeter sample and the samples at 10.7 and 12.2 m were taken from the VMS.

Journal of Geophysical Research: Biogeosciences
decrease and air-filled porosity increases ( Figure S1). During both wet and dry seasons, CO 2(g) concentrations increase and O 2(g) concentrations decrease from the surface to~4-8 m depth (Figures 4a, 4b, and 3). When water contents from 1.5-11 m are low (July to December), this inverse concentration correlation approaches 1:1, O 2(g) :CO 2(g) at all depths (Figures 4d-4f) which reflects organic carbon oxidation (Bergel et al., 2017). When water contents from 1.5-11 m are high (January to June), the inverse concentration correlation departs from the 1:1 line, indicating other processes are altering CO 2 and O 2 concentrations such as CO 2 dissolution, iron redox processes (Kim, Stinchcomb, et al., 2017), or methane oxidation for example.
The deep respiration of CO 2 within bedrock leads to a persistent concentration gradient, and thus upward diffusive flux, from the location of CO 2 production in weathered bedrock toward the base of the soil (Figure 4a). The upward CO 2 flux from bedrock to soil, F b , is calculated using the observed concentration gradient from 1-4 m depth at different times of year (between 0.11 and 0.85 mol C/m) and a diffusion coefficient which lies between 4.9 × 10 −3 and 9.8 × 10 −3 cm 2 /s (Table S2). This diffusion coefficient spans the range of diffusivities we calculated based on 5-10% air-filled porosities (see section 2). This diffusion coefficient is also consistent with diffusivities estimated in other deep vadose zones (Tokunaga et al., 2016) and is lower than diffusivities of deep soils (Davidson & Trumbore, 1995). During the 2.5-year monitoring period, we expect air-filled porosity to be close to 10% when the subsurface is driest (October, Figure 4c)) and 5% when the subsurface is wettest (April, Figure 4c). We calculate the flux of CO 2 from weathered bedrock toward the soil, F b , for each date with diffusivities bounded by 5% and 10% air-filled porosity, with the expected F b falling between these two values depending on subsurface conditions. We find that for 5% air-filled porosity (minimum diffusivity), F b ranges from 0.06-0.43 g C/m 2 day over the study period of 2.5 years, and for 10% air-filled porosity (maximum diffusivity), F b ranges from 0.12-0.87 g C/m 2 day ( Figure 5). The average F b for each of these subsurface air-filled porosity scenarios (5% and 10%) is 0.28 and 0.55 g C/m 2 day, respectively, or 102 and 201 g C/m 2 year. Measurements of F SOIL at the ground surface at the site taken throughout the year ranged from 1.08-5.98 g C/m 2 day. The average F SOIL is 3.02 g C/m 2 day, or 1,100 g C/m 2 year. Thus, the upward flux of CO 2 from bedrock to soil, F b , amounts to between 2% and 29% of the average daily F SOIL at the site. Globally, the soil efflux, F SOIL , for Mediterranean forest soils ranges from 27-2,190 g C/m 2 year with a mean of 964 g C/m 2 year (Bond-Lamberty & Thomson, 2018). This indicates that our study site is average with respect to soil efflux and suggests that in other locations where there is substantial CO 2 production in weathered rock, F b could constitute a significant fraction of F SOIL .
In addition to the flux of CO 2 from bedrock to soil calculated from CO 2 concentrations (F b ), we calculated a flux of CO 2 based on an apparent respiratory quotient (F ARQ ) (section 2). We find that the calculated ARQ for each date ranges from 0.17 to 1.46 over the 2.5-year sampling period. ARQ is greater than 1 where ΔCO 2(g) is greater than the ΔO 2(g) which occurs in environments under anaerobic conditions. ARQ is less than 1 in environments where ΔCO 2(g) is less than the ΔO 2(g) such as those where mineral oxidation consumes O 2(g) or CO 2(g) dissolves into rapidly infiltrating water. These ARQ values are addressed further in section 4. For the F b calculated under the minimum air-filled porosity scenario (5%, wet subsurface), R CO2 (or F ARQ ) ranges from 0.09-1.2 g C/m 2 day ( Figure 5). For the maximum air-filled porosity scenario (10%, drier subsurface), R CO2 ranges from 0.19-2.4 g C/m 2 day ( Figure 5). Over the 2.5-year monitoring period, F ARQ is rarely less than F b , indicating that anaerobic redox conditions exert little influence on the concentrations of major gases at the scale of analysis (meters). When water content and water fluxes throughout the profile are at a maximum (January to May), F ARQ is large and significantly larger than F b , and when water content and water fluxes are at a minimum (July to December) F ARQ is nearly equal to F b ( Figure 5).

Discussion
While a considerable amount of CO 2(g) is transported from weathered bedrock to the soil, some fraction of the CO 2 produced at depth dissolves into the water and leads to subsequent silicate and carbonate weathering. Two lines of evidence support this inference. First, O 2(g) :CO 2(g) stoichiometry departs from the 1:1 Figure 5. Evidence for dynamic CO 2 production in weathered bedrock. Top panel shows time series of cumulative seasonal precipitation for the record of observation and the relative moisture content within the weathered bedrock at 6 m, with periods of increasing water content highlighted by blue shading. The lower panel compares the seasonal patterns of CO 2 flux estimated using two methods: CO 2 concentration gradient (i.e., F b ) (gray shading) and F b corrected for the O 2 concentration gradient (i.e., F ARQ ) (hatched shading). The difference between the two flux estimates represents the amount of CO 2 dissolved in the aqueous phase, which is highest during periods with high water content, and this high water flux, in the unsaturated zone. The range of values reported represents the variation in possible θ a with the bottom and top lines representing the flux at 5% and 10%, respectively.

Journal of Geophysical Research: Biogeosciences
relationship expected for respiration at all depths from January to June (Figures 4d-4f), with the largest deviation between 6-7.5 m (Figure 4e). This coincides with the location of maximum CO 2(g) and minimum O 2(g) in the profile. The characteristics of this O 2(g) :CO 2(g) departure are consistent with CO 2 dissolution and carbonate weathering (Bergel et al., 2017). Second, we compare F b to F ARQ . F ARQ may deviate from F b due to CO 2 dissolution, CO 2(aq) advection, abiotic O 2 consumption, differences in respiratory quotient among species, and other biotic reactions (Angert et al., 2015). We observed a range of ARQ values from 0.17 to 1.46 over the 2.5 year study. While recent studies in soils have highlighted that ARQ changes with respect to substrate and environmental variables (Hicks Pries et al., 2020), we find that, when the subsurface is wet from January to July and water fluxes are expected to be at a maximum, the majority of the ARQ values observed are low (mean 0.46, standard deviation 0.18, range 0.17 to 0.84). This is much lower than the range of ARQs expected for variable organic carbon substrates (~0.7 to 0.9, Hicks Pries et al., 2020), indicating that CO 2 dissolution and CO 2(aq) advection are likely to be the dominant influence on the stoichiometric relationship between CO 2(g) and O 2(g) . This is in agreement with ARQ observations made in soils (Angert et al., 2015;Sánchez-Cañete et al., 2018). The weathered bedrock region does experience changes in temperature and moisture content. However, the changes in these variables at depth are significantly dampened compared to the soil environment (Figure 4c and Rempe & Dietrich, 2018), further supporting our attribution of low ARQ to CO 2 dissolution and CO 2(aq) advection. We find calculated F ARQ exceeds observed F b at nearly all times of year ( Figure 5). This is consistent with dissolution and downward transport of CO 2(aq) from the location of maximum R CO2 .
The largest amount of CO 2 dissolution, reflected by the largest discrepancy between F ARQ and F b, occurs when bedrock water content and water fluxes throughout the profile are at a maximum (January to May). We estimate that as much as 81% of the produced CO 2 is dissolved into water during the wet season. The alignment of high R CO2 and CO 2(aq) with high water fluxes indicates significant transport of CO 2(aq) and likely leads to strong seasonal differences in bedrock weathering. Indeed, groundwater observations at the site indicate that CO 2 -enhanced cation exchange reactions and other fast mineral weathering processes dominate during the wet season when water is rapidly transiting the unsaturated zone (Kim et al., 2014). Table 2 reports the fraction of DIC within the unsaturated, weathered bedrock zone relative to the groundwater DIC that supplies streamflow (i.e., DIC base unsat − DIC soil )/DIC well 12 ). Our calculation assumes that DIC SOIL can be calculated by charge balance and equilibration with soil CO 2(g) . It is likely that dissolved organic carbon (DOC) is an additional anion. In terrestrial systems, the majority of the DOC anion charge derives from carboxyl groups (Ritchie & Perdue, 2003) with an average charge density of 12.8 meq/g C in  Note. DIC observed within vadose zone pore waters, groundwater, and streamflow is used to calculate the contribution of CO 2 dissolution in the bedrock vadose zone to groundwater and stream DIC. Values are reported for wet (April) and dry (September), and transition from dry to wet (December) conditions. Kim et al. (2014) report that significant DIC is lost due to degassing and calcite precipitation as groundwater supplies runoff. A significant fraction of this DIC is derived from the bedrock vadose zone, and this contribution is different in wet and dry seasons. a There was not enough water in the soil lysimeters so the same contribution from 22 April 2017 was assumed. b Depth of unsaturated zone sample at certain times of the year. c Well 12 sampled on 11 October 2017. d Elder Creek sampled on 9 September 2017.
fulvic acids (Ritchie & Perdue, 2003). While the sample collected on 8 March 2017 was not analyzed for DOC, we did measure a DOC concentration of 7.1 ppm C at 40 cm depth from a nearby lysimeter from 23 April 2017. For this sample, the resulting contribution to the charge balance equation would be 0.091 meq/L, or roughly one quarter of the calculated anion contribution. This results in a smaller contribution of DIC from the soil to the weathered bedrock (0.79 mM instead of 0.88 mM) which translates to an underestimation of the DIC contribution of weathered bedrock to groundwater by less than 2%. At other times of year, it is possible that the contribution of soil DOC to buffering capacity is larger, and the DIC SOIL contribution is even lower.
Well 12 was chosen to represent groundwater for calculations in Table 2 due to its location at the base of the hillslope (Figure 1). The groundwater is thus an integrated signature of the upslope contributing area and is likely characteristic of water exiting the hillslope and entering the stream. Of the total DIC exported from the hillslope, between 63-80% directly results from the production of CO 2 within the unsaturated, weathered bedrock zone. This finding contrasts with the assumption that weathering is focused in soils or within groundwater, instead indicating that the production and transport of CO 2 within unsaturated, weathered bedrock leads to significant weathering and export of DIC from the hillslope.
Evasion of CO 2 from headwater streams, like Elder Creek which drains the base of our hillslope, is known to be significant, highly variable in space and time, and directly linked to heterogeneity in subsurface flow paths (Duvert et al., 2018). It has been suggested that the source of this CO 2 is in part biogenic (Horgby et al., 2019). Our study shows that respiration in weathered bedrock is a significant source of dissolved CO 2 to groundwater and streamwater, and therefore a likely source of CO 2 for evasion from streams. Previous observations from the study location indicate CO 2 degassing and subsequent calcite precipitation as groundwater exits the subsurface and enters the stream (Kim, Dietrich, et al., 2017). The complex flow paths and rapid transit times of water, and thus dissolved CO 2 , through weathered and fractured hillslopes may be one explanation for this variability in CO 2 evasion in streams. Indeed, recent studies investigating CO 2 dissolution in soils have shown that pulsed wetting fronts with atmospheric CO 2(aq) concentrations are capable of rapidly dissolving CO 2 generated in the subsurface environment, initiating local silicate weathering events, and transporting dissolved CO 2 downgradient (Olshansky et al., 2019). At our site, it is likely that this occurs when water is rapidly transiting the weathered bedrock vadose zone and this mechanism is supported by our F ARQ calculations. If this is the case, then it may result in periodically lower fluxes of CO 2 to the overlying soil (F b ). The range over which temperature and water content fluctuate decreases with depth ( Figure 4c), yet O 2(g) and CO 2(g) show large amplitude seasonal changes throughout the profile (Figures 4a and 4b). While soil organic carbon respiration rates in soils are strongly regulated by both temperature (e.g., Q 10 ) and water content (e.g., Birch effects; e.g., Raich & Schlesinger, 1992;Skopp et al., 1990), our records of water content and temperature (Figure 4c) indicate that it is possible these relationships will not hold when extended to deeper, weathered, and fractured bedrock. Instead, other processes may need to be considered. For example, the availability of O 2(g) may be impacted by heterogeneous flow paths or the saturation and desaturation of pore spaces in the region where a water table fluctuates. The O 2(g) minimum that occurs during the wet season cannot be explained by 1-D diffusive transport alone, suggesting that the dissolved phase and interactions between the unsaturated and saturated conditions may play a role in O 2 availability in the unsaturated zone.
There are several possible sources of organic carbon for the respiration observed in weathered bedrock. These mechanisms could include translocation and oxidation of DOC from soils with rates that vary with depth (Keller, 2019), deep root exudation (Schwinning, 2010), or petrogenic OC oxidation (Soulet et al., 2018). During much of the dry season, there are negligible downward fluxes of water, as evidenced by the soil water tension exceeding that of the VMS lysimeter vacuum pressure by the end of June. Delivery of DOC from soils may be important during the wet season when water is transported from soil to weathered bedrock, but translocation of DOC does not explain the sustained high CO 2 concentrations observed in the weathered bedrock during the summer growing season ( Figure S1). During the dry season, when water contents are dropping and O 2 is increasing, deep root exudation may play an important role. Organic carbon in the shale itself may also be a source for respiration (Soulet et al., 2018). While difficult to separate, different mechanisms of OC delivery to the subsurface will have distinct implications for long-and short-term carbon cycling in addition to differentially controlling the production of CO 2 in weathered bedrock.

Journal of Geophysical Research: Biogeosciences
We've found coupled O 2(g) and CO 2(g) dynamics are regulated by a combination of temporally variable transport of water and gas through the interconnected fracture network and CO 2(g) dissolution. With up to 29% of the surface CO 2(g) flux and between 63-80% of the hillslope DIC exports coming from unsaturated, weathered bedrock, it is clear that understanding the drivers of CO 2 production in this region is important. A recent study called attention to the failure of terrestrial biosphere models to account for the terrestrial carbon contribution to inland waters, thus overestimating the amount of carbon that could be stored in the subsurface (Butman et al., 2016). Therefore, accurately conceptualizing and quantifying the contribution of terrestrial carbon exports to the atmosphere and to rivers and streams is critically important for understanding of global carbon cycling, particularly as our climate changes. Our work supports a conceptual model of respiration and subsequent DIC and weathering that occurs on a continuum in the subsurface (Keller, 2019), in which the weathered bedrock region, or regolith, is a key component of terrestrial carbon cycling.

Conclusions
Given that rooting in bedrock is common in a variety of lithologic, climatic, and topographic settings Fan et al., 2017;Schwinning, 2010), rock respiration associated with deep rooting could be an important, but currently missing, component of the carbon cycle. Here, we demonstrate that rock respiration is dynamic and constitutes a significant process by which CO 2 is exported as a gas flux to the atmosphere, and a dissolved flux to groundwater and river systems. We directly observe phenomenon hypothesized by previous researchers, namely that respiration within weathered bedrock directly transforms the aqueous phase at depth, leading to a deep source of acidity with the potential to increase porosity, facilitate delivery of reactive fluids, and deepen weathering fronts of carbonate and sulfur-bearing minerals (Hasenmueller et al., 2017;Kim, Dietrich, et al., 2017;Stinchcomb et al., 2018;Winnick et al., 2017). Our exploration in the region of the critical zone below soils and above groundwater indicates that representing subsurface carbon as a one box or 1-D model that only considers translational flow of carbon from surface to depth is insufficient to capture the deep carbon dynamics in bedrock rhizospheres.
To quantify the drivers of rock respiration, an understanding of the sources of organic carbon and their dynamics is needed. In this study, rooting and associated water uptake occurs within unsaturated, weathered bedrock that contains between 0.12% and 0.77% of organic matter. While it is likely that deep root respiration and associated processes such as decomposition within the rhizosphere contribute to the observed CO 2 dynamics, we have not ruled out shale-derived organic carbon as a carbon source for respiration. The translocation of DOC from soils to the bedrock profile is unlikely to play an important role, because rock respiration persists during the driest periods of the year when aqueous transport is negligible. At present, a dependence of respiration on soil moisture and temperature is broadly implemented into Earth System Models (Lehmann & Kleber, 2015), which does not account for the efflux of bedrock respired CO 2 , which may be independent of soil temperature and moisture conditions. Future study is therefore needed to identify the source of organic carbon and its relationship to environmental variables, such as water content and temperature, such that respiration rates can be estimated and the fate of CO 2 determined.