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

Lateral carbon flux 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 field 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 State Integrated Hydrologic Model, BFP), which met the satisfactory standard of a Nash–Sutcliffe efficiency (NSE) value greater than 0.5. We used the calibrated model to estimate and compare the daily DOC production rates (R_p; the sum of the local DOC production rates in individual grid cells) and export rate (R_e; the product of the concentration and discharge at the stream outlet, or load).

Results showed that daily R_p varied by less than an order of magnitude, primarily depending on seasonal temperature. In contrast, daily R_e varied by more than 3 orders of magnitude and was strongly associated with variation in discharge and hydrological connectivity. In summer, high temperature and evapotranspiration dried and disconnected hillslopes from the stream, driving R_p to its maximum but R_e to its minimum. During this period, the stream only exported DOC from the organic-poor groundwater and from organic-rich soil water in the swales bordering the stream. The DOC produced accumulated in hillslopes and was later flushed out during the wet and cold period (winter and spring) when R_e peaked as the stream reconnected with uphill and R_p reached its minimum.

The model reproduced the observed concentration–discharge (C–Q) relationship characterized by an unusual flushing–dilution pattern with maximum concentrations at intermediate discharge, indicating three end-members of source waters. A sensitivity analysis indicated that this nonlinearity was caused by shifts in the relative contribution of different source waters to the stream under different flow conditions. At low discharge, stream water reflected the chemistry of organic-poor groundwater; at intermediate discharge, stream water was dominated by the organic-rich soil water from swales; at high discharge, the stream reflected uphill soil water with an intermediate DOC concentration. This pattern persisted regardless of the DOC production rate as long as the contribution of deeper groundwater flow remained low (<18 % of the streamflow). When groundwater flow increased above 18 %, comparable amounts of groundwater and swale soil water mixed in the stream and masked the high DOC concentration from swales. In that case, the C–Q patterns switched to a flushing-only pattern with increasing 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 flow 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.

1 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₂ and CH₄ (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.

Figure 1 Attributes of the Susquehanna Shale Hills Critical Zone Observatory (SSHCZO): (a) surface elevation, (b) soil depth, and (c) soil organic carbon (SOC). 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 Tsites (squares).

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; Wen and Li, 2018; Wen et al., 2018). 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 (Tang et al., 2018) models have investigated the importance 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²) 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 (Brantley et al., 2018). 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 high-frequency measurements of soil properties, hydrology, and biogeochemistry (Brantley et al., 2018). 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. Although soil respiration is an important process, this study focuses on the net production and export of DOC.

2 Methods

2.1 Study site: a small catchment with an intermittent stream

The Shale Hills catchment is a 0.08 km², V-shaped, first-order 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).

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.

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), mid-slope (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 $\sim \mathrm{0.10},$ 0.20, 0.40, 0.70, 0.90, 1.00, and 1.30 m (Lin and Zhou, 2008).

2.2 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 (Duffy et al., 2014). 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 variation 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:

$$\begin{array}{}\text{(1)}& \begin{array}{rl}& V_i{\displaystyle \frac{\mathrm{d}\left(S_{\mathrm{w},i}{\mathit{\theta }}_i{C}_{m,i}\right)}{\mathrm{d}t}}=\\ & \sum _{j=N_{i,\mathrm{1}}}^{N_{i,x}}\left(A_{ij}{D}_{ij}{\displaystyle \frac{C_{m,j}-C_{m,i}}{l_{ij}}}-q_{ij}{C}_{m,j}\right)+r_{m,i},\\ & m=\mathrm{1},\mathrm{np},\end{array}\end{array}$$

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³) of each grid block; S_w is the soil moisture (m³ water m⁻³ pore volume); θ is porosity; C is the aqueous species concentration (mol m⁻³ water); t is time (s); N is the index of elements sharing surfaces; A is the grid interface area (m²); D is the diffusion/dispersion coefficient (m² s⁻¹); l is the distance (m) between the center of two neighboring grid blocks; q is the flow rate (m³ s⁻¹); r_m is the kinetically controlled reaction rates (mol s⁻¹) 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.

2.3 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₂ 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:

$$\begin{array}{}\text{(2)}& r_{\mathrm{p}}=kAf\left(T\right)f\left(S_{\mathrm{w}}\right),\end{array}$$

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⁻¹⁰ mol m⁻² s⁻¹) (Zhi et al., 2019; Wieder et al., 2014); and A is a lumped “surface area” (m², ($\mathrm{2.5}\times {\mathrm{10}}^{-\mathrm{3}}$ m² g⁻¹) × (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₁₀-based formation: $f\left(T\right)=Q_{\mathrm{10}}^{\left|T-\mathrm{10}\right|/\mathrm{10}}$, where Q₁₀ quantifies the rate increases with T, with the superscript 10 referring to a T value of 10 ^∘C (Davidson and Janssens, 2006). Q₁₀ 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)ⁿ 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 $C_{\mathrm{d}}\left(z\right)=C_{\mathrm{0}}\exp \left(-\frac{z}{b_m}\right)$, where C_d is SOC at depth z below the surface; C₀ 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: $\equiv X+\mathrm{DOC}\leftrightarrow \equiv X\mathrm{DOC}$, where ≡X and ≡XDOC represent the functional group without and with sorbed DOC, respectively (Rasmussen et al., 2018). 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_{\mathrm{eq}}=\frac{\left[\equiv X\mathrm{DOC}\right]}{\left[\equiv X\right]\left[\mathrm{DOC}\right]}$. 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 $\mathrm{4.0}\times {\mathrm{10}}^{-\mathrm{5}}$ to $\mathrm{6.0}\times {\mathrm{10}}^{-\mathrm{5}}$ mol g⁻¹ soil was used for Shale Hills (Jin et al., 2010; Li et al., 2017a) depending on the mineralogy in different zones of the catchment.

2.4 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⁻¹, 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⁻¹, which was the average concentration from the six field sampling locations in Fig. 1.

2.5 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⁻³–10⁰ m² g⁻¹ (Zhi et al., 2019; Chiou et al., 1990; Kaiser and Guggenberger, 2003). We also set K_eq using a literature range of 10⁰–10¹ (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.

Figure 3 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.

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_{\mathrm{cs}}/\mathrm{Width}=\mathrm{0.10}$, 0.35, and 0.70 in Fig. S3). Depending on the catchment geometry 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_{\mathrm{cs}}={\int }_{\mathrm{0}}^{\mathrm{\infty }}\mathit{\tau }\left(h\right)\mathrm{d}h$ 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.

2.6.2 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⁻¹) is the precipitation rate (m d⁻¹) times the rainfall DOC concentration ($\mathrm{6.0}\times {\mathrm{10}}^{-\mathrm{4}}$ mg m${}^{-\mathrm{3}}=$ 0.6 mg L${}^{-\mathrm{1}}\times {\mathrm{10}}^{-\mathrm{3}}$ L m⁻³) and the catchment drainage area (m²). The DOC input from groundwater R_g (mg d⁻¹) is the total groundwater influx (flow rate) times the groundwater DOC concentration (1.2 mg L⁻¹).

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 $\frac{{\mathrm{CV}}_{\left[\mathrm{DOC}\right]}}{{\mathrm{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 $\frac{{\mathrm{CV}}_{\left[\mathrm{DOC}\right]}}{{\mathrm{CV}}_Q}\ll \mathrm{1}$, C–Q patterns were considered chemostatic; values of $\left|b\right|>\mathrm{0.2}$ or $\frac{{\mathrm{CV}}_{\left[\mathrm{DOC}\right]}}{{\mathrm{CV}}_Q}\ge \mathrm{1}$, indicated a chemodynamic behavior. In the chemodynamic category, values of b>0.2 indicate flushing, whereas values of $b<-\mathrm{0.2}$ indicate dilution. We used the MATLAB curve-fitting toolbox to obtain the best fit model parameters.

2.7 Sensitivity analysis

We used a sensitivity analysis to explore the influence of soil T and moisture in the DOC production kinetics. The Q₁₀ in $f\left(T\right)=Q_{\mathrm{10}}^{\left|T-\mathrm{10}\right|/\mathrm{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\left(T\right)={\mathrm{4}}^{\left|T-\mathrm{10}\right|/\mathrm{10}}$. The rate dependence on soil moisture was explored using the base case f₁(S_w)=(S_w)ⁿ (increase behavior), and three additional functions (f₂, f₃, and f₄) 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

$$\begin{array}{}\text{(3)}& f_{\mathrm{2}}\left(S_{\mathrm{w}}\right)={\left({\displaystyle \frac{\mathrm{1}-S_{\mathrm{w}}}{\mathrm{0.6}}}\right)}^{\mathrm{0.77}},\end{array}$$

the constant behavior function was

$$\begin{array}{}\text{(4)}& f_{\mathrm{3}}\left(S_{\mathrm{w}}\right)=\mathrm{0.65},\end{array}$$

and the threshold behavior function was

$$\begin{array}{}\text{(5)}& f_{\mathrm{4}}\left(S_{\mathrm{w}}\right)=\left\{\begin{array}{ll}{\left(\frac{S_{\mathrm{w}}}{\mathrm{0.7}}\right)}^{\mathrm{1.5}}& S_{\mathrm{w}}\le \mathrm{0.7}\\ {\left(\frac{\mathrm{1}-S_{\mathrm{w}}}{\mathrm{1}-\mathrm{0.7}}\right)}^{\mathrm{1.5}}& S_{\mathrm{w}}>\mathrm{0.7}.\end{array}\right.\end{array}$$

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_{\mathrm{G}}=\mathrm{3.3}\times {\mathrm{10}}^{-\mathrm{4}}$ and $\mathrm{1.0}\times {\mathrm{10}}^{-\mathrm{4}}$ m d⁻¹ 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⁻¹). 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⁻¹.

Figure 4 (a) Temporal dynamics of measured and simulated stream DOC concentrations as well as groundwater and soil water DOC. The stream DOC (bright blue line) was from the soil water (light blue line) and groundwater Q_G (dark blue line). Under low-discharge conditions (e.g., July–September), Q_G contributed a larger proportion of discharge and stream DOC was more similar to groundwater DOC. Under wet conditions, stream DOC resembled soil water DOC from Q_L. (b–g) The local soil water DOC concentration for the six sampling locations shown in 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.

3 Results

3.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 ($\sim {\mathrm{10}}^{-\mathrm{2}}$ m d⁻¹) in spring and fall compared with summer with high soil T and high ET ($\sim {\mathrm{10}}^{-\mathrm{5}}$ m d⁻¹). The model captured the temporal dynamics of daily discharge, ET, and soil T with NSE values 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 $\mathrm{1.3}\times {\mathrm{10}}^{-\mathrm{4}}$ and $\mathrm{4.0}\times {\mathrm{10}}^{-\mathrm{5}}$ m d⁻¹ 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 fluctuations. High summer ET drove the catchment to drier conditions, thereby decreasing the connectivity to the stream.

3.2 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<\mathrm{1.0}\times {\mathrm{10}}^{-\mathrm{4}}$ m d⁻¹), 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.

Figure 5 Spatial profiles in May (wet), August (dry), and October (wet after dry) of 2009: (a) soil T, (b) soil moisture, (c) hydrologically connected zones, (d) local DOC production rates r_p, and (e) soil water DOC concentration. 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.

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., a NSE value of 0.36 (SPRT), 0.42 (SPMS), 0.60 (SPVF), 0.46 (SSRT), 0.40 (SSMS), and 0.51 (SSVF). The variation in model performance at different locations may arise from the lack of detailed information on local soil properties and 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.

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).

Figure 7 (a) Relationship of daily discharge (Q) with stream DOC concentration: open circles are simulations and filled circles with a black outline are data. (b) 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.

3.3 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-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⁵ mg d⁻¹) was greater than the daily R_r from rainfall (1.6×10⁵ mg d⁻¹) or groundwater R_g (1.2×10⁴ mg d⁻¹). 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_{\mathrm{p}}+R_{\mathrm{r}}+R_{\mathrm{g}}$) or R_p. Instead, it primarily followed the discharge patterns: large rainfall events exported disproportionally 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⁶ mg over the year, which was about 1.0 % of the overall DOC production, indicating a general mass balance at the catchment scale.

Figure 9 Sensitivity analysis of temporal DOC rates for (a) soil temperature f(T) and (b) soil moisture f(S_w). A varying Q₁₀ 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.

3.4 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=-\mathrm{0.23}$ and $\frac{{\mathrm{CV}}_{\left[\mathrm{DOC}\right]}}{{\mathrm{CV}}_Q}=\mathrm{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 ($<\mathrm{1.8}\times {\mathrm{10}}^{-\mathrm{4}}$ m d⁻¹) with small water storage (0.25–0.28 m) and connectivity ($I_{\mathrm{cs}}/\mathrm{Width}<\mathrm{0.1}$) (Fig. 7b), the stream DOC was a mix of organic-poor groundwater and organic-rich 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 valley 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⁵–10⁶ mg d⁻¹) 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.

Figure 11 C–Q relationships under different (a) T, (b) f(S_w), 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⁻¹ 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.

3.5 Sensitivity analysis

3.5.1 Control of temperature, soil moisture, and sorption

Higher Q₁₀ values in f(T) led to more pronounced seasonality in R_p (Fig. 9a). The R_p for Q₁₀=4.0 was more than 4 times higher than that of Q₁₀=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₁₀ 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₁₀ 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₁₀ 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.

3.5.2 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<\mathrm{1.8}\times {\mathrm{10}}^{-\mathrm{4}}$ m d⁻¹). 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=\mathrm{1.8}\times {\mathrm{10}}^{-\mathrm{4}}$ to about $\mathrm{1.0}\times {\mathrm{10}}^{-\mathrm{3}}$ m d⁻¹. This reflects the relative 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 Total annual R_p (red) and R_e (blue) under two groundwater volume contribution conditions ($Q_{\mathrm{G}}/Q=\mathrm{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.

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₁₀ because a higher production led to higher DOC export. Annual R_p also depended on f(S_w), with the threshold function f₄(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 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.

4 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₂ 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₂ across terrestrial and freshwater ecosystems.

4.1 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₄(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 co-generate 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⁻¹), compared with 3.0 mg L⁻¹ in temperate humid catchments in Germany (Musolff et al., 2018), 4.1 mg L⁻¹ in some UK catchments with oceanic climate (Monteith et al., 2015), and 4.5 mg L⁻¹ in France (Moatar et al., 2017). It is also close to 9.5 mg L⁻¹ in boreal catchments in Sweden (Winterdahl et al., 2014), and 8.1 mg L⁻¹ measured in boreal wet and 10.5 mg L⁻¹ 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.

4.2 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, 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).

4.3 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₂ fluxes, which has been the focus of some previous work (Brantley et al., 2018; 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 (Perdrial et al., 2018; Stielstra et al., 2015). In contrast, low temperature and high soil moisture can hinder aerobic respiration and associated carbon losses as CO₂ (Smith et al., 2018), effectively accumulating DOC until the arrival of large storms. This pattern is consistent with observations that the total CO₂ 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₂, 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 (Abbott et al., 2016). 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).

4.4 Regulation of C–Q patterns

During dry periods, stream water mostly reflects the carbon-poor 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⁻¹). 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-end-member 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 (Herndon et al., 2018; 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; 2018). In the dry summer, stream water derives from rich-organic swales and riparian zones with high concentrations of soluble Fe and Mn (Herndon et al., 2018), 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.

5 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.