Temperature controls production but hydrology regulates export of dissolved organic carbon at the catchment scale

. Lateral carbon ﬂux through river networks is an important and poorly understood component of the global carbon budget. This work investigates how temperature and hydrology control the production and export of dissolved organic carbon (DOC) in the Susquehanna Shale Hills Critical Zone Observatory in Pennsylvania, USA. Using ﬁeld measurements of daily stream discharge, evapotranspiration, and stream DOC concentration, we calibrated the catchment-scale biogeochemical reactive transport model BioRT-Flux-PIHM (Biogeochemical Reactive Transport– Flux–Penn DOC concentration at high discharge. These results depict a conceptual model that the catchment serves as a producer and storage reservoir for DOC under hot and dry conditions and transitions into a DOC exporter under wet and cold conditions. This study also illustrates how different controls on DOC production and export – temperature and hydrological ﬂow paths, respectively – can create temporal asynchrony at the catchment scale. Future warming and increasing hydrological extremes could accentuate this asynchrony, with DOC production occurring primarily during dry periods and lateral export of DOC dominating in major storm events.


Introduction
Soil organic carbon (SOC) is the largest terrestrial stock of organic carbon, containing approximately 4 times more carbon than the atmosphere (Stockmann et al., 2013;Hugelius et al., 2014). Understanding the SOC balance requires the consideration of lateral fluxes in water, including dissolved organic and inorganic carbon (DOC and DIC, respectively), and vertical fluxes of gases such as CO 2 and CH 4 (Chapin et al., 2006). Both lateral and vertical fluxes influence SOC mineralization to the atmosphere (Campeau et al., 2019), although lateral fluxes are arguably less understood and integrated into Earth system models (Aufdenkampe et al., 2011;Raymond et al., 2016). Lateral fluxes from terrestrial to aquatic ecosystems are similar in magnitude to net vertical fluxes (Regnier et al., 2013;Battin et al., 2009), highlighting the importance of quantifying the controls on the lateral carbon (C) flux. In addition to its role in the global C cycle, DOC is an important water quality parameter that can mobilize metals and contaminants as well as imposing challenges for water treatment (Sadiq and Rodriguez, 2004;Bolan et al., 2011). DOC also regulates food web structures by acting as an energy source for heterotrophic organisms and interacts with other biogeochemical cycles (Malone et al., 2018;Abbott et al., 2016).
SOC decomposition and DOC production have been studied extensively (Abbott et al., 2015;Hale et al., 2015;Humbert et al., 2015;Lambert et al., 2013;Neff and Asner, 2001), yet the interactions between SOC and DOC and their response to climate change at catchments or larger scales remain unresolved (Laudon et al., 2012;Clark et al., 2010). Some regions have experienced long-term increases in DOC, potentially due to recovery from acid rain or climate-induced changes in temperature (T ) and hydrological flow (Laudon et al., 2012;Perdrial et al., 2014;Evans et al., 2012;Monteith et al., 2007), whereas others have observed decreases or no change (Skjelkvale et al., 2005;Worrall et al., 2018). Generally, the linkages among SOC processing, hydrological conditions, and DOC export or concentration remain poorly understood. Recent analyses indicate that the relationship between DOC concentration and discharge (C-Q) at stream outlets is primarily positive (Moatar et al., 2017;Zarnetske et al., 2018). Approximately 80 % of watersheds in the USA and France show a flushing C-Q pattern (i.e., the stream DOC concentration increases with discharge), whereas the rest shows dilution (decreasing DOC with discharge) or chemostatic behavior (negligible concentration change with discharge). These C-Q patterns generally correlate with catchment characteristics, including topography, wetland area, and climate characteristics, but it remains uncertain how hydrological and biogeochemical processes regulate SOC decomposition, DOC production, and DOC export (Jennings et al., 2010;Worrall et al., 2018). This gap in process understanding limits the integration of lateral carbon dynamics into projections of future ecosystem response to environmental change.
Stream DOC can be influenced by a variety of factors that control SOC decomposition and DOC production rates. DOC production generally increases as T increases; however, there may be multiple thermal optima, and the local rates can vary with SOC characteristics, soil type, and soil biota (Davidson and Janssens, 2006;Jarvis and Linder, 2000;Yan et al., 2018;Zarnetske et al., 2018). DOC production rates can exhibit low temperature sensitivity in highly weathered soils with a high clay content (Davidson and Janssens, 2006). They have also been shown to increase with soil water content in sandy loam soils (Yuste et al., 2007) and to have an optimum with a volumetric water content of approximately 0.75 in fine sands (Skopp et al., 1990). Because DOC export (or load) is the product of discharge and DOC concentration, it may differ from local DOC production rates in complex ways. For example, high T can produce a peak soil water DOC concentration but not necessarily stream concentration or export, due to temporal or spatial mismatches (D'Amore et al., 2015). These confounding factors present significant challenges to quantify the predominant mechanisms that regulate DOC production and export under varying environmental conditions.
One approach to understanding DOC production and export is the use of reactive transport models (RTM). These models integrate multiple production, consumption, and export processes, enabling the differentiation of individual and coupled processes (Steefel et al., 2015;Li, 2019;Li et al., 2017b). The use of RTMs complements statistical tools for the identification of influential factors (Correll et al., 2001;Herndon et al., 2015;Zarnetske et al., 2018). Historically, RTMs have been used in groundwater systems, where direct observations are particularly challenging (Kolbe et al., 2019;Li et al., 2009;. At the catchment scale, biogeochemical modules have been developed as add-ons to hydrological models. For example, a DOC production module was coupled to the HBV hydrological model using a static SOC pool that emphasized the influence of active-layer dynamics and slope aspect (Lessels et al., 2015). The INCA-C (Futter et al., 2007) and extended LPJ-GUESS  models have investigated the im- . The surface elevation was generated from lidar topographic data (criticalzone.org/shale-hills/data), whereas soil depths and SOC were interpolated using ordinary kriging based on field surveys with 77 and 56 sampling locations, respectively (Andrews et al., 2011;Lin, 2006). The SOC distribution in panel (c) is further simplified using the high, uniform SOC (5 % v/v) in swales and valley soils based on field survey information (Andrews et al., 2011). Swales and valley floor areas were defined based on surface elevation via field survey and a 10 m resolution digital elevation model (Lin, 2006). Additional sampling instrumentation is shown in panel (b), including six soil water sites (circles) and three soil T sites (squares). portance of land cover in determining DOC terrestrial routing and lateral transport. Terrestrial and aquatic carbon processes have also been integrated into the Soil and Water Assessment Tool (SWAT) to simulate aquatic DOC dynamics (Du et al., 2019). These modules typically simulate individual reactions without considering multicomponent reaction thermodynamics and kinetics.
In this context, the recently developed BioRT-Flux-PIHM model (BFP, Biogeochemical Reactive Transport-Flux-Penn State Integrated Hydrologic Model) fills an important gap by incorporating coupled elemental cycling, stoichiometry, and rigorous thermodynamics and kinetics (Bao et al., 2017;Zhi et al., 2019). We used the BFP to address the following question: how do hydrology and T interact to determine rates of DOC production and export at the catchment scale? We applied the BFP to a temperate forest catchment in the Susquehanna Shale Hills Critical Zone Observatory (SSHCZO). This small catchment (< 0.1 km 2 ) has gentle topography with a network of shallow depressions or swales that have high SOC and deep soils (detailed in Sect. 2). It is underlain by one type of lithology (shale) and land use (forest), providing a useful test bed to evaluate biogeochemical and hydrological functions . Previous lab and field work have identified non-chemostatic C-Q patterns of DOC at SSHCZO that are attributable to differences in the hydrologic connectivity of organic-rich soils during different flow conditions (Andrews et al., 2011;Herndon et al., 2015). SSHCZO has spatially extensive and highfrequency measurements of soil properties, hydrology, and biogeochemistry . These data facilitate detailed benchmarking of the BFP model and evaluation of processes controlling DOC production and export. We expected that T and soil moisture would drive DOC production in the soil, whereas DOC export and, thus, C-Q patterns would be most related to hydrological connectivity. Therefore, we predicted that DOC production and export might be asynchronous (i.e., not occurring at the same time) because they respond differently to changes in T and hydrology. Al-though soil respiration is an important process, this study focuses on the net production and export of DOC.

Study site: a small catchment with an intermittent stream
The Shale Hills catchment is a 0.08 km 2 , V-shaped, firstorder watershed with an intermittent stream in central Pennsylvania. It is forested with coniferous trees and is situated on the Rose Hill shale formation. The annual mean air T is 9.8 ± 1.9 • C (± SD) and the annual mean precipitation is 1029 ± 270 mm over the past decade. The watershed is characterized by large areas of swales and valley floors with deep and wet soils (Fig. 1b). These lowland soils contain more SOC (∼ 5 % v/v) than the hillslopes and uplands (∼ 1 % v/v; Fig. 1c). Soil water DOC samples were collected using lysimeters with a diameter of 5 cm installed at 10 or 20 cm intervals from the soil surface to a depth of hand-auger refusal, which varied from 30 to 160 cm depending on soil thickness. There were a total of six sampling locations (Fig. 1b), including three at the south planar sites -valley floor (SPVF), midslope (SPMS), and ridgetop (SPRT) -and three at the swale sites -valley floor (SSVF), mid-slope (SSMS), and ridgetop (SSRT). No soil water DOC samples were collected on the north side of the catchment. Stream water DOC samples were collected daily in glass bottles at the stream outlet weir. All soil water and stream water DOC samples were filtered to 0.45 µm using Nylon syringe filters and were analyzed with a Shimadzu TOC-5000A analyzer (detailed in Andrews et al., 2011). Real-time soil T (every 10 min) was measured at the ridgetop, mid-slope, and valley floor (squares in Fig. 1b) using automatic monitoring stations at depths of ∼ 0.10, 0.20, 0.40, 0.70, 0.90, 1.00, and 1.30 m (Lin and Zhou, 2008). Figure 2. A schematic representation of major processes in the catchment reactive transport model BFP (BioRT-Flux-PIHM). Stream discharge Q includes surface runoff Q S , soil water interflow (lateral flow) Q L , and groundwater flow Q G . In the vertical direction, soil pores are not saturated with water in the shallow unsaturated zone and water flows vertically until it reaches the saturated zone where water forms interflows and moves laterally to the stream. Soil water total storage S T is the sum of water in the unsaturated (S u ) and saturated zones (S s ). Some water also recharges further into deeper groundwater. Within the soil zone, SOC decomposes and releases DOC, which also sorbs onto the soil surface to become ≡ XDOC.

The BioRT-Flux-PIHM (BFP) model
BFP is the catchment reactive transport model of the general PIHM (Penn State Integrated Hydrologic Modeling System) family of code . The code includes three modules ( Fig. 2): the surface hydrological module PIHM, the land surface module Flux, and the multicomponent reactive transport module BioRT (Biogeochemical Reactive Transport). The code has been applied to simulate conservative solute transport, chemical weathering, surface complexation, and biogeochemical reactions at the catchment scale (Bao et al., 2017;Zhi et al., 2019;Li, 2019). Here, we only introduce the salient features that are relevant to this study; readers are referred to earlier publications for further details. Flux-PIHM separates the subsurface flow into active interflow in shallow soil zones and groundwater flow deeper than the soil-weathered rock interface. Note that this "deeper groundwater" is the groundwater that actively interacts with the stream and shallow layer, not necessarily the water in the deep groundwater aquifer. The PIHM module simulates hydrological processes including precipitation, infiltration, surface runoff Q S , soil water interflow (lateral flow) Q L , and discharge Q (Fig. 2). The Flux module simulates processes including solar radiation and evapotranspiration. Flux-PIHM calculates water variables (e.g., water storage, soil moisture, and water table depth) in unsaturated and saturated zones and assumes a no-flow boundary at the soil-bedrock interface with high permeability contrast. In this version of Flux-PIHM, the deeper groundwater flow Q G is a separate input to the stream and is decoupled from the shallow soil water. This is supported by field data that show negligible seasonal varia-tion in groundwater chemistry (Jin et al., 2014;Thomas et al., 2013;Kim et al., 2018). The Q G is estimated using conductivity mass-balance hydrograph separation (Lim et al., 2005).
The BioRT module takes in water calculated at each time step to simulate reactive transport processes. BFM discretizes the domain into prismatic elements and uses a finite volume approach based on mass conservation. The mass conservation governing equation for the reactive transport of a single solute m is as follows: where i and j represent the grid block i and the neighboring grid j ; the subscript x distinguishes between flow in the unsaturated zone (infiltration and recharge) and saturated zone (recharge and lateral flow); V is the total bulk volume (m 3 ) of each grid block; S w is the soil moisture (m 3 water m −3 pore volume); θ is porosity; C is the aqueous species concentration (mol m −3 water); t is time (s); N is the index of elements sharing surfaces; A is the grid interface area (m 2 ); D is the diffusion/dispersion coefficient (m 2 s −1 ); l is the distance (m) between the center of two neighboring grid blocks; q is the flow rate (m 3 s −1 ); r m is the kinetically controlled reaction rates (mol s −1 ) involving species m, which is the DOC production rate from SOC decomposition at the grid block i; and np is the total number of independent solutes.

DOC production and sorption
In the model, DOC is produced by the decomposition of SOC via the kinetically controlled reaction SOC(s) → DOC. With abundant SOC and O 2 in soils serving as electron donors and acceptors, a typical dual Monod kinetics can be simplified into zero-order kinetics with additional temperature and soil moisture dependence: where r p is the local DOC production rate in individual grids (r m in Eq. 1, m is DOC); k is the kinetic rate constant of net DOC production (10 −10 mol m −2 s −1 ) (Zhi et al., 2019;Wieder et al., 2014); and A is a lumped "surface area" (m 2 , (2.5 × 10 −3 m 2 g −1 ) × (g of SOC mass)) that quantifies SOC content and biomass abundance (Chiou et al., 1990;Kaiser and Guggenberger, 2003;Zhi et al., 2019). The functions f (T ) and f (S w ) describe the rate dependence on soil T and moisture, respectively. f (T ) follows a widely used Q 10based formation: f (T ) = Q |T −10|/10 10 , where Q 10 quantifies the rate increases with T , with the superscript 10 referring to a T value of 10 • C (Davidson and Janssens, 2006). Q 10 in the base case is set at 2.0, within the typical range of 1.2-3.8 for forest ecosystems (Liu et al., 2017). The f (S w ) has the form f (S w ) = (S w ) n in the base case, where n is the saturation exponent with a value of 1.0, which is within the typical range of 0.75-3.0 for most soils (Yan et al., 2018;Hamamoto et al., 2010). The dependence of production rates on soil T and moisture have been described using multiple forms in existing studies (Davidson and Janssens, 2006;Yan et al., 2018) and will be further explored via sensitivity analysis, as detailed in Sect. 2.6. The SOC content typically decreases with depth (Billings et al., 2018;Bishop et al., 2004), although the specific pattern may vary with soil texture, landscape position, vegetation, and climate (Jobbagy and Jackson, 2000). The depth function of SOC at Shale Hills has been observed to be exponential (Andrews et al., 2011), which is typical of many soils (Billings et al., 2018;Currie et al., 1996). To take this into account, we use the equation where C d is SOC at depth z below the surface; C 0 is the SOC level at the ground surface, and b m quantified the decline rate with depth, which is set here to a value of 0.3 (Weiler and McDonnell, 2006). DOC produced from SOC can also sorb on soils via the following reaction: ≡ X+DOC ↔≡ XDOC, where ≡ X and ≡ XDOC represent the functional group without and with sorbed DOC, respectively . This reaction is considered fast and is thermodynamically controlled with an equilibrium constant K eq that links the activity (here approximated by concentrations) of the three chemicals via K eq = [≡XDOC] [≡X] [DOC] . The DOC concentrations calculated from Eq. (1) were used to establish the concentrations of ≡ X and ≡ XDOC. The K eq value represents the thermodynamic limit of the sorption, i.e., the sorption affinity of the soil for DOC.
It depends on temperature but also on soil properties such as the clay content and the abundance of iron oxides (Kaiser et al., 2001;Conant et al., 2011). A K eq value of 10 0.2 was obtained by fitting the stream and soil water DOC data (detailed in Sect. 2.4). The sum of [≡ X] and [≡ XDOC] represents the soil sorption capacity. A value ranging from 4.0 × 10 −5 to 6.0 × 10 −5 mol g −1 soil was used for Shale Hills (Jin et al., 2010;Li et al., 2017a) depending on the mineralogy in different zones of the catchment.

Domain setup
BFP is a model with full discretization in the horizontal direction and partial discretization in the vertical direction with three layers: ground surface, unsaturated, and saturated zones. Although a new version of BFP explicitly includes a groundwater zone, it was not released in time for this work, so the groundwater fluxes were estimated separately. The study watershed was discretized into 535 prismatic land elements and 20 stream segments using PIHMgis (http://www.pihm.psu.edu/pihmgis_home.html, last access: 11 February 2020), a GIS interface for BFP. The land elements are unstructured triangles with mesh sizes varying from 10 to 100 m. The simulation domain was set up using national datasets: the USGS National Elevation Dataset for topography; the National Land cover Database for vegetation distribution; the National Hydrography Dataset for water drainage network; the North American Land Data Assimilation Systems Phase 2 (NLDAS-2) for hourly meteorological forcing; and the Moderate Resolution Imaging Spectroradiometer (MODIS) for leaf area index. In addition, extensive characterization and measurement data at Shale Hills were used to define soil depth and soil mineralogical properties such as surface area and ion exchange capacity that are heterogeneously distributed across the catchment (Andrews et al., 2011;Lin, 2006;Jin and Brantley, 2011;Jin et al., 2010;Shi et al., 2013) (http://criticalzone.org/shale-hills/data/, last access: 11 February 2020). Other soil matrix properties include conductivity, porosity, and van Genuchten parameters. Soil macropores such as cracks, fractures, and roots can generate preferential flows. Their properties are represented using the area macropore fraction, depth, and conductivities. They are parameterized based on values quantified in previous studies at Shale Hills (Shi et al., 2013;Lin, 2006), as shown in Fig. S1 and Table S1 in the Supplement.
Based on field measurements, the SOC content in swales and valley areas is relatively high (Andrews et al., 2011) and was set at 5 % (v/v solid phase) compared with 1 % in the rest of the catchment (Fig. 1c). The clay minerals were set at 23 % (v/v solid phase) along the ridgetop and 33 % at the valley floor (Jin et al., 2010;Li et al., 2017a). The input DOC concentrations in rainfall and groundwater (below soils) were set at reported medians of 0.6 and 1.2 mg L −1 , respectively (Andrews et al., 2011;Iavorivska et al., 2016), as high-frequency DOC observations were not available. The initial DOC concentration in soil water was set at 2.0 mg L −1 , which was the average concentration from the six field sampling locations in Fig. 1.

Model calibration
We used stream (daily) and soil pore water (biweekly) DOC concentration data from April to October 2009 for model calibration and the year 2008 as spin-up until a "steady state" for both water and DOC was reached. The "steady state" here refers to a state where the inter-annual difference between stored mass within the catchment is less than 5 % of the total mass. The water input is precipitation, and its output is ET and discharge. The DOC mass input is from rainfall, groundwater, and production, and the DOC output is the export load at the stream outlet. The model performance was evaluated using the monthly Nash-Sutcliffe efficiency (NSE) (Nash and Sutcliffe, 1970) that quantified the residual variance of modeling output compared to measurements. The general satisfactory range for monthly average outputs for hydrological models is NSE > 0.5 (Moriasi et al., 2007), and we used similar standards for biogeochemical solutes (Li et al., 2017a). To reproduce the DOC data, we first set the SOC surface area A using a literature range of 10 −3 -10 0 m 2 g −1 (Zhi et al., 2019;Chiou et al., 1990;Kaiser and Guggenberger, 2003). We also set K eq using a literature range of 10 0 -10 1 (Oren and Chefetz, 2012;Ling et al., 2006). Once the simulated output captured the temporal trend of data, we refined Q G based on the estimation from hydrograph separation ( Fig. S2) to capture the peaks of stream DOC concentration, especially under low-discharge periods. Because not all soils are in contact with water, the calibrated surface area represents the effective solid-water contact area, and is orders of magnitude lower than the reported SOC surface areas from laboratory experiments (Kaiser and Guggenberger, 2003). The calibrated hydrological parameters are mostly from Shi et al. (2013), except groundwater estimation. Groundwater estimates were based on Li et al. (2017a) and further refined using conductivity mass-balance hydrograph separation (Lim et al., 2005) and then by reproducing the stream DOC concentration. In other words, stream and groundwater chemistry data together helped constrain the groundwater flow.
2.6 Quantification of water and DOC dynamics 2.6.1 Hydrological connectivity Saturated soil water storage calculated from the model was used to quantify hydrological connectivity I cs /Width. With "Width" defined as the average width of catchment in the direction perpendicular to the stream (230 m), the term I cs /Width quantifies the average proportional width of the catchment connected to the stream (e.g., I cs /Width = 0.10, 0.35, and 0.70 in Fig. S3). Depending on the catchment ge-ometry and the extent of connectivity, I cs /Width may vary from 0 to 1.0. A high I cs /Width value (i.e., high hydrological connectivity) indicates that a large catchment area is connected to the stream. To determine whether two grids are hydrologically connected, the spatial distribution of saturated water storage was used to calculate connectivity following the equation I cs = ∞ 0 τ (h)dh and an algorithm in the literature (Allard, 1994;Western et al., 2001;Xiao et al., 2019). Here τ (h) is the probability of two grid blocks being connected at a separation distance of h. Two grids are considered "connected" if they are joined by a continuous flow path and have saturated storages exceeding the threshold of the 75th percentile of saturated storage (over the whole catchment). Note that I cs /Width here only quantifies the hydrological connectivity in soils and does not reflect the groundwater in shallow aquifers below the soil-bedrock interface.

DOC concentration-discharge relationships
At the catchment scale, we differentiate the DOC production rates and export rates. The production rate R p is the sum of the local DOC production rate r p in individual grid blocks (Eq. 2) across the whole catchment. The export rate R e is the product of discharge and the DOC concentration at the stream outlet. Total stored DOC is the difference between stream output and input from production, rainfall, and groundwater. The DOC input from the rainfall R r (mg d −1 ) is the precipitation rate (m d −1 ) times the rainfall DOC concentration (6.0 × 10 −4 mg m −3 = 0.6 mg L −1 × 10 −3 L m −3 ) and the catchment drainage area (m 2 ). The DOC input from groundwater R g (mg d −1 ) is the total groundwater influx (flow rate) times the groundwater DOC concentration (1.2 mg L −1 ).
C-Q patterns were quantified using two complementary approaches: the power law equation C = aQ b (Godsey et al., 2009) and the ratio of the coefficients of variation of the DOC concentration and discharge CV [DOC] CV Q (Musolff et al., 2015). The slope of the power law equation b does not account for the goodness of fit of the C-Q pattern itself. For example, a slope of b = 0 would be considered chemostatic (i.e., relatively small variation of concentration compared with discharge), although high variability in solute concentrations would in fact reflect a chemodynamic behavior (i.e., solute concentrations are sensitive to changes in discharge). We considered two general categories based on these metrics (Musolff et al., 2015): if b values fell between −0.2 and 0.2 and CV [DOC] CV Q 1, C-Q patterns were considered chemostatic; values of |b| > 0.2 or CV [DOC] CV Q ≥ 1, indicated a chemodynamic behavior. In the chemodynamic category, values of b > 0.2 indicate flushing, whereas values of b < −0.2 indicate dilution. We used the MATLAB curve-fitting toolbox to obtain the best fit model parameters. . Temporal dynamics of (a) daily precipitation, stream discharge Q, and evapotranspiration ET on the arithmetic scale; (b) stream discharge Q, soil water interflow Q L , and groundwater Q G on a logarithmic scale with soil T on an arithmetic axis (right); (c) soil water storage S T (unsaturated water storage S u + saturated water storage S s ) and hydrological connectivity I cs /Width. The yellow dots in panel (b) represent the average soil T from three sampling locations (square symbols in Fig. 1b) with the shading reflecting variation in measurement. Q was highly responsive to intense precipitation events in spring and winter. Note high soil T , high ET, low S s , and lowI cs /Width during July-August 2009. Stream discharge was primarily comprised of Q L , except in July-October when the relative contribution of Q G increased.

Sensitivity analysis
We used a sensitivity analysis to explore the influence of soil T and moisture in the DOC production kinetics. The Q 10 in f (T ) = Q |T −10|/10 10 was explored using a minimum value of 1.0 (i.e., no dependence on T ) and a maximum value of 4.0 (Davidson and Janssens, 2006) (Fig. S4a), i.e., f (T ) = 1 and f (T ) = 4 |T −10|/10 . The rate dependence on soil moisture was explored using the base case f 1 (S w ) = (S w ) n (increase behavior), and three additional functions (f 2 , f 3 , and f 4 ) representing the most commonly observed forms (Fig. S4b), including decrease behavior, constant behavior, and threshold behavior (Gomez et al., 2012;Yan et al., 2018): the decrease behavior function was the constant behavior function was and the threshold behavior function was The constants in Eqs.
(3)-(5) were selected to ensure similar averages of f (S w ) across the whole S w range such that trajectories rather than absolute values of f (S w ) were compared (Fig. S4b). The sensitivity of DOC sorption onto soils was tested using K eq values of 0 (no sorption), 10 0.5 and 10 1.0 .
The sensitivity of C-Q patterns and R e to changes in groundwater was also tested with groundwater flow contribution and DOC concentration. The groundwater flow rates were varied from negligible (Q G = 0) to 2.5 times those of the base case (Q G = 3.3×10 −4 and 1.0×10 −4 m d −1 for the wet and dry periods, respectively). The corresponding fractions (Q G /Q) of groundwater flow to the total annual discharge for the two cases were 0 % and 18.8 %, respectively. The groundwater DOC concentration (DOC GW ) was varied by 2 orders of magnitude (0.12 and 12.0 mg L −1 ). Results from these analyses were compared with the base case, in which the groundwater contributed to 7.5 % of the total annual streamflow at 1.2 mg L −1 .

Water dynamics
The total precipitation from 1 April 2009 to 31 March 2010 was 1130 mm. Stream discharge was highly responsive to intense precipitation events and was high (∼ 10 −2 m d −1 ) in spring and fall compared with summer with high soil T and high ET (∼ 10 −5 m d −1 ). The model captured the temporal dynamics of daily discharge, ET, and soil T with NSE values  Fig. 1b, including three planar (panels b-d) and three swale locations (panels e-g). The mean ±SD for each location was calculated based on measurements at different depths with 10 or 20 cm intervals from the soil surface down to a depth of hand-auger refusal. of 0.68, 0.72, and 0.62, respectively (Fig. 3a, b). The model estimated that 47.5 % of annual precipitation contributed to discharge, whereas the rest contributed to ET. The stream discharge has three components: surface runoff Q S , soil water interflow Q L (lateral flow), and groundwater flow Q G from the shallow subsurface that interacts with the stream (Fig. 2). On average, lateral flow Q L is about 90.2 % and surface runoff Q S is about 2.3 %. Following the conductivity mass-balance hydrograph separation (Lim et al., 2005), Q G was estimated to be 1.3 × 10 −4 and 4.0 × 10 −5 m d −1 for the wet and dry periods (August-September), which is equivalent to 6.9 % and 42.2 % of average stream discharge at the corresponding times, respectively. Overall Q G accounted for ∼ 7.5 % of the annual Q, similar to previously reported values (Li et al., 2017a;Hoagland et al., 2017). In the dry months from August to September, the stream was almost dry with no visible flow, and the relative contribution of groundwater to discharge was comparable to that of Q L (Fig. 3b). The unsaturated water storage S u was typically more than 10 times larger than the saturated storage S s such that the S T and S u curves almost overlapped (Fig. 3c). S s was negligible in the dry period (close to 0 m), contributing negligibly to the stream. Hydrological connectivity (I cs /Width) covaried with S s but showed significant temporal fluctua-tions. High summer ET drove the catchment to drier conditions, thereby decreasing the connectivity to the stream.

Temporal patterns of DOC concentrations
The model captured the general trend of stream DOC (NSE of 0.55 for the monthly DOC concentration; Fig. 4). Under dry conditions (e.g., Q < 1.0×10 −4 m d −1 ), Q G contributed substantially to Q (∼ 32 %-71 %; Fig. 3), and the stream DOC concentration reflected the mixing of groundwater and soil water (Fig. 4a), with a contribution from groundwater DOC of 7 %-17 %. Under wet conditions, the stream DOC concentration overlapped with the soil water DOC concentration (light blue line in Fig. 4). Only ∼ 1 %-8 % of stream DOC was sourced from groundwater at these times.
The temporal dynamics of soil water data showed relatively small temporal variation compared with stream DOC (Fig. 4b, c, d, e, f, g), and local soil pools were not always hydrologically connected to the stream. The simulated soil water DOC captured this small-variation trend with acceptable overall model performance (i.e., NSE > 0.5), although the goodness of fit was lower in some locations, e.g. The soil DOC and r p were high in swales and the valley that had a relatively high soil water and SOC content (Fig. 1c). Although water content in August was relatively low compared with May and October, high soil T led to high r p , with most DOC production and accumulation in zones that were disconnected from the stream. organic carbon content. Although the model explicitly considered spatial heterogeneities such as topography and soil properties, averaged values represented grid sizes from 10 to 100 m, and this local scale was large compared with the field sampling size (e.g., lysimeters with a diameter of 5 cm). Geochemical processes are sensitive to local properties, including SOC %, SOC surface area, and sorption sites, and the representation of these properties was based on a few measurements that were only coarsely defined as ridgetop, mid-slope, and valley floor.

Spatial patterns and mass balance
Spatial patterns vary between May (wet), August (dry), and October (wet after dry) (Fig. 5). In May, the average soil T was around 12 • C with small spatial variations (< 3 • C). Most flow-convergent areas (valley areas and swales) were well connected to the stream and had a high water content (Fig. 5b, c). The distribution resembles that of SOC (Fig. 1c) and water content (Fig. 5b), with a high r p and soil water DOC concentration in swales and valley. Low r p in relatively dry planar hillslopes and uplands led to a low soil water DOC concentration. In August, the average soil T increased to around 20 • C. The hydrologically connected zones shrank to the immediate vicinity of the stream, but r p increased 2- Figure 6. Temporal dynamics of DOC storage, influent rate (rainfall R r , groundwater R g , production R p ), and outflow rate (effluent R e ) at the catchment scale. The stored DOC mass (dark red line) was calculated as follows: (DOC influent rate − outflow rate) × time.
The temporal R e dynamics mostly followed the trend of discharge (black line, top panel), whereas R p mostly followed the trend of soil T (orange line, top panel).
fold from May. The simulated soil water DOC concentration increased by a factor of 2 across the whole catchment, especially in hillslope and uplands on the north side, because the DOC produced was trapped in low soil moisture areas that were not hydrologically connected to the stream. This indicates that DOC samples collected on the south side may not represent the DOC dynamics of the entire catchment, especially in the summer and fall dry months. In October, r p decreased as soil cooled down, but increased precipitation and decreased ET expanded the hydrologically connected zones beyond swales and valley areas (Fig. 5c), promoting the desorption and the flushing of stored DOC. The soil water DOC concentration, however, remained high because of the large store of sorbed DOC produced during the antecedent dry times. Figure 6 shows the catchment-scale DOC production and export rates and mass balance. Generally, the daily R p (5.1 × 10 5 mg d −1 ) was greater than the daily R r from rainfall (1.6× 10 5 mg d −1 ) or groundwater R g (1.2 × 10 4 mg d −1 ). During storm events, R r occasionally exceeded R p . R p was generally high in summer, despite low water storage. The export rate R e did not follow the temporal patterns of the total input rate (R p + R r + R g ) or R p . Instead, it primarily followed the discharge patterns: large rainfall events exported disproportion- Relationship of daily discharge (Q) with soil water storage S T , connectivity (I cs /Width), the catchment-scale DOC export rate R e , and the DOC production rate R p . At low Q, the stream water transitioned from organic-poor groundwater to organic-rich water from the valley floor and swales, leading to a flushing (positive) pattern. At higher Q, the stream water shifted from organic-rich soil water from swales and valley areas to lower DOC water from planar hillslopes and uplands, decreasing the stream DOC concentration and resulting in a dilution C-Q pattern. R e increased by 2 orders of magnitude with increasing Q, whereas R p varied within an order of magnitude. Figure 8. The catchment-scale DOC production rate R p and export rate R e as a function of (a) soil T , (b) soil water storage S T , and (c) hydrological connectivity (I cs /Width). Cross symbols are daily values in the base case. R p increased with soil T and decreased slightly with S T and connectivity. In contrast, R e increased with S T and connectivity but decreased with soil T . R e tended to decrease with soil T in the hot, dry summer due to low discharge during that period. ally high DOC, plummeting the DOC mass within the catchment. From the wet to dry period, as water levels dropped, DOC accumulated within the catchment (Fig. 5e, May to August). During the dry-to-wet transition, as the catchment became wetter, the contributing areas expanded to the uplands and the DOC was flushed out, reducing the overall DOC soil pool to much lower values (Fig. 5e, August-October). The DOC mass storage increased by 1.8 × 10 6 mg over the year, which was about 1.0 % of the overall DOC production, indicating a general mass balance at the catchment scale.

C-Q patterns and rate dependence
The C-Q relationships showed a slightly positive correlation at low Q followed by a negative correlation at higher Q (Fig. 7a). The simulated C-Q relationship captured this trend but overestimated the positive relationship at low Q.
The simulated C-Q relationships showed a general dilution behavior with the C-Q slope b = −0.23 and CV [DOC] CV Q = 0.22, which was consistent with the general pattern exhibited in the field data (Fig. 7a). This C-Q pattern can be explained by the dynamics of different water sources with different DOC contributing to the stream. At low discharges (< 1.8×10 −4 m d −1 ) with small water storage (0.25-0.28 m) and connectivity (I cs /Width < 0.1) (Fig. 7b), the stream DOC was a mix of organic-poor groundwater and organicrich swales and valley floor zones. As connectivity and discharge increased and the stream expanded, the contribution of organic-rich swales increased, elevating the DOC concentration to its maximum. Under even wetter conditions with connectivity exceeding 0.1, the contribution from planar hillslopes and uplands with a lower DOC concentration increased, diluting the organic-rich DOC from swales and val- Figure 9. Sensitivity analysis of temporal DOC rates for (a) soil temperature f (T ) and (b) soil moisture f (S w ). A varying Q 10 value in f (T ) had a larger influence on R p than varying f (S w ). Neither f (T ) nor f (S w ) had a significant influence on R e . Instead, R e mostly followed the temporal trend of discharge, indicating the predominant control of hydrological conditions. Figure 10. Sensitivity analysis of the sorption equilibrium constant K eq on (a) R p and R e and on (b) DOC sorbed on soils averaged at the catchment scale. High K eq led to more DOC sorbed on soils and, therefore, lower R e . However, R e showed similar temporal patterns regardless of K eq . ley areas. Daily R e correlated positively with S T , hydrological connectivity, and Q, and increased by 2 orders of magnitude as Q rose by 3 order of magnitude. The variation of daily R p with Q was small (10 5 -10 6 mg d −1 ) compared with that of R e (Fig. 7b). Values of R p depended more on soil T than on soil water storage and hydrological connectivity (I cs /Width) (Fig. 8). In contrast, R e increased with soil water storage S T but notably decreased with soil T (> 17 • C) due to the low discharge during the hot and dry summer.

Sensitivity analysis
3.5.1 Control of temperature, soil moisture, and sorption Higher Q 10 values in f (T ) led to more pronounced seasonality in R p (Fig. 9a). The R p for Q 10 = 4.0 was more than 4 times higher than that of Q 10 = 1.0 in summer, and much lower in winter with low soil T (< 10 • C). In contrast, the temporal pattern of R e almost overlapped at different Q 10 values, and it mostly followed the discharge dynamics (black line in Fig. 9). Daily R p varied only slightly (within 15 %) with different f (S w ) (Fig. S4b), while R e showed very little change (Fig. 9b). Although we varied Q 10 from 1.0 to 4.0 in f (T ), it is worth noting that varying the kinetic rate constant, SOC surface area, volume fraction, and biomass amount could have similar effects (not shown here) because they are all multiplied in Eq.
(2). Simulations showed that strong DOC sorption (K eq = 10 1.0 ) did not change R p but lowered the stream DOC concentration and resulted in smaller R e (Fig. 10a). DOC sorption had little impact on R p , but strong sorption decreased the magnitude of R e by 10 %-69 %. The sorbed DOC concentration differed by more than a factor of 3, with more sorbed DOC with larger K eq values (Fig. 10b). Large amounts of sorbed DOC persisted until early fall, when large rainfall events flushed out sorbed DOC and reduced DOC storage (Fig. 6). This means that catchments can store large quantities of DOC, although the specific amount of DOC stored depends on sorption capacity.
Varying DOC production kinetics did not change the overall C-Q patterns, although the magnitude of overall dilution changed slightly in cases with different f (T ) and K eq (Fig. 11). High Q 10 values in f (T ) led to less dilution, due to more accumulated soil DOC in the dry period (low discharge) and, thus, more DOC flushing as discharge increased in the dry-to-wet period. High K eq resulted in less dilution as the higher sorption capacity acts as a stronger buffer to compensate for the concentration variations. , and (c) sorption equilibrium constants K eq for the two extreme cases. The C-Q patterns were similar in all cases, although the extent of dilution slightly changed. This indicates potential factors other than reaction kinetics and thermodynamics that regulate C-Q patterns. Figure 12. Sensitivity analysis of groundwater on rates (R p and R e ) and C-Q relationships: (a) scenarios with a different groundwater volume contribution (%) to stream discharge and (b) scenarios with a different groundwater DOC concentration (DOC GW ). DOC GW and GW (Q G /Q) in the base case were 1.2 mg L −1 and 7.5 %, respectively. "2.5 GW" in panel (a) represents the case with 2.5 times Q G compared with the base case. Increases in the relative groundwater contribution lowered R e and shifted the C-Q pattern from an overall dilution pattern to an overall flushing pattern; changing DOC GW had negligible influence on the DOC rates and C-Q patterns.

Groundwater control on DOC export
As shown in Fig. 12, changing the groundwater volume contribution to stream (GW) had more significant impacts than changing the groundwater DOC concentration (DOC GW ), especially at low discharges (Q < 1.8 × 10 −4 m d −1 ). Increasing the GW contribution from no GW to 2.5 GW (i.e., 18.8 %) lowered stream DOC at low discharges, shifting the C-Q pattern from overall dilution (or a chevron pattern) to overall flushing (or flushing until stable). More specifically, the threshold that separated distinct phases of these segmented C-Q responses (Fig. 12a2) shifted from Q = 1.8 × 10 −4 to about 1.0 × 10 −3 m d −1 . This reflects the rela-tive groundwater contribution to streamflow for each case. In contrast, varying the groundwater DOC concentration (DOC GW ) by 2 orders of magnitude while keeping the same groundwater contribution (GW) did not change C-Q pattern. Figure 13 summarizes the annual total R p and R e in all sensitivity test scenarios. Annual R p was more sensitive to T than to S w or sorption thermodynamics. Annual R e was less sensitive to T variation, although it also increased with Q 10 because a higher production led to higher DOC export. Annual R p also depended on f (S w ), with the threshold function f 4 (S w ) (Sect. 2.6) having the highest production rates. However, R e did not follow the trend of R p (Fig. 13b). Generally, under the same hydrological conditions, a doubling Hydrol. Earth Syst. Sci., 24, 945-966, 2020 www.hydrol-earth-syst-sci.net/24/945/2020/ Figure 13. Total annual R p (red) and R e (blue) under two groundwater volume contribution conditions (Q G /Q = 7.5 % and 18.8 %) for three different variables: (a) soil T , (b) soil moisture, and (c) sorption equilibrium K eq . R p was not influenced by a deeper groundwater contribution, so there is only one curve for each variable. R p has higher dependence on temperature than on soil moisture function form and sorption capacity.
of R p only led to about a 50 % increase in R e . Higher sorption affinity (higher K eq ) did not change production rates but could reduce DOC export by about 30 % due to high DOC storage in soils. High relative groundwater inputs (18.8 % versus 7.5 %) lowered R e in all scenarios because more water came from deeper organic-poor groundwater.

Discussion
This study revealed that DOC production was primarily regulated by temperature, but the lateral export of DOC was controlled by hydrological conditions. This work contributes to the growing body of research concluding that lateral carbon flux is determined by water routing and hydrological connectivity and only secondarily by biological activity (Zarnetske et al., 2018). Although soil respiration and vertical CO 2 fluxes are closely related, this work focuses on the net production and export of DOC because it has been studied and understood to a much lesser extent than soil respiration (Tank et al., 2018). To better appreciate the relative importance of land-water-atmosphere carbon fluxes, future research needs to fully integrate lateral DOC fluxes in concert with vertical fluxes of CO 2 across terrestrial and freshwater ecosystems.

DOC production
The DOC production rate R p depends more on T than on water storage or soil moisture. This finding is expected, as DOC production is biologically mediated and, thus, influenced by temperature and the metabolic dependence on temperature (Gillooly et al., 2001). Although the local-scale soil moisture varies from 0.40 at the ridgetop to 0.70 in swales and riparian zones (Fig. 5b), the averaged catchment-scale soil moisture has relatively small variations across different seasons in this temperate humid catchment (0.46 to 0.56 on average over the whole catchment), especially compared with places where water is limited and soil moisture can drop below 0.15 (Korres et al., 2015). This small variation is due to the capability of the shale-derived, clay-rich soils at Shale Hills to hold water (Xiao et al., 2019). The influence of soil moisture on DOC production is likely higher in catchments with more pronounced seasonal changes and more fluctuations in soil moisture. This work also suggests that catchment-scale (R p ) and local-scale (r p ) production rates have different controls. The rate law used at the local scale is measured at relatively small scales, i.e., 0.1-2.0 m in soil pedons (Bauer et al., 2008;Yan et al., 2016). Our results show that even when the rate law with an optimum soil moisture was used at the local scale (f 4 (S w ) in Fig. S4b), the catchment-scale rates do not exhibit maximum rates at an "optimal" soil moisture (Fig. 8), indicating different controls at the local scale versus the catchment scale. In addition, due to the spatial heterogeneities of T , soil moisture, and SOC content, the temporal variations of R p and r p may be not consistent. The daily R p spanned less than an order of magnitude with its maximum in the dry, hot summer and its minimum in the wet, cold winter and spring (Fig. 6). Local-scale r p exhibited similar temporal dynamics but varied by more than 2 orders of magnitude, with rapid production mostly in "hot spots" (swales and riparian zones) with persistently high water and SOC content (Fig. 5). Note that the local-scale rate laws are often used directly at the catchment scale or at even larger scales (Crowther et al., 2016;Conant et al., 2011;Fissore et al., 2009;Moyano et al., 2012). This work suggests that although local-scale rate laws have been developed extensively, direct extrapolation of rates from local to catchment scales can be misleading. This speaks to the importance of understanding controls on biogeochemical transformation rates and developing reaction rate theories at the catchment scale for regional-scale and global-scale simulations.
The simulations here suggest that DOC storage depends on the sorbing capacity of soils such that clay content and the presence of organo-mineral aggregates might play a role in mediating DOC dynamics (Lehmann et al., 2007;Cincotta et al., 2019).
DOC production depends on catchment size, hydrogeologic structure, vegetation, and climatic setting. Geomorphological and ecological processes have been shown to cogenerate systematic differences in the vertical and lateral distribution of SOC and plant biomass, with a greater concentration of organic carbon in the valley floor than in hillslopes in some catchments (Piney et al., 2018;Temnerud et al., 2016;Campeau et al., 2019;Thomas et al., 2016) and enriched SOC in the uplands in other catchments (Herndon et al., 2015). These differences may explain the large variation of stream DOC in catchments under diverse climate conditions (Moatar et al., 2017). The median stream DOC at Shale Hills is relatively high (10.0 mg L −1 ), compared with 3.0 mg L −1 in temperate humid catchments in Germany (Musolff et al., 2018), 4.1 mg L −1 in some UK catchments with oceanic climate (Monteith et al., 2015), and 4.5 mg L −1 in France (Moatar et al., 2017). It is also close to 9.5 mg L −1 in boreal catchments in Sweden (Winterdahl et al., 2014), and 8.1 mg L −1 measured in boreal wet and 10.5 mg L −1 in boreal dry catchments in Norway and Finland (de Wit et al., 2016). These differences suggest that climate, vegetation, and landscape heterogeneity may together shape how much DOC can be produced as well as when, where, and to what degree a hill is connected to a stream and the export of DOC at different times.

Temporal asynchrony of DOC production and export
The contrasting temporal patterns of simulated DOC production and export reflect the asynchronous nature of the two processes at the catchment scale. Local DOC production is influenced by the seasonal pattern of soil T , whereas the export is predominantly controlled by precipitation events and antecedent conditions that modulate the degree to which DOC production zones are hydrologically connected to the stream. This differs from studies showing the synchronization of DOC production and export in temperate climatic regions at soil pedons (Michalzik et al., 2001). This may be due to the relatively short water residence time and high connectivity at the pedon scale. The temporal asynchrony between DOC production and export is therefore strongly influenced by the seasonality of temperature and precipitation, which is shaped by local climate. At Shale Hills the wet winter and spring happens to be the cold season, whereas the dry summer is hot. In the summer, the catchment essentially produces and stores DOC in soil water and soil surfaces and waits for the arrival of the next storm to export. In other words, low hydrological connectivity in the summer imposes a lag period with respect to DOC export such that the DOC we see today is often the DOC produced a while ago. As such, the catchment acts as a DOC producer in the summer and a DOC exporter in spring and winter when the soil is wetter. These findings indicate a strong climatic control over DOC production and export. In places where climate seasonality is not as strong and catchments are hydrologically connected to streams throughout the year, we can expect to see DOC export all year long and, therefore, much less asynchrony. In places with strong seasonality, a few high-flow events can dominate the DOC export of the whole year. Under the Mediterranean climate with strong seasonality, for example, antecedent moisture conditions have been observed to be essential for understanding the temporal pattern of DOC and nutrient (N ) export (Bernal et al., 2005(Bernal et al., , 2002. Hydrological connectivity and water flow paths become dominant as subsurface saturation expands across valley floors and into hillslopes (Covino, 2017;Abbott et al., 2016).

Implications for vertical and lateral carbon fluxes
This work focuses on DOC lateral fluxes and does not simulate the carbon loss through soil respiration and associated vertical CO 2 fluxes, which has been the focus of some previous work Hasenmueller et al., 2015). Soil respiration is an important pathway of carbon flux that, similar to DOC production, can be shaped by soil temperature and moisture. Generally, warm temperature and medium soil moisture provide optimal conditions for microbial respiration, leading to significant vertical losses of carbon during summer months Stielstra et al., 2015). In contrast, low temperature and high soil moisture can hinder aerobic respiration and associated carbon losses as CO 2 , effectively accumulating DOC until the arrival of large storms. This pattern is consistent with observations that the total CO 2 release and DOC production are positively correlated (Neff and Hooper, 2002). The dependence of DOC production and export on soil T and soil moisture might also hold true for soil respiration. Conversely, as part of the sorbed DOC may be respired by microbes into CO 2 , our model might overestimate the DOC accumulation in the catchment, especially in summer.
This work does not consider the transport of particulate organic carbon (POC) in soil water and stream water, although POC can play an important role in the carbon budget and biogeochemical cycles in some cases (Ludwig et al., 1996;Diem et al., 2013). In Shale Hills, DOC comprises a major fraction (between 70 % and 80 %) of the total organic carbon export (Jordan et al., 1997). Similar patterns have been reported for organic carbon export at the global scale (Alvarez-Cobelas et al., 2012). However, POC export can be significant in anthropogenically impacted areas (Correll et al., 2001;Mattsson et al., 2005) with significant disturbance . They follow a different temporal pattern from DOC, due to different sources, transport mechanisms, and sensitivity to hydrologic variations (Dhillon and Inamdar, 2014;Alvarez-Cobelas et al., 2012).

Regulation of C-Q patterns
During dry periods, stream water mostly reflects the carbonpoor groundwater. As precipitation wets the catchment, the valley floor area that is characterized by a high SOC and DOC concentration is connected to the stream (Figs. 5, 7), elevating the stream DOC. This increase in DOC concentration continues until the catchment becomes wetter and expands the connected zones to the whole valley and swales. Under these conditions, the influence of high DOC in the vicinity of the stream fades and the DOC concentration decreases. This occurs at a threshold connectivity of about ∼ 0.1 (≈ the valley width divided by the catchment width).
In other words, during wet periods when the whole catchment is hydrologically connected to the stream, the stream DOC reflects the "average" concentration across the catchment (∼ 2.5 mg L −1 ). The increase and subsequent decrease pattern (or chevron pattern) therefore indicates the presence of three end-members from different sources: the groundwater with a very low DOC concentration, the soil water at stream beds and in organic-rich swales with the highest DOC content, and the uphill soil water with a DOC level in between these two. The overall dilution (or chevron) C-Q pattern observed here with a maximum at a mid-range discharge contrasts the commonly observed flushing pattern for DOC (Moatar et al., 2017). In fact, it more closely resembles the hysteresis behavior often observed in storm and snowmelt events for metals and nutrients (Zhi et al., 2019;Duncan et al., 2017). Previous field studies have illustrated that the hydrological connectivity to the stream versus the distribution of SOC ultimately dictates the spatial and temporal dynamics of the DOC concentration in soil and stream water, leading to different C-Q relationships (dilution versus flushing) (Bernhardt et al., 2017;Bernal and Sabater, 2012;Covino, 2017). This has been illustrated by different C-Q relationships at Shale Hills (USA) and Plynlimon (UK) (Herndon et al., 2015). Stream water at Shale Hills is derived from SOC-rich swales with a high DOC concentration at low flow and from both swales and hillslopes with a low DOC concentration when discharge increases. Conversely, at Plynlimon, SOC is enriched in uplands; therefore, concentrations are high at high flow when water flows connect SOC-rich uplands. Our reactive transport modeling provides a quantitative and mechanistic approach to explain the overall C-Q patterns, which have generally been interpreted as a production/source limitation (Covino, 2017;Zarnetske et al., 2018). Our results suggest that the spatial distribution of source zones and the degree of their connectivity to the stream determine when they are flushed out. Modeling approaches such as the one presented here can help us understand the mechanisms underlying C-Q patterns, and, thus, improve our ability to predict the evolution of C-Q trajectories under changing conditions. C-Q patterns also relate to the mixture of different sources of water in the stream, which is composed of the time-varying relative contribution from the shallow soil water and relatively deep groundwater. Their DOC contribution can be affected by the vertical distribution of reacting materials (Musolff et al., 2017;Bishop et al., 2004;Seibert et al., 2009;Winterdahl et al., 2016) and the relative volume contribution of source water (soil water versus groundwater below the soil-weathered rock interface) to the stream (Zhi et al., 2019;Radke et al., 2019;Weigand et al., 2017). Within the shale bedrock, the groundwater contribution to the stream is relatively small (∼ 7.5 %) at Shale Hills. Soil water (although from a very limited swale area) dominates inflow to the stream even during the summer dry period. When the groundwater volume input increases to about 18.8 % of the streamflow by volume (2.5 times the actual case; Fig. 12), the C-Q relationships shift to an overall flushing pattern. This may provide a potential explanation for the DOC C-Q flushing pattern at sandstone-dominant Garner Run (a neighboring catchment of Shale Hills), where the groundwater contributions to the stream are typically higher (Hoagland et al., 2017;Li et al., 2018). More interestingly, when the groundwater contribution is "sufficiently" high, it can mask the signature of the swale-derived soil water such that the three-end-member chevron C-Q pattern become a two-endmember pattern with monotonic flushing pattern that is similar to the observation in Coal Creek where groundwater contributes about 20 % annually (Zhi et al., 2019). C-Q relationships have been categorized into nine patterns, including three monotonic and six segmented types (Moatar et al., 2017;Underwood et al., 2017). The shifting threshold that separates segments of C-Q responses by the relative groundwater contribution in this work (Fig. 12) suggests that the relative contribution of groundwater to streamflow may play a pivotal role in shaping the C-Q patterns. This threshold value can potentially provide a rough estimation for the relative contribution of different end-members to the stream.
The mechanisms that regulate DOC C-Q patterns -seasonally variable hydrological connectivity and groundwater contribution -are consistent with previous literature on geogenic species (Mn, Fe), isotopes, and particle fluxes at Shale Hills Kim et al., 2018;Sullivan et al., 2016;Thomas et al., 2013). For example, Mn is associated with DOC via biotic cycling and storage in plant species, and Fe is associated with DOC via aqueous complexation. Therefore, both solutes are more abundant in shallow soils. The C-Q pattern of Fe and Mn shows a dilution pattern with concentrations decreasing as discharge increases (Herndon et al., 2015;. In the dry summer, stream water derives from rich-organic swales and riparian zones with high concentrations of soluble Fe and Mn , leading to corresponding high stream concentrations. At high flows, these solutes are diluted by the influx of uphill soil water without as much DOC. This again emphasizes the key role of solute sources and hydrological dynamics in controlling stream chemistry.

Conclusions
The production and export of DOC remain central uncertainties in determining the ecosystem-level carbon balance. These uncertainties persist because of complex interacting controls on DOC production and export. Indeed, few studies have quantitatively addressed the linkages between SOC processing, hydrological conditions, and corresponding DOC processing and export at the catchment scale. We found that DOC production was the major DOC source at Shale Hills (75 % compared with 23 % from precipitation and 2 % from groundwater). The simulations showed that the temporal dynamics of DOC export rates (R e ) were more linked to hydrological flow paths and precipitation events. A sensitivity analysis further confirmed that the DOC production rates R p were primarily controlled by temperature: changing the temperature dependence altered DOC concentrations significantly, whereas the effects of changing soil moisture dependence were relatively small. Conversely, DOC export was most sensitive to changes in hydrology, rendering more than 2 orders of magnitude differences in dry and wet periods. This difference in environmental drivers led to an asynchrony between DOC production and export. During the wet period (spring and winter), the catchment was well connected and DOC production and export occurred simultaneously. During summer, DOC accumulated (often in sorbed form) in soils disconnected from the stream, and DOC export was limited and constrained to the near-stream areas. In other words, the catchment serves as a DOC producer in the dry and hot summer but as an exporter in the wet and cold winter.
This work quantitatively demonstrates the key role of hydrological flow paths and the degree of connectivity in determining the C-Q patterns exhibited at the catchment outlet. At low discharges where connectivity is limited (I cs /Width < 0.1), stream DOC was mainly sourced from groundwater or from the valley floor and swales with enriched SOC. At higher discharges, an increasing contribution of soil lateral flow from planar hillslopes and uplands with low soil water DOC decreased the stream DOC concentration, ultimately rendering a dilution C-Q pattern. Although changing DOC reaction characteristics alters the soil water DOC concentration, there is little change in the overall C-Q patterns. However, when groundwater contributes 18.8 % of total annual discharge, the stream DOC concentration increases with discharge and flushing patterns emerge. This underscores the significance of the relative contribution and chemical signature of different water sources in shaping DOC export patterns. This study provides new insights into how DOC production and export interact at multiple scales and emphasizes the importance of considering different constraints when projecting the response of lateral and vertical carbon fluxes to climate changes. Data availability. The field data have been digitized and are accessible from the national CZO data portal (http://criticalzone.org/ shale-hills/data/datasets/; Brantley and Duffy, 2012). The source code of BFP (BioRT-Flux-PIHM) and the input files necessary to reproduce the results are available from the authors upon request (lili@engr.psu.edu).
Author contributions. HW, LL, and the co-authors conceived the idea and designed the numerical experiments based on ideas generated from a workshop and monthly discussions. HW ran the simulations, analyzed the simulation results, and wrote the first draft of the paper. All co-authors participated in editing the paper.