Illuminating the “Invisible”: Substantial Deep Respiration and Lateral Export of Dissolved Carbon From Beneath Soil

Dissolved organic and inorganic carbon (DOC and DIC) influence water quality, ecosystem health, and carbon cycling. Dissolved carbon species are produced by biogeochemical reactions and laterally exported to streams via distinct shallow and deep subsurface flow paths. These processes are arduous to measure and challenge the quantification of global carbon cycles. Here we ask: when, where, and how much is dissolved carbon produced in and laterally exported from the subsurface to streams? We used a catchment‐scale reactive transport model, BioRT‐HBV, with hydrometeorology and stream carbon data to illuminate the “invisible” subsurface processes at Sleepers River, a carbonate‐based catchment in Vermont, United States. Results depict a conceptual model where DOC is produced mostly in shallow soils (3.7 ± 0.6 g/m2/yr) and in summer at peak root and microbial respiration. DOC is flushed from soils to the stream (1.0 ± 0.2 g/m2/yr) especially during snowmelt and storms. A large fraction of DOC (2.5 ± 0.2 g/m2/yr) percolates to the deeper subsurface, fueling deep respiration to generate DIC. DIC is exported predominantly from the deeper subsurface (7.1 ± 0.4 g/m2/yr, compared to 1.3 ± 0.3 g/m2/yr from shallow soils). Deep respiration reduces DOC and increases DIC concentrations at depth, leading to commonly observed DOC flushing (increasing concentrations with discharge) and DIC dilution patterns (decreasing concentrations with discharge). Surprisingly, respiration processes generate more DIC than weathering in this carbonate‐based catchment. These findings underscore the importance of vertical connectivity between the shallow and deep subsurface, highlighting the overlooked role of deep carbon processing and export.


Introduction
Dissolved inorganic carbon (DIC) and organic carbon (DOC) are important measures of stream water quality and ecosystem health.Excessive DOC can cause water browning, mobilize toxic metals (Laudon et al., 2012), and form carcinogenic disinfectant byproducts during water treatment (Leonard et al., 2022).DIC is the sum of dissolved carbonate species (CO 2(aq) , HCO 3 -, and CO 3 2-), the concentrations and fluxes of which regulate the extent of CO 2 evasion to the atmosphere.DOC and DIC regulate lateral carbon fluxes from land to rivers, an increasingly recognized flux in global carbon cycling (Marx et al., 2017;Öquist et al., 2009;Raymond et al., 2013).
DOC and DIC are generated by processes including soil respiration (heterotrophic via microbial decomposition, and autotrophic via roots) and chemical weathering and exported via distinct flow paths to streams (Campeau et al., 2017;Keller, 2019;Marx et al., 2017;Wen et al., 2022).DOC can be further oxidized to DIC via in-stream processing (Grandi & Bertuzzo, 2022).The dynamics of stream carbon observed at a catchment outlet reflect coupled hydrological and biogeochemical processes.DOC typically exhibits flushing behavior, where stream concentrations increase with increasing stream discharge (Boyer et al., 1997;Zarnetske et al., 2018).In contrast, DIC typically exhibits dilution behavior, where stream concentrations decrease with increasing discharge (Bluth & Kump, 1994;Najjar et al., 2020;Stewart, Zhi, et al., 2022;Wallin et al., 2013).These differential export patterns have been attributed to their distinct origins and abundance in the subsurface, and the interconnected processes that produce, consume, store, and transport dissolved carbon (Barnes et al., 2018;Botter et al., 2020;Keller, 2019;Stewart, Shanley, et al., 2022;Zhi et al., 2019).
In particular, the lateral fluxes of dissolved carbon have been increasingly recognized as important in integrating carbon processes from land to streams to oceans (Barnes et al., 2018;Battin et al., 2023;Hodges et al., 2021;Keller, 2019;Regnier et al., 2022).These fluxes are sensitive to changing weather patterns and climate conditions.Dry conditions can promote CO 2 production and vertical carbon export (e.g., soil CO 2 efflux), while wet conditions enhance DOC production and lateral carbon export (Wen et al., 2022).In addition, deep carbon in bedrock fractures beneath soils has been shown to represent a significant store of potentially labile organic carbon (Gross & Harrison, 2019;Rumpel & Kögel-Knabner, 2011;Sánchez-Cañete et al., 2018;Wan et al., 2018).Carbon respiration in the deep subsurface has been observed to occur at rates just as high as observed in soils (Hasenmueller et al., 2017;Tune et al., 2020) and may represent a missing piece of the global carbon budget calculation (Harper & Tibbett, 2013).The dynamics of lateral export and deep carbon processing, including their rates and dependence, however, have remained poorly understood and quantified such that it is challenging to compare them with measurements of soil CO 2 effluxes to the atmosphere.
One of the main challenges in addressing this knowledge gap is the lack of accessibility to the subsurface.Rates of reactions, such as soil respiration, and their dependence on temperature and soil moisture are often inferred from time series of soil CO 2 concentrations, CO 2 fluxes measured at the ground surface, and dissolved carbon concentrations at the stream outlet (Davidson & Trumbore, 1995;Jian et al., 2021).Carbon transformation and transport processes are sensitive to changes in hydrological regimes and may shift under changing weather and climatic conditions, yet they are rarely studied across disciplinary boundaries (Brookfield et al., 2021;Duvert et al., 2018;Grimm et al., 2003), resulting in the general lack of integrated conceptual frameworks on the timing, magnitude, and location of these processes.
Here we ask the question: when, where, and how much is dissolved carbon produced and exported?In particular, how and to what extent do production and export rates of dissolved carbon vary (a) across season (temperature and precipitation conditions) and (b) in the shallow and deep subsurface?We answer these questions using discharge and stream chemistry data from Sleepers River and a catchment-scale reactive transport model.We hypothesize that (a) DOC is primarily produced in shallow soils and depends on higher temperatures when microbial activities peak; DIC is derived primarily from soil respiration in shallow soils that depends more on temperature and from carbonate weathering in the deeper subsurface that is more regulated by hydrology due to its fast rates, long contact time, and quick approach to equilibrium; (b) the export of DOC, given its typically high concentrations in the shallow subsurface and flushing behavior, is primarily driven by rapid shallow soil flow and peaks under wet conditions; in contrast, DIC is primarily derived from the deeper subsurface with more abundant carbonate such that its export is driven by the relatively steady flow of groundwater across the year.
We test these hypotheses using a reactive transport model (RTM) that integrates and differentiates the effects of individual processes at the catchment scale (Li, 2019;Li et al., 2021), therefore illuminating hydrological and biogeochemical processes in the subsurface.As an example, at Shale Hills, a temperate, forested catchment in central Pennsylvania, reactive transport modeling has shown that DOC production in soils hinged more upon temperature than hydrology and peaked in the summer, whereas DOC export was mostly driven by hydrologic regimes (Wen et al., 2020).In addition, sorption on soil surfaces serves as a storage mechanism for DOC in the summer when land-stream connectivity is minimal.Hydrological regimes and characteristics also regulate the magnitude of vertical CO 2 effluxes and lateral carbon export to streams (Wen et al., 2022).These hydrobiogeochemical insights from RTMs show promise of unraveling the interconnected processes that affect the timing, location, and magnitude of dissolved carbon production, export, and stream dynamics.Here we use a recently developed reactive transport model, BioRT-HBV (Sadayappan et al., 2024), which couples the widely used HBV-light hydrology model (Seibert & Vis, 2012) and the BioRT model for solute transport and biogeochemical reactions (Zhi et al., 2022).

Site Description
The study site is the W-9 catchment at Sleepers River Research Watershed (SRRW), located in northeastern Vermont with a drainage area of 0.405 km 2 (Kendall et al., 1999).SRRW has been an active research watershed with intensive data collection since 1959 (Shanley, 2000).The catchment is entirely forested and has a humid, continental climate (Armfield et al., 2019;Shanley et al., 2015).Soils range from poorly drained Histosols (e.g., Cabot, Peacham, Lupton series) in riparian areas and headwater swamps to well-drained Spodosols (e.g., Tunbridge, Lyman, Berkshire series) and Inceptisols on hillslopes (e.g., Vershire, Glover, Dummerston, Buckland series; Figure 1a) (DeKett & Long, 1995;Shanley et al., 2004).Soils (approximately 60-80 cm) developed from a layer of glacial till (approximately 1-4 m) and the underlying bedrock of carbonaceous quartz-mica phyllite with impure rusty-weathering limestone beds.The vadose zone varies with location and hydrologic conditions, potentially interacting with all subsurface materials, though generally resides within 0.5-2 m of the land surface.Peak discharge often occurs during spring snowmelt in late March-early April, though additional discharge peaks occur throughout the year in response to precipitation events.Snowpack builds throughout the winter before it is depleted with increased temperatures that drive spring snowmelt (Chalmers et al., 2019).
The mean annual precipitation and air temperature are 1,320 mm/yr and 4.6°C, respectively (Armfield et al., 2019;Sebestyen et al., 2009).Discharge and precipitation vary seasonally (Figure 1b) with snowmelt being the major hydrological event, contributing nearly half of total annual stream discharge (Shanley, 2000).The magnitude of snowmelt varies from year to year; precipitation falls throughout the year (Figure 1b).Saturation excess overland flow occurs during snowmelt and extreme storm events (Dunne & Black, 1970b, 1971;Kendall et al., 1999).Groundwater supplies baseflow such that the stream at the gage is perennial (Shanley et al., 2015).
The two zones are conceptual and not assigned explicit depths in HBV-light, as these can vary between catchments.The shallow zone represents the shallow subsurface with more weathered soil minerals and soil organic matter, into which the water table rises under high flow conditions such that water transit times are relatively short (Sullivan et al., 2016).The deep zone represents the deeper subsurface with less-weathered minerals and bedrock and harbors slower, older water that supplies baseflow under low flow conditions (Frisbee et al., 2013).Daily hydrological and meteorological data from SRRW, including discharge (determined from stage-discharge relation), air temperature, precipitation, and potential evapotranspiration were used as model input (Shanley et al., 2021).Potential evapotranspiration values were not measured at SRRW, so we used values from the nearby Hubbard Brook Experimental Forest in New Hampshire (approximately 71 km away) with similar climate and vegetation as a proxy (Green et al., 2021).
BioRT uses the hydrology output from HBV-light (i.e., time series of water fluxes, soil moisture, and dynamic storage) to simulate solute transport (e.g., advective fluxes) and biogeochemical reactions.In BioRT, users can choose the reactions to be included.Typically, reactions occur in the SZ and DZ, while the Surface (SF) generates flow with the chemistry of precipitation and/or meltwater, when present.BioRT also has the capability of simulating surface reactions that mobilize solutes and sediments.BioRT outputs reaction rates and concentrations of aqueous and solid concentrations of relevant chemical species in each zone (SF, SZ, DZ) and in the stream at the daily scale.

Governing Equations in BioRT
BioRT solves the governing equations for the SF, SZ, and DZ for primary chemical species in the user-defined reaction network, following a typical reactive transport modeling approach (Lichtner, 1985;Steefel & Mac-Quarrie, 2018).The primary species make up the foundation of the chemical system, such that secondary chemical species can be expressed in terms of primary species' concentrations through equilibrium relationships and laws of mass action.The model solves differential equations for the concentrations of primary species, whereafter the concentrations of secondary species can be calculated.The governing equations for a representative primary species i in each zone are provided below.
In the Surface zone (SF): where the solute concentration in the infiltrating water C infil,i is determined by the mixing of rainwater and snowmelt water: Snow can accumulate at the land surface and is tracked with the equation: In the Shallow Zone (SZ): In the Deep Zone (DZ): where P snow is the precipitation falling as snow, P rain is the precipitation falling as rainfall, Q snowmelt is the flow from snowmelt, Q infil is the infiltrating water flow entering the shallow zone, Q perc is the water percolating (or recharge) from shallow to deep zone.Here V w,SF , V w,SZ , and V w,DZ are the water storage in SF, SZ, and DZ, respectively, and V w,snow is the water storage in the snow.

Reactions
Major carbon reactions are listed in Table 1.The term "OC" is used to broadly represent soil organic carbon that typically consists of a variety of compounds.The heterotrophic respiration represents the process where microbes decompose organic matter into DOC and DIC in the form of CO 2 (aq) (Reaction 1a).The autotrophic respiration is the process where plant roots produce CO 2 (aq) (root respiration) and DOC (root exudation) (Reaction 1b).The relative contributions from heterotrophic and autotrophic respiration are difficult to quantify and can change over time and space.We therefore lumped them as "soil respiration" (Reaction 1) (Barba et al., 2018;Hanson et al., 2000).In other words, we do not model the respiration of DOC to produce CO 2 (aq) individually, rather as  4, K M values reported in "Rate Law" column.b Eq = Equilibrium (thermodynamically controlled, reversible).c TST = Transition State Theory-Equation 5. d The gas-aqueous exchange of CO 2 (Reaction 6) represents the process of the gaseous soil CO 2 dissolving into water, or CO 2 (aq) becoming saturated and forming gaseous CO 2 ."CO 2 (*g)" acts as an immobile reactant (trapped in pore space) that produces CO 2 (aq).Similarly, CO 2 (aq) can enter the gaseous phase when the concentration of CO 2 (aq) reaches solubility limits.e Note that pH (H + ) was included in the model as a primary species for the calculation of carbonate speciation via Reactions 7 and 8. part of a "net" reaction that includes both OC and DOC as sources of CO 2 (aq) and DOC as an intermediate product that is not always fully oxidized to CO 2 (aq).The concentration of DOC therefore reflects the balance between reactions that produce and consume DOC.OC is assumed to be abundant in the SZ, following descriptions of organic-rich soils at W-9 (Shanley, 2000;Shanley et al., 2002).CO 2 (aq) can further speciate rapidly to become other carbonate species such as CO 3 2-and HCO 3 -(Reactions 7 and 8).DIC is calculated as the sum of CO 2 (aq) (or equivalently, H 2 CO 3 (aq)), CO 3 2-, and HCO 3 -.
Note that soil CO 2 produced from Reaction 1 can not only export laterally into streams and vertically into the DZ, as represented by the term (Q SZ + Q perc )C SZ,i in Equation 2, but can also export vertically to the atmosphere via soil CO 2 effluxes.In Equation 2, we do not include vertical soil CO 2 effluxes, as we do not have vertical efflux data.As such, the soil respiration rates calculated here represent the rates that produce dissolved carbon (DOC and DIC) and likely underestimate the overall respiration rates.The CO 2 that enters the stream can also degas into the atmosphere, which we also do not simulate due to the lack of stream CO 2 evasion data.
DOC can sorb onto soil surfaces depending on sorption affinity, pH (Mayes et al., 2012), and available sorption sites ("≡X" in Reaction 2 (Table 1)).Minerals at W-9 such as illite, smectite, chlorite, and oxides, among others, can sorb DOC and retain it in soils (Armfield et al., 2019;Cincotta et al., 2019).This process is thermodynamically controlled by the availability of sorption sites on soil as represented by an equilibrium constant (K eq ) (Langmuir, 1997).The shallow soil zone also includes carbonate weathering (Reaction 3).While calcite parent materials are generally absent from weathered soils, previous studies have found Ca 2+ with organo-mineral complexes and small carbonate crystals, likely neo-precipitates, in soil samples at W-9 (Cincotta et al., 2019).We therefore use "Carbonate-SZ" to represent a range of calcium-and carbonate-bearing minerals that may be present in the SZ and serve as sources of Ca 2+ and DIC to the stream, as it is not feasible to model every carbonate-containing mineral.
In the deeper subsurface, anoxic conditions and mineral-stabilized OC typically lead to much slower decomposition and production of DOC and CO 2 (aq) (Kleber & Johnson, 2010).In addition, DOC in the SZ enters the DZ via recharge and can be further decomposed into CO 2 (Tune et al., 2020(Tune et al., , 2023)).In other words, the deep respiration reaction is also a "net" reaction combining DOC production and consumption rates (Reaction 4).DOC in the DZ must be consumed more than it is produced, otherwise concentrations of DOC would be higher in the DZ than the SZ, which is the opposite of observed depth profiles in subsurface DOC concentrations (Stewart, Shanley, et al., 2022).In fact, if we set the DOC as a reaction product in the DZ, we cannot reproduce the observed DOC concentrations in the stream.We therefore represent deep respiration with both OC and DOC as carbon sources and DIC as a product.The model cannot differentiate the contributions of OC-DZ and DOC to CO 2 (aq) production, as these sources are lumped in Reaction 4, and we do not have data to differentiate the two.
Additional sources of DIC in the deeper zone include the DIC entering with recharge from the SZ and weathering of carbonate bedrock (Shanley, 2000).Carbonate bedrock composition often differs from carbonate neoprecipitates in shallow soils (Macpherson & Sullivan, 2019) such that weathering in the deep zone (Reaction 5) can have different reaction stoichiometry and rates from that in the shallow zone (Reaction 3).Including Ca 2+ adds another constraint for the calibration of carbonate weathering.
Several reactions occur in both the shallow and deep zones.The gas-aqueous exchange of CO 2 (Reaction 6) is regulated by the solubility of CO 2 (aq) prescribed by Henry's Law (Marx et al., 2017;Plummer & Busenberg, 1982).We represent this reaction with CO 2 (*g; refer to footnote in Table 1) as an immobile reactant (trapped in pore space) that can dissolve and become CO 2 (aq), one of the forms of DIC (Reactions 7 and 8).DIC is the sum of CO 2 (aq) (or equivalently, H 2 CO 3 (aq)), CO 3 2-, and HCO 3 -.Equilibrium constants in Table 1 are from literature (Benjamin, 2014;Wen et al., 2020Wen et al., , 2022;;Zamanian et al., 2016).

Reaction Rate Laws
Biotic reactions (including soil and deep respiration) typically follow Monod rate laws, whereas abiotic carbonate weathering reactions follow the rate laws of Transition State Theory (TST) (Lasaga, 1998;Michaelis & Menten, 1913;Monod, 1949).Reaction rates also depend on environmental conditions such as temperature and soil moisture.For biotic reactions, we use a Monod rate law with additional dependencies on temperature and soil moisture: Where r is the reaction rate [mol/s], k is the reaction rate constant [mol/m 2 /s], A is the surface area [m 2 ] and is calculated from the specific surface area [m 2 /g] times mass of reactive materials [g] that represent organic matter content, root abundance, and microbial abundance for microbially mediated reactions (Reactions 1 and 4), C D is the concentration of a limiting electron donor, K M is the half-saturation coefficient [mol/L] of the limiting electron donor, and f(T ) and f(S w ) are the rate dependencies on temperature and soil moisture.In Reaction 1, OC-SZ is abundant such that the corresponding K M is minimal (6.0 × 10 6 mol/L) relative to OC-SZ concentrations in the SZ.In Reaction 4, DOC is a limiting electron donor, as it is primarily sourced through recharged water from the SZ to the DZ, such that the K M is significant (5.0 × 10 3 mol/L) relative to DOC concentrations in the DZ.Given the complexity of processes that can contribute to microbially mediated DOC and DIC production, these reaction rates reflect the overall or net rates of processes that generate DOC and DIC, potentially from multiple reaction pathways.
The f(T ) follows the widely used Q 10 form (Mahecha et al., 2010): , where Q 10 is the increase in reaction rate when temperature increases by 10°C (Davidson & Janssens, 2006).Q 10 values typically range from 1.0 to 3.0 depending on factors such as climate and land cover (Zhou et al., 2009).The f(S w ) follows the form: , where S w,c is the critical soil moisture at which f(S w ) reaches a maximum, and n is the exponent representing the magnitude of reaction rate dependence on soil moisture.Values of n typically range from 0 to 3 depending on soil structure and texture (Hamamoto et al., 2010;Yan et al., 2018).Soil moisture is calculated in each Zone (SZ or DZ) based on the dynamic water storage, porosity, and depth of each zone, as detailed in Supporting Information S1 (Text S4).
Weathering reactions use the TST rate law with temperature and soil moisture dependence (Lasaga, 1998): Where k is the kinetic rate constant [mol/m 2 /s], A is the specific surface area of the solid being weathered [m 2 /m 3 ], IAP is the ion activity product, and K eq is the equilibrium constant.

Hydrology-HBV
We used hydrometeorology, discharge, and stream chemistry data for water years 2016-2017 to calibrate the HBV-light and BioRT models, as these years had sufficient data and represent two consecutive years with distinct stream discharge dynamics.In particular, the 2 years have small and large snowmelt events such that the model can capture the different discharge dynamics.HBV-light was first manually calibrated to develop a better understanding of how different parameters shape discharge dynamics.The built-in Genetic Algorithm and Powell (GAP) optimization tool was then used to obtain a parameter set with an optimized fit of the HBV-light results to discharge data (Okamoto et al., 1998;Seibert, 2000).This tool combines the Genetic Algorithm that is based on the idea of natural selection in biological evolution but computationally expensive, and the Powell method that is fast but can converge on local minima instead of a global minimum.The combined approach takes the advantages of each method while counteracting their respective disadvantages (Okamoto et al., 1998).
We also utilized the built-in Monte Carlo simulation tool to run 1,000,000 simulations with randomly generated parameter sets (Seibert & Vis, 2012).Model performance was measured with commonly used Nash-Sutcliffe Efficiency (NSE), Kling-Gupta Efficiency (KGE), and the Non-Parametric Kling-Gupta Efficiency (NPE) (Gupta et al., 2009;Nash & Sutcliffe, 1970;Pool et al., 2018) et al., 2018).Details of the HBV-light model and its calibration parameters can be found in Supporting Information S1, HBV-light manual (Seibert, 1996), and other literature (Seibert, 1999;Seibert & Vis, 2012).In addition to performance metrics, we also used field observations at Sleepers River to guide model evaluation (Shanley et al., 2015).Specifically, discharge is known to be dominated by groundwater (Q DZ ) even under high flow conditions like snowmelt and storm events (Kendall et al., 1999;Shanley et al., 2002;Stewart, Shanley, et al., 2022).Thus, we would expect Q DZ to contribute approximately 50% or more of annual discharge (Table S1 in Supporting Information S1).Furthermore, overland flow Q SF is known to occur primarily at very high flow such that its contribution would be minimal compared to Q SZ and Q DZ (Dunne & Black, 1970b, 1971).

Biogeochemistry-BioRT
We used stream chemistry data (∼weekly DOC, DIC, and Ca 2+ ) to calibrate BioRT.We also used subsurface water chemistry data from piezometers and soil lysimeters (monthly to quarterly; Figure 1a) to set initial solute concentrations in subsurface zones (Matt et al., 2021).Given consistent stream pH conditions between 6.3 and 10.
, K 1 = 10 6.35 and K 2 = 10 10.33 under 25°C and standard atmospheric pressure (Benjamin, 2014).Concentrations were used to approximate activities because waters are relatively dilute in Sleepers River.DIC concentrations were calculated as The model was first manually calibrated by adjusting reaction parameters including stoichiometric coefficients, reaction rate constants (logk), specific surface area (SSA) of solid reactants (OC for respiration and carbonate for weathering reactions), temperature dependence f(T ) via Q 10 , and soil moisture dependence f(S w ) via parameters S w,c and n.The reaction stoichiometry needed to be adjusted as the stream chemistry reflects the "effective" reactions of OC and multiple minerals (instead of single pure minerals) at the catchment scale.Modeled stream concentrations were compared to stream chemistry data, and the performance of each model was evaluated using NSE, KGE, and NPE, following standards from previous studies, where values closer to 1 indicate better performance (Knoben et al., 2019;Moriasi et al., 2007;Wen et al., 2020).
Following manual calibration, we ran a Monte Carlo analysis for the BioRT model by pairing BioRT parameter sets with four hydrology cases from HBV that represented the range of hydrological behavior from the 110 Monte Carlo cases with the best model performance.For each hydrology case, we ran 5,000 randomly sampled BioRT parameter sets, as detailed in Supporting Information S1 (Text S2 and Table S4).These combinations did not produce a sufficiently good fit (NSE > 0 for all three solutes, DOC, DIC, and Ca 2+ ).We ran a second round of 1,000 Monte Carlo cases with narrower ranges of parameter values (ranges given in Table S5 in Supporting Information S1).We also ran each HBV parameter set with the original manually calibrated BioRT parameters for comparison.These simulations resulted in 64 cases with NSE > 0 for DOC, 14 cases with NSE > 0 for DIC, and 20 cases with NSE > 0 for Ca 2+ (including the cases with manually calibrated BioRT parameters).Among these cases, only four cases had positive NSE values for all three solutes, three of which used the manually calibrated BioRT parameters and had NSE values >0.3 for all solutes.
All three HBV cases with sufficient BioRT performance (NSE > 0.3) had annual Q DZ contributions from 64% to 69%, indicating the hydrology needed to be in a narrow range for stream chemistry to be reproduced.We ran an additional 32 HBV Monte Carlo cases with discharge partitioning in this range (1%-5% Q SF , 25%-34% Q SZ , 64%-69% Q DZ ).These additional simulations led to seven total cases with NSE values >0.45 for all solutes (Table S2 in Supporting Information S1).We used these seven cases to calculate the mean and standard deviations of HBV output (e.g., Q, Q SF , Q SZ , Q DZ ) and BioRT output (e.g., concentrations and reaction rates).

Production and Export Rates at the Catchment Scale
BioRT outputs daily time series of reaction rates ("R ReactionName ") for each kinetically controlled reaction in each zone (Table 1).These reaction rates were used to calculate daily rates of production (R p ) for DOC, DIC, and Ca 2+ at the catchment scale (including SZ and DZ) using reaction stoichiometry and respective rates in Water Resources Research 10.1029/2023WR035940 the two zones (Text S3 in Supporting Information S1).For example, the net DOC production rates are the sum of DOC generation from Resp-SZ in SZ and DOC consumption rates from Resp-DZ in DZ.We also calculated the daily net export rates (R e ) for each solute as the product of concentration [mg/L] and discharge [mm/d] for each zone and for the stream (equivalent to solute fluxes or loads).All production and export rates are in units of mg/m 2 /d.

Discharge and Flow Paths
Discharge peaked during snowmelt and varied across the year, as shown by data and HBV model results (Figure 3).These HBV simulations had NSE values ranging from 0.71 to 0.82, which were lower than the optimized GAP case (NSE = 0.85) but had the best reproduction of snowmelt peaks in both years and discharge partitioning that was consistent with existing literature.Deep groundwater flow Q DZ contributed to discharge persistently across the year, whereas soil water flow Q SZ and surface flow Q SF varied significantly, peaking at snowmelt and dropping to zero in summer and fall.Surface flow was fleeting and lasted only a few days during snowmelt.Overall, Q DZ contributed 67 ± 1.7% to annual discharge, compared to 30 ± 2.7% and 3 ± 1.2% for Q SZ and Q SF , respectively.Across 110 Monte Carlo cases with highest performance metrics, Q DZ , Q SZ , and Q SF contributed 63 ± 12%, 34 ± 13%, and 3.3 ± 6% to annual discharge, respectively (Table S1 in Supporting Information S1).
During manual calibration, key parameters for reproducing discharge dynamics included the threshold temperature (TT), which controls the temperature at which precipitation falls as snow versus rain, and the snowmelt factor (CFMAX), which determines the snowmelt rate.Other important parameters included UZL, the threshold parameter that determines the timing and magnitude of Q SF generated from the shallow zone, and PERC, the maximum percolation rate of water from the shallow zone to the deep zone.These parameters correspond to important hydrological characteristics, including the annual spring snowmelt event (TT and CFMAX), generation of overland flow during snowmelt and extreme storm events (UZL), and year-round groundwater contributions to the stream (PERC) (Dunne & Black, 1970a, 1970b, 1971;Kendall et al., 1999;Shanley et al., 2002Shanley et al., , 2015)).

Stream Chemistry Dynamics
The water output of HBV-light (Figure 3) and carbon reaction network (Table 1) were used to run BioRT.Stream concentrations (modeled and observed) fall between the modeled concentrations in the SZ and DZ, indicating stream chemistry is derived from the shallow and deep waters in SZ and DZ (Figure 4; light and dark colors, respectively).
Stream DOC concentrations were dynamic, typically holding a consistent baseline value (∼0.5-1 mg/L) during low flow and rising to peak values (up to 6.2 mg/L) following precipitation and snowmelt events (Figure 4a).These temporal dynamics echo the flushing concentration-discharge (CQ) pattern for DOC (Figure 4b) with low concentrations under low discharge and high concentrations under high discharge.The model tended to underestimate DOC concentrations at lowest and highest discharge.In general, however, BioRT captured the temporal dynamics and CQ pattern well (NSE = 0.45 to 0.55, KGE = 0.43 to 0.57, NPE = 0.62 to 0.69).
Stream Ca 2+ followed a seasonal trend (Figure 4c) with highest concentrations (∼27 mg/L) under low flows in late summer and lowest concentrations (∼10 mg/L) under high flows in spring snowmelt.Stream DIC exhibited a similar temporal trend with peaks (∼17 mg/L) during late summer months and troughs (∼6 mg/L) during snowmelt (Figure 4e).These seasonal trends were consistent with the dilution CQ patterns for Ca 2+ and DIC (Figures 4d and 4f).BioRT simulated stream Ca 2+ (NSE = 0.51 to 0.60, KGE = 0.77 to 0.79, NPE = 0.74 to 0.76) and DIC dynamics well (NSE = 0.53 to 0.61, KGE = 0.79 to 0.80, NPE = 0.78 to 0.84); the CQ plots indicate that the model generally captured the range of stream concentrations.The model produced similar dilution CQ patterns but underestimated Ca 2+ and DIC concentrations at the lowest discharge, possibly due to the lack of representation of in-stream processes, as we discuss later.Such underestimation led to a slight threshold-type pattern that is not obvious in the observations.Key factors that reproduced stream DOC included the respiration rates (Resp-SZ and Resp-DZ) that regulate concentrations in the SZ and DZ (Figure S3 in Supporting Information S1).For Ca 2+ , the rate of carbonate weathering in the deep zone (Carbonate-DZ) was essential; in particular, increasing the temperature and soil moisture dependencies helped capture the seasonal fluctuations in DZ.DIC was produced by several reactions in the model (Table 1).The rates of respiration reactions (Resp-SZ and Resp-DZ coupled with CO 2 (*g) gas-aqueous exchange) and carbonate weathering (Carbonate-SZ and Carbonate-DZ) were important in controlling DIC concentrations in subsurface zones (Figure S3 in Supporting Information S1).The peak DIC concentrations in the stream data, typically occurring in warm and dry summer months, were not reproduced.

Influential Reactions and Drivers
The reactions in the SZ and DZ drove the subsurface source water and streamwater chemistry (Table 1).Note that soil respiration and CO 2 gas-aqueous exchange are coupled processes, so these reactions were always simulated together.When the simulation only included soil respiration and sorption in the SZ (Resp-SZ + Sorption-SZ), stream DOC concentrations were high and primarily responded to changes in hydrology (Figure 5a) with DOC exhibiting a dilution pattern rather than the characteristic flushing pattern observed in data (Figure 5b).Adding deep respiration (Resp-DZ) reduced the DOC concentrations in DZ and therefore in stream, reproducing the DOC flushing pattern.Lastly, the sorption of DOC onto soils (Sorption-SZ) served as a buffering mechanism, modulating the extent of DOC fluctuation most evidently by decreasing the magnitude of DOC peaks in SZ and stream (Figure 5a).
Stream DIC concentrations were relatively low when only Resp-SZ and Resp-DZ reactions were simulated, suggesting that geogenic reactions (Carbonate-SZ and Carbonate-DZ) are important.In fact, the seasonal trend in stream DIC concentrations was not captured by BioRT until Carbonate-DZ was included (Figure 5c).
Soil respiration (Resp-SZ) rates had a strong dependence on air temperature, peaking in summer under high temperatures (Figure 6a) but with no clear dependence on water storage (Figure 6b).In contrast, deep respiration (Resp-DZ) depended more on water storage than air temperature (Figures 6g and 6h), possibly reflecting the fact that deeper respiration rates hinge upon the amount of DOC transported with recharge from SZ to DZ. Carbonate weathering rates (Carbonate-SZ and -DZ) showed different relationships with air temperature and water storage.
Rates of Carbonate-SZ, representing reactions of soil inorganic carbon in the SZ, peaked in spring at high flow conditions and were slightly sensitive to water storage (Figures 6c and 6d).Carbonate-SZ rates were negative in dry and hot summer months, indicating precipitation of soil inorganic carbon.Carbonate-DZ showed a clearer dependence on air temperature than water storage with highest rates occurring with high temperatures in summer Respiration rates increased with air temperature and peaked in summer in both shallow zone (SZ, more positive) and deep zone (DZ, more negative).Carbonate weathering rates increased with air temperature in DZ.Net DOC production rates increased with air temperature, peaking in spring and summer.Net DIC production rates increased with air temperature and water storage with higher rates in spring and summer.

Water Resources Research
10.1029/2023WR035940 across a range of water storage values (Figures 6i and 6j).Overall, both DOC and DIC production rates increased with air temperature and water storage to some extent, with rates peaking in spring and summer months (Figures 6e and 6f).

Production and Export Rates for DOC and DIC
The net DOC production rates (R p ) are the sum of DOC production rates in SZ and consumption rates in DZ.The temporal dynamics of R p generally mirrored those of air temperature (Figure 7) and peaked in summer.In contrast, DOC export rates consistently peaked at spring snowmelt.DIC production rates (Figure 7c) were higher under higher temperatures and wetter conditions of spring and summer, where peaks coincided with discharge throughout the year.DIC export rates followed the temporal pattern of discharge and were less influenced by air temperature.
The production of DOC primarily occurred in the SZ through OC decomposition (Resp-SZ), while deep respiration (Resp-DZ) consumed DOC (Figure 8a).The mean DOC production rate (sum of rates in SZ and DZ) was 3.3 ± 1.3 mg/m 2 /d, where positive rates in the SZ reflect DOC production of DOC and negative rates in the DZ reflect DOC consumption.For DIC, production rates were higher and flashier in the SZ, responding rapidly to high discharge events when high shallow flow drove carbonate reactions to disequilibrium (Figure 8b).DIC production rates in the DZ were higher under warmer and wetter conditions (mean spring R p 17.6 ± 1.0 mg/m 2 / d and mean summer R p 9.1 ± 1.0 mg/m 2 /d).
To tease apart the contributions of DIC from different reactions, we also show respiration rates and carbonate rates in SZ and DZ (Figure 8c).Soil respiration rates (Resp-SZ) generally started to increase in early April and remained high through September, whereas deep respiration rates mostly peaked in wet spring, consistent with high dependence on water storage (Figure 6h).Carbonate reactions differ even more between SZ and DZ.Carbonate-SZ rates were hydrologically driven, peaking sporadically throughout the year during high discharge, indicating dissolution was mostly driven by water flow.Carbonate-SZ rates were notably negative in warm summer, indicating carbonate precipitation rather than dissolution.Rates of DIC production in DZ were lower than the rates in SZ but not negligible.Grouping these reactions into biogenic respiration and geogenic carbonate weathering (Figure 8d) indicates that overall, the biogenic reactions were more driven by air temperature variations, whereas geogenic reactions were more driven by hydrology (discharge events).The biogenic production rates (Resp-SZ and -DZ) were higher than geogenic production rates of DIC (Carbonate-SZ and -DZ).Reaction rates were higher in the SZ such that overall biogenic and geogenic production rates mirrored the trends of SZ reactions more than DZ reactions (Figures 8c and 8d).
DOC export rates from the SZ followed temporal patterns of discharge and shallow soil flow (Q SZ ) (Figure 8e); DOC export rates from the DZ were lower than from the SZ, largely due to lower DOC concentrations in the DZ.Mean DIC export rates from the DZ were high and relatively consistent throughout the year due to high DZ concentrations and high groundwater (Q DZ ) contributions.DIC export rates from the SZ were lower and flashy, corresponding to peaks in production rates, stream discharge, and shallow soil water (Q SZ ) (Figure 8f).Overall, DOC production and export were dominated by the SZ, highlighting the importance of shallow flow paths for DOC.In contrast, although DIC was similarly produced in the SZ and DZ, it was exported primarily from deep flow paths.

Discussion
This work used catchment-scale reactive transport modeling and stream data to illuminate the poorly understood processes of dissolved carbon production in and lateral export from the "invisible" subsurface as well as their timing and magnitude across seasons and hydrological regimes.Such understanding can help determine and quantify how dissolved carbon dynamics will change in the future climate.Modeled discharge (Q) is shown on the right y-axes of panels A and B for reference.DOC was primarily produced in and exported from the SZ.DIC was produced more in the SZ but exported more from the DZ, largely due to the higher concentrations of DIC in DZ and higher contributions of deeper flow throughout the year.Net DIC production was dominated by SZ reactions but included contributions from biogenic (Resp-SZ and Resp-DZ) and geogenic (Carbonate-SZ and Carbonate-DZ) reactions.Biogenic reactions showed a stronger seasonal trend with higher rates in warm summer, while geogenic reactions peaked throughout the year.DIC export rates from the SZ were generally flashy and more sporadic.Results show that DOC was produced and exported mostly from the SZ.DOC was produced by soil respiration in the SZ, often attributed to abundant OC and O 2 in shallow soils (Barnes et al., 2018;Jobbágy & Jackson, 2000;Wen et al., 2020).DOC production rates peaked with higher air temperatures, with more than 95% of annual production occurring during spring and summer (Figure 9), consistent with our hypothesis and existing literature (Gillooly et al., 2001;Wen et al., 2020Wen et al., , 2022)).As DOC was transported to the DZ with recharge, more than half of DOC produced in the SZ (2.5 out of 3.7 g/m 2 /yr) was further oxidized by deep respiration (Resp-DZ), thereby reducing its concentrations in the DZ.In fact, without deep respiration, the model could not reproduce the widely observed flushing pattern for DOC (Figure 5) (Botter et al., 2020;Musolff et al., 2017;Zarnetske et al., 2018).More than 60% of DOC was exported from the SZ compared to the DZ (1.0 vs. 0.6 g/m 2 /yr), mostly driven by high discharge events such as snowmelt (Figures 7 and 8).Approximately 40% of DOC was exported during annual spring snowmelt, underscoring the pivotal role of hydrologic events in laterally exporting DOC to streams (Wen et al., 2020;Zarnetske et al., 2018;Ågren et al., 2007).

DIC
In contrast, DIC was primarily produced in the SZ but exported from the DZ.DIC was derived from both biogenic respiration and geogenic carbonate weathering in the SZ and DZ.DIC production rates in the SZ were about 1.5 times those in the DZ (4.8 compared to 3.2 g/m 2 /yr) (Figures 8 and 9a).The rates were higher under warmer and wetter conditions with >70% DIC produced in spring and summer (Figures 6 and 9c).More than 80% of DIC was exported from the DZ (Figure 9a) with export rates closely mirroring Q DZ contributions to discharge (Figures 3 All rates calculated as mean of seven simulations.Both DOC and DIC were produced more in the SZ than the DZ, but DOC was primarily exported from the SZ, whereas DIC was predominantly exported from the DZ.DIC was primarily produced via biogenic respiration; DIC production rates peaked in the spring and summer.Production rates for DOC peaked in the summer.Export rates of both DOC and DIC peaked in the spring.and 7).The higher DZ export rates arose from higher DIC concentrations in the DZ and the larger proportion of deep flow (Q DZ ) contributing to the stream (Figures 3 and 4).
Production rates for DIC in the SZ were sensitive to both air temperature and water storage (Figure 6).Surprisingly, carbonate weathering rates were higher in the SZ than the DZ (Figure 8), potentially reflecting the faster shallow soil water flow that, when present, often drives carbonate dissolution to occur (Wen et al., 2024).
Alternatively, this may also reflect the presence of other calcium-and carbonate-bearing materials in the SZ beyond calcite that necessitated higher weathering rates to produce an adequate amount of Ca 2+ and DIC to reproduce observed concentrations in shallow soils and in the stream under high flow conditions.Production rates for DIC in the DZ were most sensitive to air temperature (Figure 6i).This is somewhat counterintuitive, as air temperature may have a minimal direct impact on processes in the deep subsurface.However, the deep subsurface responds to inputs from the shallow soils, as water recharges and carries Ca 2+ and dissolved carbon from the shallow subsurface, which is more directly impacted by air temperature (Fontaine et al., 2007;Tune et al., 2020).This highlights the importance of vertical connectivity between the shallow and deep subsurface, even though their dominant processes and rates may differ dramatically (Brantley & Lebedeva, 2011;Xiao et al., 2021).The different seasonal trends of biogenic and geogenic production rates (Figure 8d) further corroborate that temperature sensitivity in the DZ is more a reflection of SZ processes instead of direct air temperature dependence.

Comparison to Other Places
Mean DOC export rates (1.6 ± 0.2 g/m 2 /yr, equivalent to fluxes or loads) from SRRW are similar to those modeled for Shale Hills in Pennsylvania (∼2.3 g/m 2 /yr (Wen et al., 2020),), lower than those calculated and modeled for Bruntland Burn in the Scottish Highlands (∼5.5 g/m 2 /yr (Dick et al., 2014)), and more than twice the mean fluxes calculated for >100 Rocky Mountain streams in the United States (0.61 g/m 2 /yr (Kerins & Li, 2023)), likely reflecting differences in both DOC sources and the amount of water flowing through shallow flow paths across catchments.For example, Shale Hills has relatively similar soils as Sleepers River, while the Bruntland Burn site is dominated by organic-rich soils and blanket peat, likely leading to the much higher DOC export rates reported there (Dick et al., 2014;Herndon et al., 2015).In addition, these sites have similar humid climates but higher contributions of shallow flow to the stream compared to SRRW, likely facilitating higher DOC export rates (Dick et al., 2014;Wen et al., 2020).Mean DIC export rates from SRRW (8.3 ± 0.7 g/m 2 /yr) are more than twice the DIC fluxes reported for the Västrabäcken catchment in Sweden (3.2 g/m 2 /yr (Öquist et al., 2009)) and mean of >100 Rocky Mountain streams (3.3 g/m 2 /yr (Kerins & Li, 2023)), potentially due to the carbonate geology and higher streamflow in a more humid climate at SRRW.Thus, headwater catchments like W-9 at SRRW with organic-rich soils and carbonate lithology may be important cornerstones of both DOC and DIC fluxes from inland waters.

Deep Zone Processes Shape Depth Profiles and CQ Patterns
Results here underscore the importance of DZ biogeochemical processes in regulating the characteristic depth profiles of solutes and the CQ patterns in streams.The relationship between the two has been explained by the Shallow and Deep hypothesis, which states that solute export behavior (i.e., CQ patterns) is shaped by the shifting dominance of source waters under varying hydrological conditions (Stewart, Shanley, et al., 2022;Zhi et al., 2019;Zhi & Li, 2020).More specifically, high discharge conditions are typically dominated by shallow source waters and mostly reflect shallow water chemistry, while low flows are typically dominated by deeper source waters and mostly reflect deeper subsurface chemistry.This hypothesis has been supported by two common observations.One is the observation that biogenic solutes such as DOC and geogenic solutes such as Ca 2+ commonly exhibit flushing and dilution patterns, respectively.This has been shown from individual catchment sites to regional and global scales, despite tremendous variations in landscape structure across sites (Botter et al., 2020;Godsey et al., 2019;Musolff et al., 2017;Stewart, Zhi, et al., 2022;Zhi et al., 2019;Zhi & Li, 2020).The second observation is that biogenic solutes such as DOC typically have higher concentrations in shallow zones and geogenic solutes such as Ca 2+ (or DIC being both biogenic and geogenic) have higher concentrations in deeper zones (Stewart, Shanley, et al., 2022;Zarnetske et al., 2018).
The Shallow and Deep hypothesis has been linked to more abundant organic matter and mineral weathering in shallow and deep zones, respectively.Sensitivity analyses here show that deep respiration plays an essential role in shaping such commonly observed depth profiles and CQ patterns.Cases with only soil respiration (only Resp-SZ, no deep respiration Resp-DZ) result in similar DOC concentrations in the SZ and DZ, and a chemostatic CQ pattern with almost constant concentrations (Figure 5) despite the abundant organic matter in the SZ and weatherable minerals in the DZ.In other words, if the DZ does not further oxidize recharged DOC from the SZ, DOC in the DZ eventually reaches similar concentrations as in the SZ.The characteristic depth profile of DOC decreasing with depth and the corresponding flushing CQ pattern can only be produced when deep respiration (Resp-DZ) consumes DOC from the SZ.Similarly, deep respiration produces DIC in the DZ in addition to carbonate weathering (Carbonate-DZ), which elevates DIC concentrations at depth and leads to the dilution CQ pattern.As a result, despite faster DIC production in the SZ, concentrations of DIC are higher in the DZ.
The Shallow and Deep hypothesis suggests that vertical subsurface structure and depth profiles are the first-order control on stream CQ patterns, rather than lateral and landscape structure.Growing evidence has lent support to the importance of vertical structure in shaping stream CQ patterns.This contrasts the concepts of biogeochemical hot spots and hot moments or ecosystem control points that underscore the disproportionately larger impacts of relatively small patches of a landscape on hydrological and biogeochemical processes in some environments (Bernhardt et al., 2017;Krause et al., 2017).It is possible that more complex CQ patterns, including hysteresis and threshold behaviors, originate from comparative controls of vertical structure and lateral structure (Herndon et al., 2015;Knapp et al., 2022;Sullivan et al., 2019).For example, the segmented CQ pattern of DOC at the Shale Hills CZO in Pennsylvania could emerge from shifting hydrologic connectivity in the catchment, where lower flows reflect the organic-rich riparian zone while higher flows reflect additional inputs of water from organicpoor hillslopes (Herndon et al., 2015;Wen et al., 2020).

More Biogenic Than Geogenic DIC in Carbonate-Based Sleepers River
DIC was derived from carbonate weathering and respiration reactions in both shallow and deep zones (Figure 8).Previous works have noted the general dominance of DIC in the deeper subsurface compared to shallow soils, but the relative contributions of geogenic versus biogenic DIC sources vary and are challenging to quantify (Campeau et al., 2017;Doctor et al., 2008;Duvert et al., 2020).Here we show that the production of biogenic DIC through respiration occurred more rapidly (5.8 ± 0.4 g/m 2 /yr) than geogenic production via carbonate weathering (2.2 ± 0.2 g/m 2 /yr, Figure 9).Despite slower rates of geogenic relative to biogenic production, the inclusion of carbonate weathering was essential for reproducing DIC concentrations in the subsurface and stream (Figure S3 in Supporting Information S1).Results here show that ∼70% of DIC was biogenic, consistent with a previous study at W-9 suggesting that the average contribution of soil CO 2 to groundwater DIC was higher than 60% (Doctor et al., 2008).A recent modeling study in a catchment with karst geology also found that ∼60% of riverine DIC was produced from biogenic sources (Wen et al., 2024).Similarly, a stable carbon isotope analysis that differentiated sources of biogenic and geogenic DIC in a catchment with dolomite geology showed that 83%-94% of the total DIC flux was biogenic (Duvert et al., 2020).This indicates that even in carbonate-based catchments and streams, biogenic processes can provide more DIC than carbonate weathering processes.

Implications for Future With Climate Change and Human Disturbance
The timing and magnitude of dissolved carbon concentrations in streams depend on the temperature and hydrologic conditions that drive carbon transformation and transport, which may be altered, reduced, or amplified in a changing climate.Precipitation, temperature, evapotranspiration, and discharge at W-9 are projected to increase (Pourmokhtarian et al., 2017).In the northeast United States, temperature is expected to increase more in summer compared to winter (Hayhoe et al., 2006), which may escalate DOC and DIC production in summer and further elevate DOC and DIC concentrations in soil water overall and in the stream during large summer storms (Sebestyen et al., 2009;Shanley et al., 2002).
Precipitation has also been projected to induce more frequent and intense storm events (Melillo et al., 2014).These changes may generate more flushing events that elevate stream DOC concentrations (Sebestyen et al., 2009).DOC export is primarily driven by hydrological processes, as exhibited by the co-occurrence of DOC and discharge peaks, the flushing CQ pattern (Figure 4), and the net export rates of DOC (Figure 8) (Ruckhaus et al., 2023).DIC production and export may similarly increase, although the extent of increase may be smaller, as DIC production and export rates remained relatively consistent between the 2 years simulated here, with only slightly higher DIC export in the wetter water year 2017 (Figure 7).Increasing water temperatures may decrease solubility of CO 2 in water and increase CO 2 evasion (Schelker et al., 2016).In climates becoming drier, increased temperatures and evapotranspiration may ultimately escalate solute concentrations as discharge dwindles, though this may also reduce the total lateral export fluxes (Kerins & Li, 2023;Li et al., 2022;Stewart, Zhi, et al., 2022).
Smaller snowpack and earlier snowmelt may reduce and shift snowmelt peaks, which may reduce or amplify DOC loads depending on how flow paths are altered (Sebestyen et al., 2009).For example, sustained shallow soil flow may extend the period of high DOC export; whereas if the snowmelt peak is smaller and not significantly broader, DOC export may decrease, as in water year 2016 (Figures 1 and 7).Rain-on-snow events, however, may add episodic pulses of stream discharge and DOC export, though the magnitude may depend on antecedent snowpack.In fact, smaller snowpack has been shown to elevate discharge peaks during rain-on-snow events (Juras et al., 2021).Smaller snowpacks may also reduce insulation of the soil and induce colder conditions and slower reaction rates, therefore limiting DOC production (Seybold et al., 2022).
Other long-term disturbances, such as changing soil properties, may also affect DOC production and export (Sullivan et al., 2022).The sorption of DOC onto soils (Sorption-SZ), for example, can serve as a storage and release mechanism for DOC (Wen et al., 2020) and minimize concentration fluctuations under varying flow conditions, while low sorption capacity amplifies concentration fluctuations (Figure 5).Sorption therefore may be key to modulating concentrations and export loads during large discharge events (Wen et al., 2020).Sorption capacity may be sensitive to changes in temperature and soil solution, such that increasing temperatures and precipitation may increase desorption and liberation of DOC from soils (Georgiou et al., 2022).Future work exploring the influence of sorption could help elucidate the role of soil properties on long-term temporal trends in stream DOC (Adler et al., 2021;Cincotta et al., 2019).

Model Limitations and Future Research Directions
This work shows that catchment-scale reactive transport models such as BioRT-HBV are useful for unraveling coupled hydrological and biogeochemical processes in the "invisible" subsurface.The model structure of BioRT, inherited from HBV-light, simplifies catchments into two major subsurface zones and three major flow paths.Conceptually, the designation of subsurface zones in BioRT emphasizes the importance of subsurface vertical structure (Xiao et al., 2021), not the lateral landscape structure such as riparian versus upland areas in a catchment.This entails that the model represents the "average" lateral structure without differentiating distinct units across landscape positions.
Although HBV-light has been used and tested extensively across catchments globally (Bhattarai et al., 2018;Nonki et al., 2021;Poméon et al., 2017), one common challenge in hydrology modeling is equifinality, where multiple parameter sets can lead to similar model performances.The Monte Carlo analysis for HBV-light revealed that many parameter sets performed equally well in reproducing discharge, challenging the identification of cases that represent "real" conditions (Underwood et al., 2023).Similarly, even with a limited number of reactions, the complexity of representing reaction kinetics and thermodynamics through rate laws leads to many parameters in BioRT, as is common in RTMs (Maher et al., 2006;Saaltink et al., 2003;Stolze et al., 2023;Wen et al., 2020;Xu et al., 2022;Zhi et al., 2019).One might expect this to result in equifinality issues, where many BioRT parameter sets perform well when compared to data.Our experience here is the opposite.The BioRT Monte Carlo analysis (over 20,000 BioRT simulations) here tested a wide range of HBV cases with different BioRT parameter combinations, and none performed as well as the manually calibrated BioRT case.This was unexpected, as 1,000,000 HBV cases led to over 100 cases that fit discharge data.This difference in calibration for hydrology and biogeochemistry models is beyond our understanding, but we speculate that it may be due to the high interdependence between solute concentrations and various reaction rate parameters, leading to a very small parameter space that can reproduce observed concentrations.Thus, the equifinality issue observed in calibrating HBV-light (and commonly found with other hydrological models) may be partially counterbalanced by the additional constraints from observed stream and subsurface water chemistry.The best performing HBV simulations that reproduced stream chemistry had a narrow range of flow partitioning, with deeper groundwater flow at ∼64%-69% of annual discharge.Although still limited by computational constraints and model structure, the BioRT calibration process indicates that stream chemistry can provide additional constraints for streamflow generation (Beven, 2020).
It is important to note that the calculation of soil respiration rates that produce DOC and DIC does not include the vertical fluxes of CO 2 back to the atmosphere, which will inevitably underestimate the overall soil respiration rates.The soil respiration rates here are therefore the rates of dissolved carbon production, not vertical soil CO 2 Water Resources Research 10.1029/2023WR035940 STEWART ET AL.
effluxes.Such vertical fluxes are usually measured at the ground surface and have been observed to be higher in warm summer conditions when microbial activity is high but low water flow minimizes lateral DIC export (Wen et al., 2022).As such, the magnitude of underestimation in this study may be highest in summer compared to other seasons.Annual soil CO 2 efflux from Shale Hills, a site with relatively similar soils as W-9, has been estimated at 870 g C/m 2 /yr (Hodges et al., 2019;Shi et al., 2018).Other published rates of soil respiration, as measured by soil CO 2 efflux, are typically much higher (530-874 g C/m 2 /yr) than the biogenic production rates of DIC calculated in this study (Bae et al., 2015;Davidson et al., 1998;Fahey et al., 2005).However, the calculated DIC export rates are consistent with observed lateral fluxes from other studies, as reported in Section 4.1.
BioRT currently includes biogeochemical reactions in the subsurface without in-stream processes.In-stream production of CO 2 (stream respiration) often escalates DIC concentrations under warm, dry conditions, which may contribute to the consistent model underestimation of stream DIC from July to September (Figure 4) when warm water and low flow potentially amplify aquatic ecosystem respiration and in-stream CO 2 production.
Existing literature has used isotopic analyses to quantify in-stream carbon processes, but this data was not available for W-9 (Campeau et al., 2018;Dawson et al., 2001;Duvert et al., 2020;Perdrial et al., 2013;Solano et al., 2023).Similarly, as stream CO 2 outgassing rates at W-9 have not been measured, we do not simulate outgassing of CO 2 from streams here.This can lead to underestimation of DIC production rates (Butman & Raymond, 2011;Horgby et al., 2019;Wallin et al., 2013).In addition to model constraints, these limitations highlight the need for co-located measurements of different forms of carbon, including gaseous, aqueous, and particulate carbon, in order to thoroughly understand and quantify carbon processing and export in both vertical and lateral directions.

Conclusion
Here we ask: when, where, and to what extent is dissolved carbon produced and laterally exported at the catchment scale?Overall, model results highlighted the importance of deep zone processes and the differential controls on DOC and DIC.DOC is primarily produced by shallow zone processes whereas DIC is driven by processes in both shallow and deep zones.The shallow and deep zones produce comparable amounts of DIC, but the deep zone exports the majority of DIC.
Results here showcase the ability of BioRT to answer these questions by revealing the major dynamics of coupled hydrology, transport, and carbon reactions while keeping a relatively small number of processes and parameters in a parsimonious model structure.In particular, DOC was produced mostly in warm summer (2.7 ± 1.1 g/m 2 /yr) and in the shallow subsurface (mean production rate R p , 3.7 ± 0.6 g/m 2 /yr) but was consumed in the deeper subsurface ( 2.5 ± 0.2 g/m 2 /yr).This led to lower DOC concentrations in deeper zone, the commonly observed flushing pattern (increasing concentrations with discharge), and much lower export from the deeper zone (R e , 0.1 ± 0.1 g/m 2 /yr) compared to the shallow zone (1.0 ± 0.2 g/m 2 /yr).In contrast, DIC was produced in both shallow (R p , 4.8 ± 0.7 g/m 2 /yr) and deep subsurface (3.2 ± 0.4 g/m 2 /yr), with higher rates in wet spring and warm summer.DIC was exported more from the deeper subsurface (7.1 ± 0.4 g/m 2 /yr) than the shallow subsurface (1.3 ± 0.3 g/m 2 /yr), attributed to higher DIC concentrations at depth and more deep water flow.Surprisingly, DIC originated more from soil respiration than carbonate weathering, highlighting the importance of biogenic DIC production.In the northeast United States, climate change is expected to manifest as increased frequency and intensity of storm events, warming temperatures, and altered timing and magnitude of snowmelt (Hayhoe et al., 2006;Melillo et al., 2014;Pourmokhtarian et al., 2017).These changes will be influential in shaping the future of dissolved carbon production and export, as well as stream water chemistry and aquatic ecosystem health.

Figure 1 .
Figure 1.(a) Map of W-9 catchment at Sleepers River Research Watershed in Vermont, United States; (b) Time series of specific discharge (Q, mm/d, discharge normalized by drainage area), snow water equivalent (SWE, cm), precipitation (mm/ d), and air temperature (°C) for water years 2016-2017 indicate W-9 is seasonally snow-dominated.Peak discharge often occurs during spring snowmelt in late March-early April, though additional discharge peaks occur throughout the year in response to precipitation events.Snowpack builds throughout the winter before it is depleted with increased temperatures that drive spring snowmelt(Chalmers et al., 2019).

Figure 2 .
Figure 2. (a) Conceptual diagram of reactive transport model BioRT-HBV structure; (b) Reaction network implemented in BioRT, including soil respiration, carbonate weathering, deep respiration, DOC sorption (DOC ⇌ ≡XDOC), gas-aqueous CO 2 exchange, and DIC speciation.OC is soil organic carbon; ≡XDOC is DOC sorbed onto soil surfaces.Note that the shallow zone (SZ) and deep zone (DZ) correspond to upper zone and lower zone in HBV-light, such that Q SF , Q SZ , and Q DZ correspond to the flows defined in HBV-light as Q 0 (fast flow), Q 1 (intermediate flow), and Q 2 (slow flow), respectively.The BioRT model incorporates input precipitation chemistry and a user-defined reaction network into the SZ and DZ, such that distinct chemical signatures are calculated in each zone for each flow path, which eventually emerge in the stream.

Figure 3 .
Figure 3. (a) Discharge data (gray circles) and model output (lines with shaded region, mean ± one standard deviation for seven simulations) for discharge (black) and actual evapotranspiration (ET, red) for water years 2016-2017 at W-9 catchment, Sleepers River Research Watershed.(b) Discharge partitioning from model output including mean ± one standard deviation for surface flow (Q SF ), shallow soil water flow (Q SZ ), and deeper groundwater flow (Q DZ ).Deeper groundwater Q DZ persistently provided water flow, whereas soil water flow and surface flow tended to vary significantly over time, peaking at snowmelt and dropping to zero in summer and fall.Surface flow was generally fleeting, lasting only a few days during snowmelt.

Figure 4 .
Figure 4. Left: Time series of observed stream concentrations (symbols) versus modeled concentrations (lines with shaded regions, mean ± one standard deviation for seven simulations) in stream, shallow zone (SZ), and deep zone (DZ) for (a) dissolved organic carbon (DOC), (c) calcium (Ca 2+ ), and (e) dissolved inorganic carbon (DIC) at W-9 catchment, Sleepers River Research Watershed.Right: Corresponding concentration-discharge (CQ) patterns (observed data vs. modeled mean) for (b) DOC, (d) Ca 2+ , and (f) DIC.Light and dark gray lines correspond to modeled concentrations in SZ and DZ, respectively; intermediate gray lines correspond to modeled stream concentrations.Note that the upper y-axis limit for DOC concentration in panel A is 15 mg/L to show the stream and DZ concentrations, but peak modeled concentrations in the SZ reach ∼30 mg/L.Data and model results indicated high variation across the year, with DOC typically exhibiting high concentrations under high discharge conditions, and Ca 2+ and DIC showing the opposite pattern.

Figure 5 .
Figure 5. Left: Time series of reactive transport model BioRT output for (a) dissolved organic carbon (DOC) and (c) dissolved inorganic carbon (DIC) with different reactions included (Resp-SZ, Resp-DZ, Sorption-SZ, Carbonate-SZ, and Carbonate-DZ).Right: Corresponding CQ patterns for (b) DOC and (d) DIC.Note that soil respiration and CO 2 gas-aqueous exchange are coupled processes, so these reactions were always simulated together.The sorption reaction (Sorption-SZ) attenuated concentration fluctuations of DOC, especially for peak concentrations.Resp-DZ consumed DOC from the SZ and was essential in reproducing the flushing pattern of DOC (b).DIC was produced by respiration and weathering reactions in the SZ and DZ, all of which were needed to reproduce temporal stream DIC dynamics.SZ = shallow zone, DZ = deep zone, Resp = respiration, Sorption = DOC sorption onto soils, Carbonate = carbonate weathering.

Figure 6 .
Figure 6.Model output showing rate dependence on air temperature and water storage by season for (a, b) Soil respiration (Resp-SZ), (c, d) Shallow carbonate weathering (Carbonate-SZ), (g, h) Deep respiration (Resp-DZ), and (i, j) Deep carbonate weathering (Carbonate-DZ); (e) Net dissolved organic carbon (DOC)production rate dependence on air temperature; (f) Net dissolved inorganic carbon (DIC) production rate dependence on air temperature.Reaction rates are in mmol of how much reactant solid changed (organic carbon (OC) for respiration, carbonate for weathering)/m 2 /d.Net production rates are in mg C (in DOC or DIC)/m 2 /d.Model results are from one example simulation (MC27).Respiration rates increased with air temperature and peaked in summer in both shallow zone (SZ, more positive) and deep zone (DZ, more negative).Carbonate weathering rates increased with air temperature in DZ.Net DOC production rates increased with air temperature, peaking in spring and summer.Net DIC production rates increased with air temperature and water storage with higher rates in spring and summer.

Figure 7 .
Figure 7. (a) Modeled stream discharge and observed air temperature; (b) Net rates of dissolved organic carbon (DOC) production (R p ) and export (R e ); (c) Net rates of dissolved inorganic carbon (DIC) production (R p ) and export (R e ).Model results from one example simulation (MC27).All rates are in mg C/m 2 /d.Both DOC and DIC production rates peaked in warm spring and summer, while their export rates peaked in wet spring, especially during snowmelt.

Figure 8 .
Figure 8. Net production rates (R p ) of (a) dissolved organic carbon (DOC) and (b) dissolved inorganic carbon (DIC) in the shallow zone (SZ) and deep zone (DZ); (c) Net production rates of DIC by different reactions; (d) Net production rates of DIC by biogenic (Resp-SZ + Resp-DZ) and geogenic (Carbonate-SZ and Carbonate-DZ) reactions; Net export rates (R e ) of (e) DOC and (f) DIC in SZ and DZ.Model results are from one example simulation (MC27).All rates are in mg C/m 2 /d.Modeled discharge (Q) is shown on the right y-axes of panels A and B for reference.DOC was primarily produced in and exported from the SZ.DIC was produced more in the SZ but exported more from the DZ, largely due to the higher concentrations of DIC in DZ and higher contributions of deeper flow throughout the year.Net DIC production was dominated by SZ reactions but included contributions from biogenic (Resp-SZ and Resp-DZ) and geogenic (Carbonate-SZ and Carbonate-DZ) reactions.Biogenic reactions showed a stronger seasonal trend with higher rates in warm summer, while geogenic reactions peaked throughout the year.DIC export rates from the SZ were generally flashy and more sporadic.Resp = respiration, Sorption = DOC sorption onto soils, Carbonate = carbonate weathering.
Figure 8. Net production rates (R p ) of (a) dissolved organic carbon (DOC) and (b) dissolved inorganic carbon (DIC) in the shallow zone (SZ) and deep zone (DZ); (c) Net production rates of DIC by different reactions; (d) Net production rates of DIC by biogenic (Resp-SZ + Resp-DZ) and geogenic (Carbonate-SZ and Carbonate-DZ) reactions; Net export rates (R e ) of (e) DOC and (f) DIC in SZ and DZ.Model results are from one example simulation (MC27).All rates are in mg C/m 2 /d.Modeled discharge (Q) is shown on the right y-axes of panels A and B for reference.DOC was primarily produced in and exported from the SZ.DIC was produced more in the SZ but exported more from the DZ, largely due to the higher concentrations of DIC in DZ and higher contributions of deeper flow throughout the year.Net DIC production was dominated by SZ reactions but included contributions from biogenic (Resp-SZ and Resp-DZ) and geogenic (Carbonate-SZ and Carbonate-DZ) reactions.Biogenic reactions showed a stronger seasonal trend with higher rates in warm summer, while geogenic reactions peaked throughout the year.DIC export rates from the SZ were generally flashy and more sporadic.Resp = respiration, Sorption = DOC sorption onto soils, Carbonate = carbonate weathering.

Figure 9 .
Figure 9. (a) A conceptual figure summarizing mean annual production and export rates (g/m 2 /yr) of dissolved organic carbon (DOC) and dissolved inorganic carbon (DIC) in the shallow zone (SZ) and deep zone (DZ) from water years 2016-2017; (b) Mean annual DIC production rates via different sources, where biogenic DIC is from respiration (Resp-SZ and Resp-DZ) and geogenic DIC is from carbonate weathering (Carbonate-SZ and Carbonate-DZ); (c) Mean production rates by season; (d) Mean export rates by season.All rates calculated as mean of seven simulations.Both DOC and DIC were produced more in the SZ than the DZ, but DOC was primarily exported from the SZ, whereas DIC was predominantly exported from the DZ.DIC was primarily produced via biogenic respiration; DIC production rates peaked in the spring and summer.Production rates for DOC peaked in the summer.Export rates of both DOC and DIC peaked in the spring.
The C ppt,i , C snow,i , C infil,i , C SF,i , C SZ,i, and C DZ,i are concentrations of solute i in precipitation, snow, infiltrating water, SF, SZ, and DZ, respectively.The reaction rates R SF,i , R SZ,i , and R DZ,i are those involving solute i in the three zones.Surface zone reactions were not implemented in this study, as existing data do not suggest substantial reactions at the ground surface.If a solute participates in more than one reaction, the overall R term is the summation of multiple reaction rates (Text S3 in Supporting Information S1).Water storage (V w ) and water fluxes (Q) are drainage-area-normalized values and have units of mm (volume per unit drainage area) and mm/d of water respectively.The primary species here
STEWART ET AL.