BioRT-Flux-PIHM v1.0: a biogeochemical reactive transport model at the watershed scale

Watersheds are the fundamental Earth surface functioning units that connect the land to aquatic systems. Many watershed-scale models represent hydrological processes but not biogeochemical reactive transport processes. This has limited our capability to understand and predict solute export, water chemistry and quality, and Earth system response to changing climate and anthropogenic conditions. Here we present a recently developed BioRT-Flux-PIHM (BioRT hereafter) v1.0, a watershed-scale biogeochemical reactive transport model. The model augments the previously developed RT-Flux-PIHM that integrates land-surface interactions, surface hydrology, and abiotic geochemical reactions. It enables the simulation of (1) shallow and deep-water partitioning to represent surface runoff, shallow soil water, and deeper groundwater and of (2) biotic processes including plant uptake, soil respiration, and nutrient transformation. The reactive transport part of the code has been verified against the widely used reactive transport code CrunchTope. BioRT-Flux-PIHM v1.0 has recently been applied in multiple watersheds under diverse climate, vegetation, and geological conditions. This paper briefly introduces the governing equations and model structure with a focus on new aspects of the model. It also showcases one hydrology example that simulates shallow and deep-water interactions and two biogeochemical examples relevant to nitrate and dissolved organic carbon (DOC). These examples are illustrated in two simulation modes of complexity. One is the spatially lumped mode (i.e., two land cells connected by one river segment) that focuses on processes and average behavior of a watershed. Another is the spatially distributed mode (i.e., hundreds of cells) that includes details of topography, land cover, and soil properties. Whereas the spatially lumped mode represents averaged properties and processes and temporal variations, the spatially distributed mode can be used to understand the impacts of spatial structure and identify hot spots of biogeochemical reactions. The model can be used to mechanistically understand coupled hydrological and biogeochemical processes under gradients of climate, vegetation, geology, and land use conditions.

1 Introduction

Watersheds are the fundamental Earth surface units that receive and process water, mass, and energy (Li, 2019; Li et al., 2021; Hubbard et al., 2018). Watershed processes include land-surface interactions that regulate evapotranspiration and discharge as well as water partitioning between shallow soil lateral flow going into streams vs. downward flow and recharge into the deeper subsurface (Fig. 1). Complex biogeochemical interactions occur among soil, water, roots, and microbes along water's flow paths, regulating gas effluxes (e.g., CO₂) and solute export (Fatichi et al., 2019; van der Velde et al., 2010; Benettin et al., 2020; Adler et al., 2021).

These hydrological and biogeochemical processes determine how Earth surface systems respond to hydroclimatic forcing and human perturbations. Understanding these processes remains challenging due to the complex process interactions. An example is the concentration–discharge (C–Q) relationships of solutes at stream and river outlets. Similar C–Q relationships have been observed for some solutes across watersheds under distinct geological and climatic conditions (Godsey et al., 2009; Basu et al., 2010; Moatar et al., 2017; Zarnetske et al., 2018; Godsey et al., 2019), whereas different solutes have shown contrasting patterns in the same watershed (Miller et al., 2017; Herndon et al., 2015; Musolff et al., 2015; Stewart et al., 2021; Zhi et al., 2020b). A general theory that can explain contrasting C–Q observations (e.g., flushing vs. dilution behaviors) under diverse watershed characteristics and external conditions remains elusive. The lack of mechanistic understanding presents major roadblocks to forecasting water quality and Earth system dynamics in the future.

One of the outstanding challenges is the lack of modeling tools that mechanistically link hydrological and biogeochemical processes at the watershed scale. Model development has advanced mostly within the disciplinary boundaries of hydrology and biogeochemistry (Li, 2019). Watershed hydrologic models focus on solving for water storage and fluxes (Fatichi et al., 2016). Reactive transport models (RTMs) have traditionally focused on transport and multi-component biogeochemical reactions mostly in groundwater systems with limited interactions with climate and other surficial conditions (Steefel et al., 2015; Li et al., 2017b). Some integration crossing disciplinary boundaries did occur in recent years. For example, SWAT (Soil and Water Assessment Tool) has a version that couples with the groundwater model MODFLOW and the surface water and groundwater quality model in RT3D (Bailey et al., 2017; Ochoa et al., 2020). CATHY (Catchment Hydrology) includes processes of pesticide decay (Gatel et al., 2019; Scudeler et al., 2016). Hydrologiska Byråns Vattenbalansavdelning (HBV) and the Hydrological Predictions for the Environment (HYPE) have modules that simulate processes relevant to nutrients and contaminants (Lindström et al., 2005, 2010). While these models can simulate processes such as leaching of nutrients from agriculture lands (Lindström et al., 2005, 2010; Bailey et al., 2017), they do not explicitly solve the multi-component reactive transport equations. In other words, reactions are often represented rudimentarily without honoring kinetics and thermodynamics theories in soil biogeochemistry and geochemistry. For example, nutrient leaching is calculated based on empirical equations without explicitly solving reactive transport equations. Reaction rates are represented using first-order decay (Gatel et al., 2019), assuming constant reaction rates that do not change with environmental conditions. Biogeochemical processes however are highly variable with seasonal dynamics and depend on local environments such as substrate availability, soil temperature, and soil moisture (Li et al., 2017a; Davidson and Janssens, 2006). These models therefore cannot capture the temporal variations in environmental factors that regulate soil biogeochemical reactions and stream and water chemistry.

To fill this model capability need, we augmented the watershed model RT-Flux-PIHM (Bao et al., 2017) into BioRT-Flux-PIHM (BioRT hereafter). Compared to RT-Flux-PIHM, BioRT has two additions. One is the capability of simulating biotic processes including plant uptake of nutrients, soil respiration, and other microbe-mediated redox reactions. Examples include soil respiration that produces CO₂ and dissolved organic carbon (DOC) and nutrient transformation reactions such as nitrification and denitrification. The other is the addition of an optional deeper layer below shallow soil to enable the simulation of interacting deep-water and shallow soil water flow (Fig. 1). Here the deep water is loosely defined as the water below the soil zone, typically in less weathered, fractured subsurface that harbors relatively old and slow-moving groundwater contributing to streams. It is a fundamental component of the hydrologic cycle and water budget. The groundwater–surface-water interactions also modulate land–atmospheric energy exchanges and soil moisture dynamics (Keune et al., 2016). Evidence has been mounting in recent years that deeper water below the shallow soil interacts with streams, introduces water with distinct chemistry, sustains base flow in dry times, and buffers climate variability (Gurdak, 2017; Green, 2016; Taylor et al., 2013). Stream chemistry often reflects distinct chemistry from shallow soil water and deeper groundwater at different times, i.e., the so-called Shallow and Deep Hypothesis (Zhi et al., 2019; Zhi and Li, 2020; Botter et al., 2020). Including the deep-water component thus is essential for understanding mechanisms and predicting dynamics of water quality under changing climate and human conditions.

Figure 1 A conceptual diagram for processes at the watershed scale. Land-surface interactions include processes such as energy balance, solar radiation, and evapotranspiration; hydrological processes partition water between surface runoff, shallow soil water, and deeper water. Soil biogeochemical reactions include abiotic reactions such as weathering (e.g., mineral dissolution and precipitation, ion exchange, surface complexations) and biotic processes such as plant uptake and leaching of nutrients, soil respiration, and other microbe-mediated reactions. These processes are represented in three modules: the Flux module for land-surface interactions, the PIHM module for catchment hydrology, and the BioRT module for biogeochemical reactions. Conceptually the shallow zone is the shallow soil and weathered zone that is more conductive to water flow (e.g., soil lateral flow or interflow). The deep zone refers to the less weathered zone that often harbors the old and slow-flowing groundwater. Reactions can occur in both shallow and deep zones.

This paper introduces the new developments in BioRT. The paper starts with a brief overview of water- and energy-related processes. It then introduces governing equations and reaction kinetics used in BioRT, followed by three examples that illustrate the new capabilities. The examples include the surface water and groundwater interactions, nitrate transformation and transport, and the production and export of DOC. The model can be set up in both spatially lumped or spatially explicit modes. The source code and the examples here are archived on the Zenodo website (https://doi.org/10.5281/zenodo.3936073) and the GitHub website (https://github.com/Li-Reactive-Water-Group/BioRT-Flux-PIHM, last access: 23 September 2021).

2 Model overview

BioRT-Flux-PIHM integrates three modules (Fig. 1). The Flux module is for land-surface processes including surface energy balance, solar radiation, and evapotranspiration (ET) (Shi et al., 2013). The hydrology module PIHM simulates water processes including precipitation, interception, infiltration, recharge, surface runoff, subsurface lateral flow, and deep-water flow (Qu and Duffy, 2007). The BioRT module simulates solute transport and reactions. The abiotic reactions included in RT-Flux-PIHM (Bao et al., 2017) are mineral dissolution and precipitation, aqueous and surface complexation, and ion exchange reactions. The newly added reactions include plant uptake of nutrients, soil respiration, and microbe-mediated redox reactions (e.g., carbon decomposition and nutrient transformation).

The land-surface and hydrology modules solve for soil temperature and water storage, from which water fluxes are calculated for surface runoff and shallow- and deep-water fluxes. The BioRT module uses the calculated soil temperature, water storage, and water fluxes to simulate solute transport (advective and diffusive–dispersive) and biogeochemical reactions in both shallow and deep zones (see governing equations in later sections). The reactions include kinetically controlled (e.g., microbe-mediated redox reaction, mineral dissolution, and precipitation) or equilibrium-controlled ones (e.g., ion exchange, surface complexation (sorption), and aqueous complexation). Users can define the types of reactions to be included and the form of reaction kinetics in input files. The output of BioRT includes the spatial distribution and time series of aqueous and solid concentrations, from which we can infer reaction rates.

The model can be set up running in either spatially lumped or spatially explicit modes. When running in spatially explicit mode, the simulation domain can be structured as prismatic grids based on topography. Each grid is partitioned into surface and shallow and deep subsurface layers. The surface layer calculates water flow above ground (surface runoff). The shallow zone is loosely defined as the highly permeable subsurface that contrasts the deep zone that is broadly defined as the lower permeability zone beyond the shallow zone. In many places, the shallow zone is the soil zone that is most conductive to water flow (e.g., lateral flow) and is responsive to hydroclimatic forcing. The deep subsurface zone is the less weathered layer that harbors the old groundwater that contributes to stream flow. Note that these definitions differ from those in the hydrology community, which often refers to the shallow soil lateral flow as groundwater, in a way that distinguishes it from the surface runoff (Dingman, 2015). These source waters from different depths of the subsurface often have distinct solid and water chemistry and are dominant at different hydrological conditions at different times of the year, as has been observed in many catchments and watersheds (Brantley et al., 2018; Zhi et al., 2019; Zhi and Li, 2020; Sullivan et al., 2016). The model is flexible for taking inputs from online data portals or local measurements, and it can accommodate low data availability (see the following Sect. 5 for data need and domain setup).

Figure 2 Hillslope view of the shallow and deep zones and relevant water flows. Streams received water primarily from three water flows: the surface runoff (q_overland), lateral flow from the shallow zone $\left(q_{\text{lateral}}^{\text{sl}}\right)$, and the lateral flow that has been recharged and eventually come out from the deeper zone $\left(q_{\text{lateral}}^{\text{dp}}\right)$. The symbols “h”, “e”, and “q” denote water storage, evapotranspiration, and water flow, respectively. The superscript letters “sl” and “dp” refer to shallow and deep zones, respectively. The subscript letters “u” and “s” refer to unsaturated and saturated layers, respectively. Detailed equations are in Sect. 3.1 and 3.2. The terms “infil”, “u2s”, and “recharge” refer to infiltration, unsaturated to saturated zones, and recharge.

3 Governing equations and processes

3.1 Water equations

Flux-PIHM simulates surface runoff and a lumped subsurface flux into streams without distinguishing shallow soil water and deeper groundwater flow paths. Mounting evidence has shown that the shallow soil water and deeper groundwater have distinct chemistry and are dominant at different times of the year (Xiao et al., 2021; Zhi and Li, 2020; Zhi et al., 2019). This means that a lumped subsurface flow cannot describe the dynamics of stream chemistry. We therefore added a deeper groundwater zone to simulate deeper water flows that interact with streams (Fig. 2). Each prismatic element now has three zones in the vertical direction: surface (or above ground) and shallow and deep zones in the subsurface.

In each prismatic element i, the shallow zone includes unsaturated and saturated water storage. The unsaturated zone receives water from the surface via infiltration and flows vertically to the saturated zone. The saturated zone flows both vertically to the deep zone (recharge) and laterally to neighboring grids j or the stream (lateral). The code solves the following equations in the shallow zone:

$$\begin{array}{}\text{(1)}& {\displaystyle }{\mathit{\theta }}_i^{\text{sl}}{\displaystyle \frac{\mathrm{d}{h}_{i,\mathrm{u}}^{\text{sl}}}{\mathrm{d}t}}=q_{i,\text{inf}}-q_{i,\text{u2s}}^{\text{sl}}-e_{i,\mathrm{u}}^{\text{sl}},\text{(2)}& {\displaystyle }{\mathit{\theta }}_i^{\text{sl}}{\displaystyle \frac{\mathrm{d}{h}_{i,\mathrm{s}}^{\text{sl}}}{\mathrm{d}t}}=q_{i,\text{u2s}}^{\text{sl}}-q_{i,\text{rechg}}-e_{i,\mathrm{s}}^{\text{sl}}+\sum _{\mathrm{1}}^{N_{ij}}{q}_{ij}^{\text{sl}},\end{array}$$

where ${\mathit{\theta }}_i^{\text{sl}}$ [m³ pore space m⁻³ total volume] is the shallow-zone porosity in element i; $h_{i,\mathrm{u}}^{\text{sl}}$ and $h_{i,\mathrm{s}}^{\text{sl}}$ [m] are the unsaturated and saturated water storage in the shallow zone, respectively. The storages h here are the height of soil column with equivalent saturated water, not the height of the pure water (100 % volume) column. That is why porosity is in the equation. For saturation zones, this height is needed to quantify the depths of water tables, and it determines the direction of water flow between neighboring grids. The q_i,inf [m s⁻¹] is the infiltration rate from the surface to the shallow zone; $q_{i,\text{u2s}}^{\text{sl}}$ [m s⁻¹] is the vertical flow from the unsaturated layer to the saturated layer in the shallow zone; q_i,rechg [m s⁻¹] is the recharge rate from the shallow zone to the deep zone; $e_{i,\mathrm{u}}^{\text{sl}}$ and $e_{i,\mathrm{s}}^{\text{sl}}$ [m s⁻¹] are evapotranspiration (ET) from the unsaturated and saturated layer in the shallow zone, respectively; $q_{ij}^{\text{sl}}$ [m s⁻¹] is the lateral flows in the shallow saturated layer between the element i and its neighbor element j; N_ij (≤3) is the number of neighbor elements j. For a prismatic element i, a boundary cell has one or two neighbors; a non-boundary cell has three neighbors. ET is calculated by the Penman potential evaporation scheme (detailed equations in Shi, 2012). A similar set of water equations for the deep zone are found in the Supplement (Eqs. S1 and S2).

Figure 3 Biotic and abiotic reactions relevant to the transformation of soil organic matter (SOM). SOM can become stabilized (recalcitrant) through sorption on clay and separation from reactants. Labile OM can decompose into inorganic forms, releasing C, N, and P that further transform into various forms (taken from Li, 2019, with permission from the Mineralogical Society of America).

Infiltration and vertical fluxes from the unsaturated to saturated layer in the shallow zone are based on the Richards equation, in which hydraulic water head H (i.e., the summation of water storage h and elevation head z) and hydraulic conductivity K determine the fluxes:

$$\begin{array}{}\text{(3)}& {\displaystyle }{q}_{i,\text{inf}}=K_{i,\text{inf}}{\displaystyle \frac{H_{i,\text{sur}}-H_{i,\mathrm{u}}^{\text{sl}}}{d_{i,\text{inf}}}},\text{(4)}& {\displaystyle }{q}_{i,\text{u2s}}^{\text{sl}}=K_{i,\mathrm{V}}^{\text{sl}}{\displaystyle \frac{H_{i,\mathrm{u}}^{\text{sl}}-H_{i,\mathrm{s}}^{\text{sl}}}{\mathrm{0.5}{d}_i^{\text{sl}}}},\end{array}$$

where d_i,inf and $d_i^{\text{sl}}$ [m] are the thickness of infiltration layer and shallow-zone depth, respectively; K_i,inf [m s⁻¹] is the hydraulic conductivity of the infiltration layer, the top 0.1 m of the subsurface that has different conductivity from the rest of subsurface; $K_{i,\mathrm{V}}^{\text{sl}}$ [m s⁻¹] is the hydraulic conductivity in the vertical direction (i.e., weighted average of macropore K_i,macV and soil matrix K_i,satV, Eq. S7); H_i,sur [m] is the surface hydraulic water head ($=h_{i,\text{sur}}+z_{i,\text{sur}}$); $H_{i,\mathrm{u}}^{\text{sl}}$ and $H_{i,\mathrm{s}}^{\text{sl}}$ [m] are the shallow hydraulic water heads in the unsaturated and saturated layer, respectively. The lateral flow in the shallow saturated layer is calculated using Darcy's law:

$$\begin{array}{}\text{(5)}& q_{ij}^{\text{sl}}=K_{ij}^{\text{sl}}{\displaystyle \frac{H_{i,\mathrm{s}}^{\text{sl}}-H_{j,s}^{\text{sl}}}{d_{ij}}},\end{array}$$

where d_ij [m] is the distance between the centers of elements i and j; $K_{ij}^{\text{sl}}$ [m s⁻¹] is the harmonic mean of shallow hydraulic conductivity in the horizontal direction between elements i ($K_{i,\mathrm{H}}^{\text{sl}}$) and j ($K_{j,\mathrm{H}}^{\text{sl}}$). The interaction between the shallow saturated zone and stream channel also follows Eq. (5), except that the adjacent head is replaced by the level of stream water. Similar to the shallow zone, hydrological equations in the deep zone are detailed in the Supplement (Eqs. S1–S8).

3.2 Reactive transport equations

The governing advection dispersion reaction (ADR) equation for an arbitrary solute m in an element i is as follows (Bao et al., 2017):

$$\begin{array}{}\text{(6)}& \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 _{\mathrm{1}}^{N_{ij}}\left(A_{ij}{D}_{ij}{\displaystyle \frac{C_{m,j}-C_{m,i}}{d_{ij}}}-q_{ij}{A}_{ij}{C}_{m,j}\right)\\ & +R_{m,i},m=\mathrm{1},\mathrm{\dots },nm,\end{array}\end{array}$$

where V_i [m³ total volume] is the total volume of element i; S_w,i [m³ water m⁻³ pore space] is soil water saturation; θ_i [m³ pore space m⁻³ total volume] is the porosity; C_m,i [mol m⁻³ water] is the aqueous concentration of species m; N_ij is the number of fluxes from neighbor element j for element i, and N_ij is 2 for the unsaturated zone (infiltration, recharge) with only vertical flows and 5 for saturated zone with flux from (or to) the unsaturated zone, from (or to) the deeper zone, and fluxes between i and three neighbor elements j in lateral flow directions for non-boundary grids; A_ij [m²] is the grid area shared by i and its neighbor grid j; D_ij [m² s⁻¹] is the hydrodynamic dispersion coefficient (i.e., sum of mechanical dispersion and effective diffusion coefficient) normal to the shared surface A_ij; d_ij [m] is the distance between the center of i and its neighbor elements j; q_ij [m s⁻¹] is the flow rate across A_ij; R_m,i [mol s⁻¹] is the total rate of kinetically controlled reactions in element i that involve species m; nm is the total number of independent primary species to be solved for reactive transport equations. Equation (6) states that the change in solute mass (the left term) is driven by dispersive transport, advective transport, and reactions (i.e., the first, second, and third right-hand side terms, respectively).

3.3 Biogeochemical processes and reaction kinetics

3.3.1 Biogeochemical processes

Here we discuss representative biogeochemical processes that involve plants and microbes that can be included in BioRT. BioRT differs from general water quality models that primarily target a few contaminants (e.g., N, P, metals). The framework of BioRT is flexible, and the users can define reactions and solutes of interest in the input files. For abiotic reactions such as mineral dissolution and surface complexation or ion exchange, readers are referred to Bao et al. (2017). Generally speaking, shallow soils contain more weathered materials and soil organic matters (SOMs) including roots, leaves, and microbes (Fig. 3). SOM can be decomposed partially into organic molecules that dissolve in water, i.e., dissolved organic carbon (DOC). It can also become oxidized completely into CO₂, which can emit back to the atmosphere in gas form (Davidson, 2006) or transport and enter streams in the form of dissolved inorganic carbon (DIC). With coexisting cations (e.g., Ca, Mg), DIC can precipitate out and become carbonate minerals (e.g., CaCO₃).

Organic matter (OM) decomposition releases organic nitrogen (R-NH₂), which can further react to become ${\mathrm{NH}}_{\mathrm{4}}^{+}$ and other nitrogen forms (N₂, N₂O, NO, ${\mathrm{NO}}_{\mathrm{2}}^{-}$, NO₂) (Fig. 3). The gases can emit back to the atmosphere. Denitrification requires anoxic conditions and occurs less commonly in shallow soils owing to the pervasive presence of O₂ (Sebestyen et al., 2019); denitrification can become important under wet conditions and in O₂-depleted groundwater aquifers. Phosphorous (P) can be in organic forms in organic matter, sorbed on fine soil particles, dissolved in water, or in solid forms as P-containing minerals. Transformation of nutrients occurs through various bio-mediated or abiotic reactions. A representative P-containing mineral in the Earth's crust is apatite Ca₅(PO₄)₃ (F, Cl, OH). Once liberated via rock dissolution, P is biologically assimilated and locked in organic forms. These organic forms have very low solubility, allowing them to bind on and be transported together with soil particles in the form of orthophosphate or pyro-diphosphate.

3.3.2 Reaction kinetics in natural soils

Rate dependence on temperature and soil moisture

Reactions such as soil respiration and plant uptake typically depend on environmental conditions (temperature or soil moisture). For example, in shallow oxic soils where organic carbon and O₂ are often abundant, the rate law for carbon decomposition can be simplified to the following form assuming microorganism concentrations are relatively constant.

$$\begin{array}{}\text{(7)}& r=kAf\left(T\right)f\left(S_{\mathrm{w}}\right)f\left(Z_{\mathrm{w}}\right),\end{array}$$

where the reaction rate r [mol s⁻¹] depends on rate constant k [$\mathrm{mol}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{s}}^{-\mathrm{1}}$], the surface area A [m²] is a lumped parameter that quantitatively represents SOM content and biomass abundance, f(T) and f(S_w) describe the temperature and soil moisture dependence, respectively, f(Z_w) can be included to account for the depth distribution of SOM (Seibert et al., 2009), and Z_w [m] is the water table depth. An example for the depth distribution is $f\left(Z_{\mathrm{w}}\right)=\exp (-Z_{\mathrm{w}}/b_m)$ (Weiler and McDonnell, 2006), with b_m as the depth coefficient describing the gradient of SOM content over depth. Users can choose to include either one or all of these dependencies in input or database files.

Figure 4 Reaction rate dependence. (a) Function form of soil temperature dependence and (b, c) soil moisture dependence for reaction rates. The f(T) takes the Q₁₀ form (Eq. 8). The soil moisture factor f(S_w) depends on S_w,c and n and soil water saturation S_w (Eq. 9). The soil moisture function can represent three types of behaviors: the threshold behavior (b, $\mathrm{0}<S_{\mathrm{w},\mathrm{c}}<\mathrm{1}$), increase behavior (red in c, $S_{\mathrm{w},\mathrm{c}}=\mathrm{1}$), and decrease behavior (blue in c, $S_{\mathrm{w},\mathrm{c}}=\mathrm{0}$). Values of n=1 lead to a linear threshold dependence of S_w, while n<1 and n>1 lead to concave and convex dependences, respectively.

The temperature dependence follows a Q₁₀-based form (Lloyd and Taylor, 1994; Davidson and Janssens, 2006) as follows:

$$\begin{array}{}\text{(8)}& f\left(T\right)=Q_{\mathrm{10}}^{(T-\mathrm{20})/\mathrm{10}},\end{array}$$

where Q₁₀ is the relative increase in reaction rates when temperature increases by 10 ^∘C. Values of Q₁₀ (Fig. 4a) can vary from 1.0 to 3.0, depending on climatic conditions, substrate availability, and ecosystem type (e.g., grassland, forest) (Davidson et al., 2006; Liu et al., 2017). The mean values are in the range of 1.4 to 2.5 (Zhou et al., 2009; Bracho et al., 2016). The Q₁₀ value can be specified in the input file.

The soil moisture dependence function f(S_w) is coded in the following form:

$$\begin{array}{}\text{(9)}& \left\{\begin{array}{ll}{\left(\frac{S_{\mathrm{w}}}{S_{\mathrm{w},\mathrm{c}}}\right)}^n,& S_{\mathrm{w}}\le S_{\mathrm{w},\mathrm{c}}\\ {\left(\frac{\mathrm{1}-S_{\mathrm{w}}}{\mathrm{1}-S_{\mathrm{w},\mathrm{c}}}\right)}^n,& S_{\mathrm{w}}>S_{\mathrm{w},\mathrm{c}}.\end{array}\right.\end{array}$$

Here S_w,c [0 to 1] is the critical soil moisture (saturation) at which rates are highest, and n is the exponent reflecting the dependence of rates on soil moisture. A typical n value is 2 (Yan et al., 2018) with a range between 1.2 and 3.0 (Hamamoto et al., 2010), depending on soil structure and texture. As shown in Fig. 4b, the form indicates an intermediate critical soil moisture S_w,c at which f(S_w) reaches its maximum. When $S_{\mathrm{w}}\le S_{\mathrm{w},\mathrm{c}}$, f(S_w) increases with S_w; when $S_{\mathrm{w}}>S_{\mathrm{w},\mathrm{c}}$, f(S_w) decreases with S_w (Fig. 4b). Under the extreme conditions of S_w,c equaling 0 or 1, f(S_w) monotonically increases or decreases (Fig. 4c). The two parameters, S_w,c and n, determines the shape of the curve. They can be specified in input or database files. One can also choose not to have temperature or soil moisture dependence by choosing parameters that would lead to the value of the exponent being zero

Rate dependence on substrates: Monod kinetics and the biogeochemical redox ladder

Deeper groundwater aquifers often experience anoxic conditions that lead to processes such as denitrification or methanogenesis. This can also happen in wetlands or wet soils. Under such conditions, the rates of microbe-mediated redox reactions depend not only on temperature and soil moisture as discussed above, they also depend on concentrations of electron donors and non-oxygen electron acceptors (e.g., nitrate, iron oxides, sulfate) that are often limited under anoxic conditions (Bao et al., 2014; Li, 2019). The order of redox reactions typically follows the biogeochemical redox ladder, which is based on how much microbe can harvest energy by reducing different types of electron acceptors. Monod reaction rate laws are often used for quantifying rates of these redox conditions. These rate laws are detailed in Sect. S2 in the Supplement. Users can combine these Monod rate laws and the temperature and soil moisture dependence described above if needed.

3.4 Plant-related processes: root uptake of nitrate as an example

Nutrient uptake by plants is complex and remains poorly understood. A variety of plant uptake models exists with varying degrees of complexity (Neitsch et al., 2011; Fisher et al., 2010; Cai et al., 2016). These models are mostly based on a plant growth module or a supply and demand approach that often requires detailed phenological and plant attributes including growth cycle, root age and biomass, nutrient availability, and carbon allocation, in addition to local temperature and soil moisture (Dunbabin et al., 2002; Buysse et al., 1996; Fisher et al., 2010). Without a detailed mechanistic understanding, we assume a simple and operational approach. In Example 2, which we show later, for example, nitrate uptake was modeled with dependence on ${\mathrm{NO}}_{\mathrm{3}}^{-}$ concentration, soil temperature and moisture, and rooting density (McMurtrie et al., 2012; Yan et al., 2012; Buljovcic and Engels, 2001).

$$\begin{array}{}\text{(10)}& {\displaystyle }{r}_{\text{uptake}}=k_{\text{uptake}}{C}_{{\mathrm{NO}}_{\mathrm{3}}^{-}}f\left(T\right)f\left(S_{\mathrm{w}}\right)f_{\text{root}}\left(d_{\mathrm{w}}\right),\text{(11)}& {\displaystyle }{f}_{\text{root}}\left(d_{\mathrm{w}}\right)=\exp \left(\left(-d_{\mathrm{w}}+\mathit{\delta }/\mathit{\lambda }\right)\right),\end{array}$$

where k_uptake [L s⁻¹] is the nitrate uptake rate, f_root(d_w) is the normalized rooting density term in the range of 0 to 1 as a function of water depth to the groundwater (d_w). The rooting term (Eq. 14) was exponentially fitted (δ=0.013, λ=0.20) based on field measurements of root distribution along depth (Hasenmueller et al., 2017). It is common to observe a root density decrease exponentially in forests (López et al., 2001). Other forms of user-tailored plant uptake rate law can be added if needed.

4 Numerical scheme and model verification

The system of differential equations for water storage (e.g., Eqs. 1 and 2, and S1 and S2) is assembled into a global system of ordinary differential equations (ODEs). It is solved in CVODE (short for C-language Variable-coefficients ODE solver, https://computing.llnl.gov/projects/sundials/cvode, last access: 10 January 2021), a numerical ODE solver in the SUite of Nonlinear and Differential/ALgebraic equation Solvers (SUNDIALS) (Hindmarsh et al., 2005). In BioRT, the transport step is first solved with water by the preconditioned Krylov (iterative) method and the generalized minimal residual method (Saad and Schultz, 1986). All primary species in element i are then assembled in a local matrix and solved iteratively using the Crank–Nicolson and Newton–Raphson methods in CVODE (Bao et al., 2017).

Model verification

The BioRT module has been verified against CrunchTope under different transport and reaction conditions (Figs. S1–S7 in the Supplement). CrunchTope is a widely used subsurface reactive transport model (Steefel and Lasaga, 1994; Steefel et al., 2015) and is often used as a benchmark to verify other reactive transport models. Verification was performed under simplified hydrological conditions with 1-D column and constant flow rates such that it focuses on advection, diffusion, dispersion, and biogeochemical reactions. Specifically, three cases were verified. The phosphorus case that involves kinetics-controlled apatite dissolution and thermodynamics-controlled phosphorous speciation was first tested for solution accuracy of the bulk code that was inherited from the original RT-Flux-PIHM. Soil carbon and nitrogen processes were further verified for solution accuracy of the augmented BioRT module. Table S7 shows an average percent bias and Nash–Sutcliffe efficiency (NSE) of 1.1 % and 0.98, indicating a robust performance for a variety of solutes under different transport and reaction conditions. Note NSE ranges from −∞ to 1, with NSE=1 being the perfect fit (Moriasi et al., 2007).

Figure 5 Model structure, input, and output of BioRT-Flux-PIHM. The Flux-PIHM takes in watershed characteristics including topography (digital elevation model, DEM), land cover, shallow- and deep-zone properties, and meteorological forcing time series and solves for water storage and ground and soil temperature. BioRT takes in water- and temperature-related output from Flux-PIHM and additional inputs such as precipitation chemistry and shallow- and deep-water chemistry and biogeochemical kinetics parameters and solves for aqueous and solid concentrations in the shallow and deep zone and stream water. NADP stands for the National Atmospheric Deposition Program.

5 Model structure, data needs, and domain setup

5.1 Model structure

The model takes meteorological forcing time series as input and solves for water storage and soil temperature, along with other hydrologic and land-surface states and fluxes (Fig. 5). BioRT reads in the model output of water and temperature from Flux-PIHM and solves biogeochemical reactive transport equations. At the timescale of months to years that are typical for BioRT simulations, alterations in solid phase properties, including porosity, permeability, and reactive surface area, are considered negligible such that hydrological parameters remain constant with time.

5.2 Data needs

The code sets up the model domain based on watershed characteristics including topography, land cover, and shallow- and deep-zone properties (Fig. 5). When the model is used in a spatially distributed form, the model domain is set up using elevation, land cover, soil, and geology maps supplied by the user. A useful data portal is the Geospatial Data Gateway (https://datagateway.nrcs.usda.gov, last access: 22 May 2019). Another geospatial data source is HydroTerre (http://www.hydroterre.psu.edu/, last access: 22 May 2019), where users can obtain data on elevation, land cover, geology, and soil (Leonard and Duffy, 2013). Meteorological forcing data can be downloaded from the North American Land Data Assimilation Systems Phase 2 (NLDAS-2, https://ldas.gsfc.nasa.gov/nldas/v2/forcing, , last access: 22 May 2019). The vegetation forcing, i.e., leaf area index (LAI), can be obtained from MODIS (Moderate Resolution Imaging Spectroradiometer, https://modis.gsfc.nasa.gov/data, last access: 22 May 2019). Other vegetation properties (e.g., shading fraction, rooting depth) can be adopted from, for example, the Noah vegetation parameter table embedded in the Weather Research and Forecasting model (WRF; Skamarock and Klemp, 2019). Local measurements from meteorological stations and field campaigns (e.g., land cover, soil, geology) can be used in the model. Initial water and solid phase chemistry can be based on measurements or general knowledge of the simulated sites. The form of reaction rate laws can be defined in the input files and calibrated to reproduce field data. Reaction thermodynamics, mostly equilibrium constants, are from the geochemical database EQ3/6 by default (Wolery, 1992). These reaction parameters can be modified when necessary. The model outputs include aqueous and solid concentrations of shallow and deep zones and stream water.

Table 1 Example applications of BioRT-Flux-PIHM.

Note: transition state theory (TST) is a classic kinetic rate law for mineral dissolution and precipitation (Brantley et al., 2008) (Eq. S15); SOC stands for soil organic carbon.

5.2.1 Domain setup: spatially lumped and distributed domains

The domain can be set up at different spatial resolutions. A simple domain can be set up with two land grids representing two sides of a watershed connected by one river cell. This setup uses averaged properties without the need for spatial data. Alternatively, a complex domain can be set up to track “hot spots” of biogeochemical reactions using many grids with explicit representation of spatial details (e.g., topographic map, river network, land use map, soil and geology map, mineral distribution). The model domain can be set up using PIHMgis (http://www.pihm.psu.edu/pihmgis_home.html, last access: 22 May 2019), a standalone GIS interface for watershed delineation, domain decomposition, and parameter assignment (Bhatt et al., 2014). The same processes (e.g., hydrology, reaction network) can be setup in both types of spatial configurations. Auto-calibration is not built into the model, but a global calibration coefficient approach can be used to reduce parameter dimension and facilitate manual calibration. A typical model application requires 20 to 30 hydrological parameters to be calibrated. These parameters include land-surface parameters (e.g., canopy resistance, surface albedo), soil and geology parameters (e.g., hydraulic conductivity, porosity, Van Genuchten, macropore properties) (Shi et al., 2013). Reaction-related parameters (e.g., reaction rate constant, mineral surface area, Q₁₀, S_w,c, and n) are additionally needed for calibration, the number of which depends on the numbers of reactions involved in a particular system.

6 Model applications

The original RT-Flux-PIHM has been applied to understand processes related to the geogenic solutes of Cl and Mg at the Shale Hills watershed and for Na in a watershed on Chimborazo in the Ecuadorian Andes (Table 1). The new BioRT-Flux-PIHM has been demonstrated for understanding the dynamics of DOC and nitrate at Shale Hills and Coal Creek. This section will present one hydrology and two biogeochemical examples in the Susquehanna Shale Hills Critical Zone Observatory (SSHCZO), a small headwater watershed in central Pennsylvania, USA. The mean annual precipitation is approximately 1070 mm and the mean annual temperature is 10 ^∘C (Brantley et al., 2018). Soil carbon storage and respiration and nitrogen budget and fluxes have been studied in detail (Andrews et al., 2011; Hasenmueller et al., 2015; Weitzman and Kaye, 2018). Modeling work has been conducted to understand hydrological dynamics (Shi et al., 2013; Xiao et al., 2019), transport of the non-reactive tracer Cl, and the weathering-derived solute Mg (Bao et al., 2017; Li et al., 2017a).

Figure 6 (a) Hydraulic conductivity (K_satH) contrast controls the proportion of deep groundwater (Q_G). The cases of high ($K_{\text{satH}}^{\text{dp}}=\mathrm{2.6}$ m d⁻¹, red) and low conductivity ($K_{\text{satH}}^{\text{dp}}=\mathrm{0.22}$ m d⁻¹, blue) led to 26 % and 5.2 % of annual Q_G contribution to discharge (Q), respectively. (b) Deep-groundwater fraction as a function of K_satH contrast between the deep and shallow zone. The upper limit of the deep/shallow K_satH contrast was set to 1 as most watersheds have smaller K_satH in the deep zone than in the shallow zone. The two red and blue dots correspond to the two cases in (a).

6.1 Example 1: shallow and deep-water interactions

The model was set up using the spatially lumped mode with two grids and one river grid characterized by average land cover, soil, and rock properties based on previous work. The model assumed the dominant soil type (Weikert soil) at Shale Hills. The porosity of the deep zone was set to a tenth of the shallow soil porosity based on measurements of the groundwater aquifer (Brantley et al., 2018; Kuntz et al., 2011). In a headwater catchment like Shale Hills where the deep groundwater is most likely sourced from recharge, the deep-groundwater contribution to the stream can be primarily controlled by the hydraulic conductivity (K_satH) contrast between the deep and shallow zones (i.e., $K_{\text{satH}}^{\text{dp}}/K_{\text{satH}}^{\text{sl}}$). This is because the K_satH contrast determines the partitioning of infiltrating water between the shallow lateral flow and the downward recharge to the deep zone and then deep-groundwater flow. Two cases of high (red) and low (blue) $K_{\text{satH}}^{\text{dp}}$ were set up to showcase the control of K_satH contrast on deep groundwater (Fig. 6a). By changing the deep zone $K_{\text{satH}}^{\text{dp}}$ from 2.6 to 0.22 (m d⁻¹), the annual deep-groundwater (Q_G) contribution to discharge (Q) decreased from 26 % to 5.2 %, although the total stream discharge remains the same. This indicates that the changing $K_{\text{satH}}^{\text{dp}}$ mostly changes the flow partitioning between the shallow soil flow and deeper groundwater flow into streams.

Several additional cases were further tested to examine the relationship between deep-groundwater fraction (%) of discharge and K_satH contrast. Figure 6b shows that the deep-groundwater fraction rapidly increases with the increasing ratio of $K_{\text{satH}}^{\text{dp}}/K_{\text{satH}}^{\text{sl}}$, reaching a limit when the K_satH contrast is sufficiently high. The deep-groundwater contribution to the stream reaches ∼40 % when $K_{\text{satH}}^{\text{dp}}$ and $K_{\text{satH}}^{\text{sl}}$ are equal. In natural systems, we do see places, for example, karst formations, where groundwater contributes to more than 40 % (Hartmann et al., 2014; Husic, 2018). These places may have higher deeper conductivity than shallow soils due to the development of highly conductive conduits.

6.2 Example 2: nitrate dynamics in a spatially implicit domain

This example focuses on nitrate (${\mathrm{NO}}_{\mathrm{3}}^{-}$), a dominant dissolved N form in water (Weitzman and Kaye, 2018). The N processes at Shale Hills include atmospheric N deposition, soil N leaching, stream export, denitrification, and plant uptake (Fig. 7). Based on field measurements, the atmospheric deposition at the site is the dominant N input; N export via discharge is only a small fraction (2.5 %) of atmospheric N input. Most deposited N is tightly cycled by plants or lost to the atmosphere via denitrification.

Figure 7 Modeled nitrogen processes in Example 2. Atmospheric N deposition is the major N input; denitrification and plant uptake are the major N loss and sink. Export via discharge only occupies a small fraction.

The soil N leaching process was represented using a lumped reaction could represent the total rates of reactions including the decomposition of soil organic matter (SOM), nitrification, and rock weathering that generates ${\mathrm{NO}}_{\mathrm{3}}^{-}$. Its rate was assumed to depend on soil temperature and moisture and follows the equation r_leach=kAf(T)f(S_w), where r_leach [mol s⁻¹] is the leaching rate, $k={\mathrm{10}}^{-\mathrm{9.7}}$ [$\mathrm{mol}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{s}}^{-\mathrm{1}}$] is the leaching rate constant (Regnier and Steefel, 1999), and A [m²] is the surface area that represents the contact area between substrates and N-transforming microbe. The surface area was calculated based on SOM volume fraction [m³ m⁻³], specific surface area (SSA, [m² g⁻¹]), substrate density [g cm⁻³], and element volume [m³].

Figure 8 Stream nitrate dynamics and fluxes at Shale Hills in Example 2. (a) Three simulation scenarios with different processes are demonstrated here: transport only (dashed line, tran), transport+N reaction (gray line, tran+N react), $\text{transport}+\text{N reaction}+\text{plant uptake}$ (thick black line, $\text{tran}+\text{N react}+\text{upta}$), where N reactions include both nitrate leaching and denitrification (see Fig. 7); (b) nitrate concentration in precipitation and shallow and deep water; (c) nitrate fluxes and budget. Note that nitrate leaching was ignored in (b) due to its minimal flux as N deposition from rainfall was the dominant input (Weitzman and Kaye, 2018).

Denitrification converts ${\mathrm{NO}}_{\mathrm{3}}^{-}$ to N₂ gas under anaerobic conditions. Here this process was modeled by the Monod rate law with DOC as the electron donor, ${\mathrm{NO}}_{\mathrm{3}}^{-}$ as the electron acceptor, and with an inhibition term f(O₂) (Eq. S13). The reaction rate is

$$\begin{array}{}\text{(12)}& \begin{array}{rl}{r}_{\text{denitrification}}=& kA{\displaystyle \frac{C_{\text{DOC}}}{K_{m,\text{DOC}}+C_{\text{DOC}}}}{\displaystyle \frac{C_{{\mathrm{NO}}_{\mathrm{3}}^{-}}}{K_{m,{\mathrm{NO}}_{\mathrm{3}}^{-}}+C_{{\mathrm{NO}}_{\mathrm{3}}^{-}}}}\\ & \times f\left({\mathrm{O}}_{\mathrm{2}}\right)f\left(T\right)f\left(S_{\mathrm{w}}\right),\end{array}\end{array}$$

where $k={\mathrm{10}}^{-\mathrm{10}}$ [$\mathrm{mol}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{s}}^{-\mathrm{1}}$] is the denitrification rate constant (Regnier and Steefel, 1999) and half-saturation constants $K_{m,\text{DOC}}=\mathrm{15}$ [µM] and $K_{m,{\mathrm{NO}}_{\mathrm{3}}^{-}}=\mathrm{45}$ [µM] (Regnier and Steefel, 1999). For soil N leaching and denitrification, the SSA were, respectively, tuned as $\mathrm{1.6}\times {\mathrm{10}}^{-\mathrm{6}}$ and $\mathrm{7.5}\times {\mathrm{10}}^{-\mathrm{5}}$ [m² g⁻¹] to reproduce observed stream nitrate dynamics. The calibrated values were orders of magnitude lower than the lab-measured SSA of natural materials (e.g., SOM, 0.6∼2 m² g⁻¹) (Rutherford et al., 1992). Such discrepancies between calibrated effective reactive surface area (i.e., solid–water contact area) and lab-measured absolute surface area are consistent with observations of field–lab discrepancies in the literature (Li et al., 2014; Heidari et al., 2017; Wen and Li, 2017, 2018). The uptake rate constant was calibrated by constraining the partitioning of N transformation flux between denitrification and plant uptake by the ratio of 1:5, a value estimated from field measurements of gaseous N outputs (3.53 $\mathrm{kg}-\mathrm{N}{\mathrm{ha}}^{-\mathrm{1}}{\mathrm{yr}}^{-\mathrm{1}}$) and plant N uptake (18.3 $\mathrm{kg}-\mathrm{N}{\mathrm{ha}}^{-\mathrm{1}}{\mathrm{yr}}^{-\mathrm{1}}$) (Weitzman and Kaye, 2018). The uptake rate constant in the deep zone (>2 m in depth) was considered negligible (Hasenmueller et al., 2017). Groundwater nitrate was initialized as 0.43 mg L⁻¹, the average of measured groundwater concentration during 2009/10.

Figure 9 Attributes of Shale Hills in the spatially distributed mode in Example 3: (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); soil depths and SOC were interpolated using ordinary kriging based on field surveys (Andrews et al., 2011; Lin, 2006). The SOC distribution in (c) was further simplified using the high, uniform SOC (5 % $v/v$) in swales and valley soils based on field survey (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).

Temporal nitrate dynamics

Three cases were set up to understand and quantify the effects of different processes in determining nitrate dynamics (Fig. 8a). The transport-only case (dashed line, tran) simulates nitrate input from precipitation (at 1.4±0.96 mg L⁻¹, based on the 2009 data of NADP PA42 site) and N transport but without any reactions. It overestimated stream nitrate data (0.33±0.39 mg L⁻¹) throughout the year. The transport+N reactions case (gray line, tran+N react) has denitrification and soil N leaching processes but not plant uptake. These two reactions lowered the nitrate concentration slightly, as these two processes compensate each other in adding and removing nitrate from water. The $\text{transport}+\text{N reactions}+\text{uptake}$ case (thick black line, $\text{tran}+\text{N react}+\text{upta}$) have all processes. It significantly lowered the nitrate concentration, especially in April–May and October–December. Nitrate peaks from May to July, exhibiting comparable levels of high nitrate concentration (Fig. 8b). It is noticeable that the three cases almost overlapped at these overestimated short nitrate peaks, suggesting nitrate-rich precipitation may not be routed into the subsurface where denitrification and plant uptake could occur.

Although precipitation from April to August accounted for 70 % of the total simulation period, larger storm events in October contributed more to nitrate export. Deeper groundwater had higher nitrate concentration than shallow water because most plant uptake occurred in the shallow zone. The nitrate fluxes into the deeper zone however only contributed 26 % of stream nitrate export at the annual scale, due to the relatively small groundwater contribution (9.5 %) to the stream. Denitrification and plant uptake largely occurred during wet spring with leaves growing. Denitrification peaks often appeared after major storm events because wet conditions facilitate denitrification. Comparing the three outfluxes (Fig. 8c), nitrate export via discharge (red) was negligible compared to denitrification (blue) and plant uptake (green). At the annual scale, stream export accounted for 9.5 %, whereas denitrification and plant uptake took up 15 % and 75 % of deposited ${\mathrm{NO}}_{\mathrm{3}}^{-}$, respectively. In other words, as Nitrate enters this system via precipitation, plant uptake can play a significant role in reducing nitrate level, indicating the precipitated nitrate is tightly cycled in the system.

Figure 10 (a) Temporal dynamics of stream DOC concentration; spatial profiles of (b) shallow soil moisture, (c) hydrologically connected zones, (d) local DOC production rates r_DOC, and (e) shallow-water DOC concentration in May (wet), August (dry), and October (wet after dry) of 2009. The soil DOC and r_DOC were high in swales and valley with relatively high shallow-water and SOC content. August had the highest shallow-water DOC concentration compared to May and October because most DOC accumulated in zones that are disconnected to the stream.

6.3 Example 3: DOC production and export in a spatially distributed domain

This example showcases the application of BioRT-Flux-PIHM in a spatially distributed mode. This work has been documented with full details in Wen et al. (2020). Here we only introduce some key features and capabilities in the spatially distributed mode. The Shale Hills catchment was discretized into 535 prismatic land elements and 20 stream segments through PIHMgis based on the topography (Fig. 9a). The heterogeneous distributions of soil depth and solid organic carbon within the domain (Fig. 9b and c) were interpolated through ordinary kriging based on field surveys (Andrews et al., 2011; Lin, 2006). Other soil and mineralogy properties such as hydraulic conductivity, van Genuchten parameters, and ion exchange capacity were spatially distributed following intensive field measurements (Jin and Brantley, 2011; Jin et al., 2010) (https://criticalzone.org/shale-hills/data/, last access: 22 May 2019).

6.3.1 Temporal and spatial patterns of DOC production and export

The model outputs followed the general trend of stream DOC measurements with the model evaluation index NSE of 0.55 for monthly DOC concentration (Fig. 10a). NSE values greater than 0.5 are considered good performance for the monthly water quality model (Moriasi et al., 2015). The model reproduced high DOC values (∼15 mg L⁻¹) in the dry periods (July–September). The model enabled the identification of reaction hot spots. In May when soil water is relatively abundant, the valley and swales with deeper soils (Fig. 10b) are generally wetter compared to the hillslope and ridgetop, and are hydrologically connected to the stream (Fig. 10b and c). The distribution of local DOC production rate r_DOC and DOC concentration followed that of SOC (Fig. 10c) and water content (Fig. 10b). Low r_DOC in relatively dry planar hillslopes and uplands resulted in low soil water DOC. The average stream DOC (∼5 mg L⁻¹) reflected soil water DOC in the valley and swales.

In August, the hydrologically connected zones with high water content shrank to the vicinity of the stream and riverbed. With high temperature in summer, r_DOC increased 2-fold from May across the whole catchment while still exhibiting the highest values in the SOC-rich regions. Soil water DOC concentration increased by a factor of 2 because the produced DOC was trapped in low soil moisture areas that were not hydrologically connected to the stream. On the north side with low water content (Fig. 10b), the soil water DOC (∼7 mg L⁻¹ on average) accumulated more than on the south side (∼5 mg L⁻¹ on average). The high shallow-water DOC (∼10 mg L⁻¹) in the stream vicinity dominated the stream DOC in August.

In October, precipitation wetted the catchment again. The hydrologically connected zones expanded beyond swales and the valley to the upland hillslopes (Fig. 10c). The increase in hydrological connectivity zones favored the mixing of shallow-water DOC sourced from upland hillslopes (low DOC), swales, and valley (high DOC) into streams rather than only from the stream vicinity with high DOC in the dry August, leading to a drop in stream DOC.

Figure 11 Relationships between daily discharge (Q) and stream DOC concentration. With increasing Q, stream flow first shifted from the dominance of groundwater with low DOC at very low discharge to the predominance of organic-rich soil water from swales and valley at intermediate discharge. As the discharge increased further, the stream water switched to the dominance of high flow with lower DOC water from planar hillslopes and uplands, resulting in a dilution C–Q pattern at high discharge (from Wen et al., 2020, used with permission).

6.3.2 C–Q patterns

The DOC C–Q relationship showed a non-typical pattern with flushing first and transitioning into a dilution pattern, with an overall C–Q slope $b=-\mathrm{0.23}$ (Fig. 11). At low discharges ($<\mathrm{1.8}\times {\mathrm{10}}^{-\mathrm{4}}$ m d⁻¹) in the summer dry period, the stream DOC mainly came from the organic-rich swales and valley floor zones with high soil water DOC (Fig. 10e). With discharge increasing in wetter periods (i.e., spring and fall), the contribution from planar hillslopes and uplands with lower DOC concentration increased (Fig. 10e), leading to the dilution of stream DOC.

7 Discussion

BioRT-Flux-PIHM brings the reactive transport modeling capabilities to the watershed scale, enabling the simulation of subsurface shallow and deep flow paths and biogeochemical reactions influenced by hydroclimatic conditions and land-surface interactions. The expanded model capability of simulating bio-mediated processes such as plant uptake, soil respiration, and microbe-mediated redox reactions enables the simulation of carbon and nutrient cycling in the shallow subsurface. The inclusion of the deep-groundwater zone allows the exploration of the effects of subsurface structure on hydrological partitioning between shallow soil lateral flow and deep groundwater and their relationships with stream discharge. Although not shown here, the model can also simulate deeper groundwater coming from regional aquifers across the outer boundary. This can be particularly useful for watersheds of higher stream orders, where a large proportion of deep water may come from nearby regional aquifers.

The advantage and disadvantages of simple vs. complex models have long been debated in the modeling community (Fatichi et al., 2016; Li et al., 2021; Wen et al., 2021). The computational cost of solving a spatially distributed, nonlinear, multi-component reactive transport model is high, posing challenges for the application of ensemble-based analysis. With additional reactions and transport processes, the model includes more functions (such as reaction kinetic rate laws) and parameters (e.g., reaction rate constants, surface area) than hydrological models. The complexity brings in issues of equifinality, uncertainty, and data demands (Beven, 2006; Kirchner et al., 1996). These issues will persist even though reactive transport models will be constrained by additional chemical data.

It is in this spirit of balancing the cost and gain that we present both spatially distributed and lumped modes for the BioRT model (Li et al., 2021). Compared to the distributed version, the spatially implicit model requires less spatial data and is computationally inexpensive. It can assess the average dynamics of water and solute dynamics and focus on the interacting processes without resolving spatial details. The lumped approach can accommodate basins with low data availability, and it can be easier for students to learn. In contrast, spatially explicit representations enable the exploration of the hot spots (e.g., swales and riparian zones with high soil water DOC concentrations in Fig. 10e) and their contribution to stream chemistry at different times. Spatial heterogeneities in watershed properties (e.g., soil types and depth, lithology, vegetation, biomass, and mineralogy) are ubiquitous in natural systems and are challenging to resolve. A general understanding of the linkage between local catchment features and catchment-scale dynamics (e.g., stream concentration dynamics and solute export pattern) is often lacking. The spatially distributed model provides a tool to explore these questions. Ultimately, the choice of the model complexity level depends on research questions that the model is set to answer and the available data. At the end, we all need to balance cost and gain when deciding to use a simple or complex model, striving to be “simple but not simplistic” (Beven and Lane, 2019).

8 Summary and conclusion

This paper introduces the watershed-scale biogeochemical reactive transport code BioRT (short for BioRT-Flux-PIHM). The code integrates processes of land-surface interactions, surface hydrology, and multi-component biogeochemical reactive transport. The new development enables the simulation of (1) biotic reactions including plant uptake, soil respiration, and microbe-mediated redox reactions and (2) surface water interactions with groundwater from a deeper subsurface that still interacts with streams. BioRT has been verified against the widely used reactive transport code CrunchTope for soil carbon, nitrogen, and phosphorus processes. It has been applied to understand carbon, nitrogen, and weathering processes at Shale Hills in Pennsylvania, Coal Creek in Colorado, and the Chimborazo watershed in the Andes in Ecuador. Here we showcase the modeling capability of surface–groundwater interactions and reactive transport processes relevant to nitrate and DOC in Shale Hills in two simulation modes. One is in a spatially lumped mode using averaged properties and another is in a spatially distributed mode with a consideration of spatial heterogeneity. Results show that the deep-groundwater flow that interacts with the stream is primarily controlled by the hydraulic conductivity contrast between the shallow and the deep zones. Biogeochemical reactions in shallow soil primarily determine the stream water chemistry under high-flow conditions. The spatially lumped method with two lumped grids can capture the temporal dynamics of average behavior and mass balance; the spatially distributed running mode can be used to understand the spatial dynamics and to identify hot spots of reactions. The code can be used for biogeochemical reactive transport simulations in watersheds under diverse climate, land cover, and geology conditions.