A watershed (also drainage basin, river basin, or catchment) is defined as “… the area that topographically appears to contribute all the water that passes through a specified cross section of a stream (the outlet)” (Dingman 2015). In this chapter, I choose to use the term “watershed” as it is a broadly used one; it should be understood more as small watersheds or catchments. Watersheds are the fundamental units that support river networks, the blood vessels at Earth's surface ultimately draining into the ocean.

Watersheds are complex hydro-biogeochemical reactors. They receive water, mass, and energy, transport them to distinct compartments, and transform them into various forms (Fig.1). Hydrological processes partition precipitation to the atmosphere, to the ground surface, and to the subsurface, eventually entering streams. Similarly, plants translate sunlight, water, CO2, and nutrients into organic matters (leaves, stems, roots) that fall and deposit in soil. As water routes through soils, it interacts with roots, microbes, and reactive gases (i.e., CO2 and O2), releases solutes, and ultimately transporting them out of watersheds. The water flow (discharge) and solute concentrations measured at stream outlets therefore reflect convoluted signature of ecohydrological and biogeochemical coupling. The process interactions and feedbacks are dictated not only by hydroclimatic forcing but also the architecture of watersheds, in particular the above-ground characteristics such as land cover and surface topography, as well as the below-ground structure including soil depth, soil type, geology, and root architecture. The external forcing and internal idiosyncrasies dictate the magnitude, timing, and spatial distribution of water flow and chemistry (Chorover et al. 2011; Brooks et al. 2015), giving rise to non-linear emergent behaviors that are unique of watershed reactive transport processes.

Understanding feedbacks and non-linearity requires process-based models that integrate multiple interacting processes. Such integration however does not come easily with traditional boundaries of hydrology versus biogeochemistry relevant disciplines. As will be discussed later in this chapter, relevant model development has been advancing along two separate lines: hydrology models that solve for water storage and fluxes at the watershed scale and beyond (Fatichi et al. 2016), and reactive transport models (RTMs) that center on aqueous and solid concentration changes arising from transport and multi-component biogeochemical reactions typically in “closed” groundwater systems without much interactions with “open” watersheds directly receiving precipitation and sun light (Steefel et al. 2015; Li et al. 2017b). This comes along with a history of hydrologists often trained as physicists studying fluid mechanics, and biogeochemists typically grow up as geologists, chemists, or environmental engineers.

There are however considerable needs to reach beyond disciplinary boundaries and integrate the two lines to develop watershed reactive transport models, not only to gain fundamental understanding of watersheds as complex systems, but also to solve today's pressing global environmental challenges. Natural systems such as watersheds do not set up artificial disciplinary boundaries; processes relevant to different disciplines all occur simultaneously. Research questions at the interfaces or “ecotones” of different disciplines can stimulate answers that shed light on puzzling observations arising from non-linear emerging behaviors. In addition, as the pace of climate change and human perturbation accelerates, a wide spectrum of water-related hazards (e.g., hurricanes, flooding, droughts, melting glaciers) loom (IPCC 2013; Fan et al. 2014; Zhao et al. 2016; Roque-Malo and Kumar 2017). Large hydrological events bring out excessive water and disproportionally large pulses of “stored” contaminants that deteriorate aquatic and ecosystem health (Raymond and Saiers 2010; Huntington et al. 2016). On the other hand, droughts induce water-borne diseases (Perez-Saez et al. 2017) and impair water quality (Ejarque et al. 2018). Problems related to excessive nutrient export, including eutrophication and hypoxia in rivers, lakes, and coastal areas worldwide, have lingered for decades and will continue to do so, calling for advanced tools for prediction and for management (Royer et al. 2006; Seitzinger et al. 2010; Van Cappellen and Maavara 2016; Van Meter et al. 2017). As such, it is essential to possess the capabilities of forecasting not only for water quantity (flow), which is well underway with models such as the National Water Model (Cosgrove et al. 2016) (https://water.noaa.gov/about/nwm), but also for water quality (chemistry).

Although solute transport and reaction modules have been developed as add-ons to hydrological models at the watershed scale (e.g., Arnold and Soil 1994; Donigian Jr et al. 1995; Santhi et al. 2001), they typically have relatively crude representation without rigorous incorporation of reaction thermodynamics and kinetics theory in geochemistry and biogeochemistry. Reactive transport models have been brought to the watershed scale only recently (Yeh et al. 2006; Bao et al. 2017). The goal of this chapter is to illustrate this development and their promises in facilitating the understanding of watershed processes. To do so it is important to bring together the hydrology and biogeochemistry communities on fundamental concepts and processes. I therefore first lay out concepts in watershed hydrology and biogeochemistry, as surface hydrological processes may sound foreign to geochemists; similarly, biogeochemical processes may appear alien to hydrologists. I will then illustrate examples of hydrological and biogeochemical coupling and the use of watershed reactive transport models in offering mechanistic understanding. The chapter ends at thoughts of existing knowledge gaps and future directions that can potentially be tackled with watershed RTMs.


Watersheds differ from groundwater systems in that they are “open” systems receiving energy, rainfall and snow from the atmosphere. This brings in significant temporal dynamics or non-stationary that depends on hydroclimatic conditions. In addition, watersheds are spatially heterogeneous not only with distinct subsurface zones of soil, weathered rock, and parent rocks, but also with topography, land cover (e.g., barren soil versus forests, grassland, urban), and river networks. Flow complexity therefore originates from temporal dynamics, variations in subsurface properties such as porosity and permeability (e.g., highly permeable soils versus tight rocks), as well as slopes, land cover, and aspect (sunny or shady sides) (Buttle 1998; Bishop and Seibert 2015). This section does not intend to cover the full spectrum of watershed processes. Rather, I will introduce fundamental aspects that are relevant to later discussions on biogeochemical processes.

Water balance

A myriad of processes determine where precipitation (P, rain + snow) goes and how it partitions into different parts of a watershed. It can evaporate through soils or transpire through plants, returning to the atmosphere. It can be intercepted by vegetation leaves or fall on the ground; it can flow directly on the ground surface to the stream or infiltrate into soils. Within soil it can penetrate vertically through the unsaturated zone, form interflow toward the stream or channel into macropores. Some soil water may cross the soil–weathered bedrock interfaces and flow downward into groundwater aquifers, which may eventually re-enter rivers and streams, although not necessarily in the same watershed.

Here I will only coarsely delineate these processes without getting into the intricacies of these processes. Considering a watershed as a control volume or a water processing box, we have precipitation (P) as the input and evaporation (E), and transpiration (TR) and discharge (Q) as the output (Fig. 1). This leads to a simple water balance equation (e.g., Kirchner 2009):


where S, the water storage, is in units of depth (e.g., mm of water, or volume of water per unit area, L), and P(t), ET(t) (= E(t) + TR(t)), and Q(t) are the rates of precipitation, evapotranspiration, and discharge, respectively, in units of depth per time (e.g., mm of water per hour, L/t, essentially area-normalized water volume per time). All terms in Equation (1) are a function of time and represent average values over an entire watershed. The equation assumes that water loss into aquifers is negligible. Equation (1) essentially follows the mass conservation principle and states that changes in water storage in a watershed depends on the temporal dynamics of water input (P) and outputs (ET and Q). Solving this equation needs constitutive equations rsuch as storage-discharge relationships. An extensively used example is the power law form Q = aSb (Horton 1936) (Brutsaert and Nieber 1977; Wittenberg 1999). Whereas precipitation is driven by meteorological forcings, a plethora of watershed features determine the partitioning between ET and Q, as explained below.

Evapotranspiration (ET).

Evapotranspiration is the summation of evaporation (from open water, bare soil, and vegetated surfaces), transpiration (from within plant leaves), and sublimation from ice and snow surfaces. ET is regulated by land surface interactions that dictate how land surface and atmosphere exchange energy and water. Another important concept is the climatic potential ET (PET), defined as “the rate at which evapotranspiration would occur from a large area completely and uniformly covered with growing vegetation with access to an unlimited supply of soil water and without advection or heat-storage effects.” (Dingman 2015). Although somewhat ambiguous in that different types of vegetation have different transpiration demand, PET is considered as a measure of “drying power” of a watershed influenced by both climate and land cover.

A widely used relationship, Budyko equation, relates the long-term ratio of evaporative index ET/P to the aridity index PET/P. Various forms of Budyko equation exist (Pike 1964; Budyko 1974; Fu 1981; Gentine et al. 2012). An example is as follows:


Here β is dimensionless and modifies the curvature of ET as a function P and PET (Fig. 2). This equation says (empirically) that the long-term average of ET is a function of water supply (P) and available energy estimated as the evaporative demand for water by the atmosphere (PET). Under arid conditions when PPET, ET converges toward P, implying that ET is limited by water supply (Fig. 2, water limit line). Alternatively, under humid conditions when PPET, ET is limited by atmospheric demand and collapses toward PET (Fig. 2, energy limit line). Budyko's work and other follow-up studies have shown that under steady state conditions without significant groundwater inputs, losses, or storage changes, hydrological systems often operate close to either energy or water constraints. In addition to climate influence, vegetation modulates ET via dynamical relations between soil moisture in the rooting zone, evapotranspiration rates, and rainfall interception, therefore regulating Budyko's original energy-based prediction of the evaporative fraction (Porporato et al. 2004; Gentine et al. 2012). ET and PET can be estimated based on ground-based meteorological observations and remote sensing data from the MODerate Resolution Imaging Spectroradiometer (MODIS) (Mu et al. 2007). Precipitation data are available from various databases (e.g., WorldClim (Hijmans et al. 2005), NLDAS-2 (the North American Land Data Assimilation Systems phase 2 (https://ldas.gsfc.nasa.gov/nldas/NLDAS2forcing.php)) or measurements from local weather stations.

Discharge (Q).

Assuming negligible water loss to deep aquifers, the net difference between P and ET is discharge (stream flow) (Fig. 1). The generation of stream flow is an integrated outcome of surface runoff, soil interflow, and shallow groundwater contributing to streams. Surface runoff (also known as overland flow) occurs when excess water flows over the Earth's surface because soil is fully saturated, or rain arrives more quickly than the rate of infiltration, or impervious areas (roofs and pavement) cannot infiltrate most of water. The generation of shallow soil interflow depends on multi-phase flow dynamics between water and air in the unsaturated zone (via Richard's equation), and the formation of perched water table and saturated zone where water is sufficient to flow toward a stream (Darcy's law). The formation of such flow depends on soil properties as described by the water retention curve (Van Genuchten 1980):


where θ(h) is the water retention ([L3 L−3]), |h| is the water pressure (or water head, L), θs and θr are the saturated and residual water content ([L3 L−3]), respectively; by definition is the same as porosity; α and n are parameters defining the shape of the water retention curve and are generally determined by soil properties. Two example water retention curves are shown in Figure 3. Generally sandy soils tend to leak water quickly and have a “flatter” shape, whereas the water holding capacity of clay (θr) is higher than sandy soil. The porosity of clayey soil (θs) is much lower than that of sandy soil, therefore leading to a smaller water content range (θs − θr) that drives flow.

Note that stream water can also come from the groundwater below the soil–weathered bedrock interface, which often forms base flow under dry conditions. Note that this may be different from regional groundwater in deep subsurface that is not connected to the stream. The contribution of groundwater to streamflow vary depending on both climate and watershed characteristics such as topography (e.g., slope), land cover, and lithology (Beck et al. 2013; Welch and Allen 2014). Discharge essentially sums up three major components: surface runoff, shallow soil interflow, and shallow groundwater flow that is connected to streams. Their contributions to stream water vary depending on hydrological regimes. Under dry conditions, stream flow is often dominated by groundwater as the base flow. Under wet conditions, stream flow is typically composed primarily of soil interflow and / or some surface runoff. These end member water sources have very different water chemistry such that their proportions determine the stream water composition at different times, as will be discussed in the section Watershed response to hydrological changes: Concentration discharge relationships

Water storage (S).

Equation (1) prescribes that water storage, or water content, is determined by the dynamics of water input (P) and output (ET and Q). Direct measurements of water storage is challenging and often require large field campaigns (Lin et al. 2006; Trompvan Meerveld and McDonnell 2006). As a result, it is often indirectly estimated based on discharge and / or tracer measurements (Seibert et al. 2011). These observations have led to several major understandings. First, although it is challenging to accurately quantify the total storage, we know that it can be grouped into dynamic storage that is responsive to transient hydrological conditions, and passive storage that is much less responsive (McNamara et al. 2011). Dynamic water store leaves a watershed rapidly as “young” water whereas passive water storage lingers for substantially longer time (Barnes and Bonell 1996; Dunn et al. 2010; Benettin et al. 2015a; Sprenger et al. 2018). Second, it has been shown that dynamic water storage is often dramatically smaller (often > an order of magnitude) than the total storage inferred from average soil porosity (Birkel et al. 2011). Such immense reduction has been attributed to catchment characteristics that lead to low connectivity, a measure of the linkage between disparate regions of the hillslope to the stream via subsurface water flow (Jencso et al. 2009; Bracken et al. 2013). Tracer or isotope signatures often bear signatures of the relative size of the two different pools. A much smaller variation of water isotope signature in streams compared to that of rainfall indicate significant damping and large mixing volumes and water storage. Kirchner (2009) demonstrated that changes in water storage often strongly relate to changes in stream discharge, such that a storage–discharge relationship can be used directly to reconstruct the time series of precipitation and evapotranspiration in some watersheds. This develops from a tradition in catchment hydrology using estimates of storage change in linear or nonlinear reservoirs from water balance calculations and is the foundation for streamflow simulation in many conceptual rainfall–runoff models (Bergstrom 1976).

Modeling watershed hydrological processes

Rainfall-runoff models have been used extensively as conceptual, lumped models to estimate streamflow hydrograph for rainfall events and are well documented in textbooks (Beven 2011; Dingman 2015). The development of process-based watershed models date back fifty years (Crawford and Linsley 1966; Freeze and Harlan 1969). Equation (1) is essentially a simple one box model based on water balance and can be solved for water storage and water fluxes. A slightly more complex model with two boxes (an upper and lower box) has also been formulated based on the conceptualization of the two types of water storage to represent the dynamics of “fast and small” reservoir versus the passive, “slow and large” reservoir (Benettin et al. 2015a; Kobierska et al. 2015; Kirchner 2016b). This two-box model has been shown to better capture the temporal non-stationary of streamflow than the one-box model (Kirchner 2016b).

The hydrology community has developed and utilized models in a spatially distributed manner (distributed models) for more than five decades (Freeze and Harlan 1969; James 1972; Jarboe and Haan 1974; Bergstrom 1976; Abbott et al. 1979, 1986; Beven 1989; Quinn et al. 1991; Singh 1995; VanderKwaak and Loague 2001; Gan et al. 2006; McDonnell et al. 2007; Qu and Duffy 2007; Kumar et al. 2009; Therrien et al. 2010; Fatichi et al. 2016). With the integration of surface energy balance, the recent introduction of land surface processes into hydrological models marks a new advance toward more accurate representation of evapotranspiration (Maxwell and Miller 2005; Shi et al. 2013; Davison et al. 2014). These hydrology models have been brought to continental scales for water forecasting and earth system modeling (Clark et al. 2016; Fan et al. 2019). For example, the prediction of stream flows in the continental U.S. is now accessible through the National Water Model (https://water.noaa.gov/about/nwm).

These spatially distributed models have been criticized for over parameterization, equifinality, or parameter non-uniqueness, i.e., different sets of parameters can produce the same model output (Beven and Freer 2001). Along similar lines, model parameters calibrated for one watershed are often not transferrable to other watersheds, or even to periods other than calibration data in the same watershed (van der Linden and Woo 2003; Heuvelmans et al. 2004; Li et al. 2012; Smith et al. 2016). This issue of parameter transferability prevents the prediction of stream flow in ungauged river basins where data are not available (Sivapalan 2003; Hrachowitz et al. 2013). The hydrology community generally recognizes these challenges (Sivapalan 2017). In many cases however spatially distributed models are our only resort, particularly when spatial and temporal patterns are important drivers of water flow (Fatichi et al. 2016).

Here I use the Penn State Integrated Hydrology Model (PIHM) as an example of a distributed surface hydrology model to illustrate the governing equations (Qu and Duffy 2007). The land surface interaction processes are simulated in the Flux module using the Noah Land Surface Model (LSM) that is coupled to PIHM (Chen and Dudhia 2001; Ek et al. 2003; Shi et al. 2013). A variation of the PIHM code, FIHM, has more complex subsurface structure representation (Kumar et al. 2009). A suite of modules for processes of interests to multiple disciplines (Duffy et al. 2014) have been developed, including LE for landscape evolution (Zhang et al. 2016b), Cycles for agroecosystem functioning (Stöckle 2003; Kemanian 2010), and RT for biogeochemical reactive transport (Bao et al. 2017), all with Flux-PIHM being the backbone of the hydrology code.

PIHM applies the semi-discrete finite volume scheme by integrating the partial differential form of governing equations over a three dimensional control volume, thereby converting them to ordinary differential equations (ODEs) of individual prismatic elements. It solves for five water storages in each prismatic element i: above-ground storage in vegetation canopy, snow, and ground surface (water on land surface that forms surface runoff), and below-ground storage in unsaturated and saturated zones (Qu and Duffy 2007). Here I show the equations of water storage for the last three items. The equation for the surface water storage in element i is as follows:


where Ssf is the water storage above the land surface (L). The rates of net precipitation, i.e., precipitation not intercepted by canopy (pnet), infiltration from land surface to unsaturated zone (I), evaporation from surface water (Esf), and lateral surface flow from element i to j (qsf,ij) are all normalized by the base area of the finite volume [L3/(L2t)]; Ni,1~3 is the index of the neighboring elements of i. The qsf,ij is calculated based on a diffusion wave approximation of the 2D St. Venant equation (Gottardi and Venutelli 1993).

The governing equations for the water storages in the unsaturated and saturated zones of element i are as follows:


where Su and Ss are the water storage in the unsaturated and saturated zones, respectively; R is the recharge rate downward from the unsaturated into the saturated zone (L/t); Eu is the evaporation from the unsaturated zone; Et is the transpiration (L/t); qbedrock is the downward flow rate into bedrock [L/t]; qs,ij is the lateral flow rate [L/t] from element i to j in the saturated zone, which is calculated using Darcy's law.

This system of ODEs for individual elements is assembled to form a global ODE system, which is solved by an ODE implicit solver. Note that if we add Equations (4)(6) and the unpresented equations for snow and tree canopy, the vertical water flux terms between different vertical zones (surface, unsaturated, saturated) within an individual element (e.g., infiltration, recharge) will be cancelled and we are left with an equation for total water storage with changes depending on the total input and output of the element i. If we add these equations for all elements in the whole watershed, the lateral fluxes between elements (qij) cancel out and we will have a mass balance equation for the whole watershed that is essentially Equation (1).

An example of hydrological processes in the Shale Hills catchment

Here I use an example calculation using Flux-PIHM in the Shale Hills, a first order catchment in the Susquehanna Shale Hills Critical Zone Observatory (SSHCZO) in central Pennsylvania, U.S.A (Brantley et al. 2018; Li et al. 2018). The mean annual temperature is 10°C with a mean annual precipitation of ~ 1000 mm. Extensive field studies have characterized the topography, hydrological properties and mineral composition (Lin 2006; Jin et al. 2010, 2011b; Ma et al. 2010; Brantley et al. 2013b).

Figure 4 shows the partition of precipitation into several major components. From April 1st to Dec. 31st 2009, the precipitation was 0.9 meter. Assuming no deep groundwater draining bypassing the stream, the model estimates that ~40% of precipitation contributes to stream discharge, and ~ 60% to ET (Fig. 4A). SSHCZO is hydrologically responsive with stream discharge closely following intensive precipitation (Fig. 4B). In general, large rainfall events lead to more discharge and less ET. A large rainfall event occurred on Oct. 24th resulting in the highest discharge in 2009. Surface runoff showed up as short-lived pulses during rainfalls, followed by soil water lateral flow (interflow). Annually this lateral flow contributes to ~ 70–80% of the discharge, whereas surface runoff contributes ~ 10–20%. Groundwater from below the bedrock is not calculated (in this version of the code) and is assumed constant across the year contributing a total of ~ 10% to the stream, based on field observations (Jin et al. 2014; Sullivan et al. 2016). Although the precipitation is relatively invariable over the year, water storage is lower in the summer due to the higher ET as a result of higher temperature and the presence of tree leaves. This is consistent with the observed water level drop at groundwater monitoring wells (Sullivan et al. 2016).


Imagine that you walk in a forest on a beautiful autumn afternoon. The sun light and specks of blue sky sneak in between green and orange leaves. You hear the satisfying crunch as you step on the thick layer of dry leaves and branches on the ground. You notice that the shallow soil is the territory of living things and organic matter (OM, roots, fallen leaves, and microbe). The stabilization of OM can occur via sorption on clays or separation from water and microbe (Rasmussen et al. 2018). With exposure to water and other reactive chemicals, soil OM (SOM) can also decompose, followed by the release and reuse of organic carbon, nitrogen, and phosphorous, often called biogenic solutes (Schlesinger and Bernhardt 2013) (Fig. 5), as well as other elements. SOM decomposition is a redox reaction where OM is oxidized and electron acceptors (e.g., O2, NO3, Fe(III) oxides, sulfate) are reduced. These reactions have long been recognized as complex and mostly microbe-mediated (von Lützow et al. 2007; Conant et al. 2011; Lehman and Kleber 2015).

Organisms in shallow soils modulate chemical weathering and production of geogenic solutes by transferring water and nutrients among solid, aqueous, and biological reservoirs through root growth and exudation, and through associated fungal and microbial activities. The partial pressure of soil CO2 (pCO2), a product of soil respiration, is typically 1–2 orders higher than the atmospheric level in vegetated lands (Romero-Mujalli et al. 2018) and therefore plays a major role setting the acidity of soil waters (Drever 1989, 1997; Heidari et al. 2017). The presence of O2 drives oxidative weathering of Fe(II)-bearing minerals such as biotite and pyrite, creating fractures and porosity that allow water infiltration and weathering of other minerals (Molins and Mayer 2007; Buss et al. 2008; Bazilevskaya et al. 2013).

Water plays a central role in all these processes. It serves as a medium for the transport of reactants and products, sets the soil–water–microbe contact area, and influences the levels of reactive gases (e.g., CO2, O2). The thin soil mantling the Earth's surface is hydrologically more responsive to meteorological events (rainfall, snow, and temperature) compared to weathered and parent bedrocks at depth that harbor groundwater. Much remains to be learned about feedbacks among water flow, transport, weathering, and microbe-mediated processes, and how they couple to the long-term evolution of the critical zone and to global elemental cycles. Interested readers are referred to literature on the global cycle of C (Berner 1999; Falkowski et al. 2000; Wang et al. 2010), N (Galloway et al. 2004; Fowler et al. 2013), and P (Filippelli 2008; Oelkers et al. 2008). Here I briefly introduce important microbe-mediated biogeochemical reactions and chemical weathering.

Biogeochemical reactions

Carbon (C).

Carbon exists in both inorganic and organic forms. Inorganic carbon includes CO2 gas, dissolved or total inorganic C (DIC, TIC, in water), and carbonate rocks (solids). Organic C occurs in ecosystems (plants, animals, other living things), in water (dissolved organic C (DOC)), and in soil (soil organic C (SOC)). SOM represents the largest terrestrial pool of organic carbon (Stockmann et al. 2013). It can decompose partially into smaller organic molecules that dissolve in water, i.e., DOC, or oxidize completely into CO2 gas or dissolved inorganic carbon (DIC). With coexisting divalent cations (e.g., Ca, Mg), DIC can also precipitate and become carbonate minerals. Hence soil C decomposition can release CO2 back into the atmosphere and changes CO2 level (Davidson and Janssens 2006), or releases DOC and DOM to surface water. These processes occurs in soils and also as dissolved carbon transport laterally from land to ocean. They represent a small term in regional carbon budgets compared with gross primary production and respiration and are traditionally thought as negligible in the global carbon cycle. Recent literature however has revealed that these processes are important in assessing net changes in terrestrial C storage and carbon cycle (Battin et al. 2009; Regnier et al. 2013), challenging the previous paradigm and suggesting the need to include them in global carbon models.

Nitrogen (N).

The SOM decomposition releases organic nitrogen (R-NH2), after which complex biotic and abiotic reactions ensue. As shown in Figure 5, organic N can become mineralized and transform into ammonia, which can further become oxidized into forms of high oxidation states through a chain of reactions. Inorganic N is notorious for having 7 oxidation states with ammonium (NH4+ and NH3) at its lowest oxidation state of −3 up to nitrate (NO3, N2O5) at +5, with various forms in between (N2, N2O, NO, N2O3 (NO2), NO2). Some of the gaseous forms emit back to the atmosphere and are in fact greenhouse gases (Saha et al. 2017; Maavara et al. 2018). Denitrification requires anoxic conditions and does not occur as much in soils owing to the pervasive presence of O2; it can become prevalent however under extremely wet conditions and in O2-depleted groundwater systems. The significance of nitrogen resides in its role not only as a limiting element for the growth of organisms but also as a component in fertilizers that feed crops and ultimately human population. Its export from agriculture lands into aquatic systems cause water quality problems such as eutrophication worldwide (Rabalais et al. 2001; Davies and Neal 2007; Conley et al. 2009).

Phosphorous (P).

Earth's biological systems have depended on P since the beginning of life (Nealson and Rye 2003). Phosphate is central to the functioning of the adenosine triphosphate (ATP), the keystone for metabolism and the most abundant biomolecule in nature (Schlesinger 1997); the phosphate ester bridges bind the helix strands of DNA; all cell membranes are made up with phospholipids. In soils, P can be in organic form (e.g., leaves), sorbed (on fine soil particles), dissolved in water, or in solid forms as P-containing minerals (Fig. 5). The transformation between different forms occurs through various bio-mediated or abiotic reactions. Unlike C and N that are abundant in the atmosphere, all P on Earth comes from the minuscule quantities of P-containing rock on Earth's crust (0.09 wt%) (Filippelli 2002). The most abundant P-containing mineral is apatite Ca5(PO4)3(F, Cl, OH). Once liberated via rock dissolution, P is mostly locked in organic forms (e.g., inositol phosphates, phospholipids, nucleic acids, phosphoproteins). It is barely soluble so it binds on and transports together with soil particles in the form of orthophosphate or pyro-diphosphate. To feed the ever-increasing population in Earth, human has mined rock and manufactured P-containing fertilizer. Excessive leaching from agriculture soils, together with deforestation and soil erosion, have accelerated P delivery from the deep subsurface to the surface water and the ocean, exacerbating the issue of P limitation on Earth.

Writing microbe-mediated reactions in light of microbial energetics.

Microbe-mediated redox reactions involve organic carbon as electron donors that become oxidized and electron acceptors that are reduced. This process produces energy that sustains and grows microorganisms. Electron acceptors in natural environments include, for example, O2, nitrate, iron-containing and manganese-containing minerals, and sulfate (Fierer et al. 2003; Dunn et al. 2006; Philippot et al. 2009). This order represents the descending amount of energy that microorganisms can derive from mediating these reactions, often called biogeochemical redox ladder (Van Cappellen and Gaillard 1996; Thullner et al. 2005; Borch et al. 2010). That is, per unit organic carbon oxidized to the same product, microorganisms can glean more energy by reducing O2 compared to reducing NO3.

To take into account of the mass and energy balances of these reactions, one will need to consider microbial energetics and thermodynamics and microorganisms as part of reaction products (VanBriesen 2002; Xiao and VanBriesen 2006; Xiao, 2008). Microorganisms cannot grow if the energy derived from these redox reactions are not sufficient (Jin and Bethke 2003; Jin and Bethke 2007). The amount of energy needed for microbial cell maintenance is relatively similar, such that microorganisms using redox reactions that are high in the biogeochemical redox ladder (e.g., oxidation) have more left over energy to channel into cell synthesis and growth. Therefore, values of fs, the fraction of energy used for microbial cell growth, decrease as we go down the biogeochemical redox ladder from aerobic oxidation to sulfate reduction. That is, less microbial cells are produced as we go down the redox ladder with the same amount of electron donor.

Details of how to write such reaction equations are given in Chapter 2 of Rittmann and McCarthy (2001). Table 1 shows a few example redox reaction equations with microbial cells as a product and acetate as the electron donor (Li et al. 2009, 2010; Cheng et al. 2016) using the chemical formula C5H7O2N for microbe. Note that R is the overall reaction equation incorporating the half reactions of electron donor (Rd), electron acceptor (Ra) and cell synthesis (Rc). The values of fe and fs are the fractions of energy that partition into energy consumption reaction (anabolic pathway that consumes energy for cell maintenance) Re = RaRd and cell synthesis reaction (the catabolic pathway that generate new microbial cells) Rs = RcRd, respectively. The overall reaction R = feRe+ fsRs= fe(RaRd) + fs(RcRd), where fe + fs = 1.0. If your do not care about microbe production and cell growth, the reaction can be written simply as R=RaRd. In that case, the aerobic reaction is simply 14O2+18CH3COO=14HCO3+18H+.

Microbe-mediated reaction kinetics.

SOM is often conceptualized and modeled as pools with different decomposition rates and turnover times (Ostle et al. 2009; Thornton et al. 2009). An extensively used three-pool model includes a readily degradable (labile) pool with residence times less than 5 years; a slowly degrading pool with residence times of decades; and a relatively stable pool, with residence times between 103–105 years (Trumbore 1993; Trumbore et al. 1995; Marin-Spiotta et al. 2009). The kinetics of microbe-mediated reactions can be described by the general dual Monod rate law, reflecting the need for both electron donor and acceptor in these reactions (Monod 1949):


Here µmax is the rate constant (mol/t/microbe cell), CC5H7O2N is the concentration of microorganisms (cells / L3), CD and CA are the concentrations of electron donor and acceptor (mol/L3), respectively. The Km,D and Km,A are the half-saturation coefficients of the electron donor and acceptor (mol/m3), respectively; they are the concentrations at which half of the maximum rates are reached for the electron donor and acceptor, respectively. If an electron donor or acceptor is not limiting, it means that CDKm,D or CAKm,A, so that the term CDKm,D+CDorCAKm,A+CA is essentially 1, lending to a rate that only depends on the abundance of microorganisms or one of the chemicals.

In natural subsurface where multiple electron acceptors coexist, the biogeochemical redox ladder dictates the sequence of redox reactions. That is, aerobic oxidation occurs before denitrification, which in turn occurs before iron reduction, because microorganisms harvest more energy in aerobic oxidation than denitrification. Inhibition terms are introduced to account for the sequence of redox reactions as follows:


Here KI,H is the inhibition coefficient for the inhibiting chemical H. The inhibition term is 1 (not inhibiting) only when CHKI,H. In a system where oxygen and nitrate coexist, which is common in agriculture lands, aerobic oxidation occurs first before denitrification. The denitrification rates can be represented by:


Here CNO3 is the concentration of nitrate, KI,O2 is the inhibition coefficient of O2, or the O2 concentration at which it inhibits denitrifcation. This rate law ensures that denitrification kicks in substantially only when O2 is depleted to a level that CO2KI,O2, such that the term KI,O2KI,O2+CO2 approaches 1.0. If an electron acceptor that is lower in the redox ladder than nitrate also exist, then multiple inhibition terms are needed. For example, for iron oxide, we write the following:


Here KI,O2 is the inhibition coefficient of NO3 or the NO3 concentration at which it inhibits iron reduction. The additional nitrate inhibition term means that iron reduction occurs at significant rates only when both oxygen and nitrate are low compared to their corresponding inhibition coefficients.

Rates in natural soils.

The dual-Monod and inhibition term are important under conditions where electron donors and acceptors are limited. In shallow soil, O2 is prevalent except under wet conditions with little pore space for air. Anoxic conditions can develop in local environments such as dead-end pores where water is saturated for a long time and not easily flows out. Under conditions organic carbon and O2 are abundant, the rate law is simplified to the following form assuming microorganism concentrations are relatively constant:


where μmax (mol/L2/t) here depends on the original μmax (mol/L2/t) in previous equations but also on the concentrations of microbe. That is, here we lump the rate constant with microbial abundance. The A is the SOC surface area (L2) as an approximation of SOC content, and f(T) and f(Sw) describe the temperature and soil moisture dependence, respectively. For temperature dependence, a Q10-based form (Friedlingstein et al. 2006; Regnier et al. 2013) is commonly used: f(T)=Q10|T20|/10, where Q10 is the relative increase in reaction rates when temperature increases by 10 °C (Davidson and Janssens 2006). The f(Sw) accounts for the nonlinear dependence of rates on soil moisture. A simple form of f(Sw) = (Sw)ε where ε is the saturation exponent (a typical ε value is 2) is often used. More complex forms of f(Sw) considering both water limitation under dry conditions and O2 limitation under wet conditions have been proposed (Yan et al. 2018). It has also been suggested that SOC decomposition depends strongly on the depth distribution of SOC (Seibert et al. 2009), which is sometime accounted with a depth function:


where Zw is the water table depth (m). An example is f(Zw) = exp(−Zw/bm) (Weiler and McDonnell 2006). Here bm is the declining coefficient describing the gradient of SOC content over depth.

Chemical weathering

Chemical weathering is the process that transforms primary minerals (such as silicates), those formed in the cooling and solidification of molten mass, into secondary minerals such as clays, essentially turning rocks into soils (White and Brantley 1995). In this process, primary minerals dissolve and leach out cations whereas secondary minerals precipitate. At a short time scale (days to years), these reactions change water chemistry. Over long time scales, these reactions alter solid phase mineralogical compositions and physical properties. Soils are the relatively weathered materials so primary minerals are more abundant at depth, the specific distribution of which depend on the position of reaction fronts, or the transition between weathered and unweatherred zones (Brantley et al. 2013b). Still, even the relatively slow-dissolving clay continues to leach out cations when solubilized with infiltrating meteoric water at disequilibrium (Heidari et al. 2017). Chemical weathering has been a rich research topics in geochemistry that has accumulated an impressive literature documenting reaction thermodynamics (Johnson et al. 1992; Wolery 1992) and kinetics, i.e., rate laws that prescribe rate dependence on aqueous chemistry (Blum and Stillings 1995; Oelkers 2001; Brantley 2008). Here I only briefly describe the widely-used Transition State theory (TST) rate law and list a few representative reactions in Table 2. Interested readers are referred to other chapters in this book and earlier RiMG volumes for in depth coverage of the topic.

Silicates and carbonate are among the most common primary rocks on Earth's surface (Moosdorf et al. 2010). Note that in Table 2, carbonate is listed as both primary and secondary minerals, as it can precipitate and form as a secondary mineral. An example of silicate weathering is K-feldspar dissolving out solutes such as Al3+ and SiO2 (Reaction (4) toward the right), which further precipitate and form kaolinite (Reaction (5) toward the left). The net reaction is 2 KAlSi3O8(s) + 2H+ ↔ Al2Si2O5(OH)4(s) + 2 K+ + 4 H2O + 4 SiO2(aq), essentially K-feldspar transforming into kaolinite while leaching out K and SiO2. Note that all reactions except pyrite dissolution and ferrihydride precipitation are reversible, meaning the reactions can go either direction depending on aqueous chemistry and their reaction thermodynamics. The widely accepted Transition State Theory (TST) based rate laws prescribe this reversibility as follows (Lasaga 1998):


Here r is the mineral reaction rate (mol/L3/t), A is the mineral surface area per unit volume (which depends on mineral abundance (L2/L3)), k (kH,kH2O, kOH) are rate constants under acidic, neutral, and alkaline conditions (mol/L2/t), the activities (a) are for H+, water, and OH, n (nH,nH2O, nOH) are the exponents of activities, IAP is the ion activity product, and Keq is the equilibrium constant. The saturation index IAP/Keq quantifies how far the aqueous phase is from equilibrium. At equilibrium, IAP/Keq, equals 1.0. When the saturation index is less than one, the mineral dissolves. Note that in shallow soils with unsaturated water, the surface area A needs to be modified with a dependence on soil moisture, as minerals only dissolve when solubilized in water. An additional term such as f(Sw) = (Sw)ε is necessary to account for that effect, similar to the way we consider soil moisture effects on SOC decomposition. The equilibrium constants determine how much a mineral can dissolve in water (thermodynamics), whereas the rate constants control how fast reactions reach equilibrium.

The exception in Table 2, pyrite dissolution and ferrihydrite precipitation, are redox reactions that involve changes in oxidation states of involved solutes. These reactions cause acid mine drainage and the yellowish water that is enriched with Fe(OH)3 (Druschel et al. 2004). These reactions are irreversible and follow different rate laws (Williamson and Rimstidt 1994; Rimstidt and Vaughan 2003).


A brief history of reactive transport modeling

Multi-component Reactive Transport Models (RTMs) originated in the 1980s and have been extensively used in the subsurface geochemistry community (Chapman 1982; Chapman et al. 1982). RTMs couple flow and transport calculation within a full biogeochemical thermodynamic and kinetic framework (Steefel et al. 2015), therefore explicitly tracing spatial and temporal patterns of geochemical species in fluid and solid phases. Built upon the theoretical framework of reaction thermodynamics and kinetics (Lichtner 1985, 1988), RTM development advanced rapidly in the 1990s as illustrated by the emergence of various RTM codes that have become extensively used in the past decades (Ortoleva et al. 1987; Yeh and Tripathi 1989; Steefel and Lasaga 1994; Bethke 1996; Lichtner et al. 1996; Van Cappellen and Wang 1996; Xu et al. 1999; White and Oostrom 2000; Mayer et al. 2002; Hammond et al. 2014).

RTMs have since been utilized as an integration and interpretation tool across diverse environments involving both porous and fractured media (as reviewed in MacQuarrie and Mayer 2005; Steefel et al. 2005; Sprocati et al. 2019). They have simulated a plethora of processes, including tracer transport, mineral dissolution and precipitation, ion exchange, surface complexation, as well as biotic processes such as microbe-mediated redox reactions, and biomass growth and decay. RTMs have been applied to understand processes on topics including chemical weathering (Bolton et al. 1996; Maher et al. 2009a; Brantley and Lebedeva 2011; Moore et al. 2012; Heidari et al. 2017), biogeochemical cycling (Regnier et al. 1997; Dale et al. 2008; Krumins et al. 2013; Ng et al. 2017), environmentally bioremediation (Li et al. 2010; Yabusaki et al. 2011; Druhan et al. 2012), natural attenuation (Mayer et al. 2001; Liu et al. 2008), geological carbon sequestration (Xu et al. 2003; Atchley et al. 2013; Brunet et al. 2013; Navarre-Sitchler et al. 2013; Tutolo et al. 2015), nuclear waste storage (Saunders and Toran 1995; Soler and Mader 2005), and energy production (Audigane et al. 2007; Qiao et al. 2015).

RTMs have been primarily applied in groundwater and deep subsurface systems at spatial scales from pores (Kang et al. 2006; Li et al. 2008; Fang et al. 2011; Molins et al. 2014; Scheibe et al. 2015; Sprocati et al. 2019) to field scales at tens to hundreds of meters (Li et al. 2011; Beaulieu et al. 2011, 2015; Navarre-Sitchler et al. 2013). Regional scale RTMs have recently been linked to global vegetation models to understand the role of climate change in controlling weathering over periods of 102 to 103 years (Goddéris et al. 2006, 2013; Roelandt et al. 2010). Only recently RTM has been introduced to the watershed scale, where hydrological conditions are transient (Yeh et al. 2006; Bao et al. 2017).

Reactive transport equations at the watershed scale

A watershed reactive transport model needs to couple hydrological processes, land-surface interactions with multi-component reactions to capture the dynamics of water and biogeochemical interactions. The governing equations of reactive transport processes are in fact not that different from those for groundwater systems in that it is determined by advective and dispersive / diffusive transport and reactions. The major difference is that now water flow is dictated by watershed hydrology that is open to changes in meteorological conditions and other watershed-relevant characteristics. The watershed reactive transport code BioRT-Flux-PIHM (Bao et al. 2017; Li et al. 2017a; Zhi et al. 2019) follows the semi-discrete finite volume approach discussed earlier for PIHM, for an arbitrary solute m in an arbitrary prismatic element i, an example reactive transport equation is as follows (Bao et al. 2017):


where Vi is the total volume of the element i; Sw,i is the soil water content (L3 water/L3 porous medium volume) that can be calculated based on water storage from a hydrological model; Cm,i is the aqueous concentration of solute m, mol/L3; j is the index of neighbor elements of i sharing grid interfaces with i; Aij is the grid interface area (L2) shared by i and its neighbor grid j; Dij is the combined dispersion/diffusion coefficient (L2/t) normal to the shared surface Aij; Lij is the distance between the center of i and its neighbor element j; qij is the flow rate across the shared interface Aij, L3/t; Rm is the total rate of kinetically controlled reactions that involve solute m, mol/s; ntot is the total number of solutes that are kinetically controlled.

Note that Equation (14) can be written for any solutes that are of interests with their corresponding rate laws represented in Rm. It describes processes that lead to mass changes of solute m. The first two terms in the right hand side describes the conservative (non-reactive) transport and the last term describes reaction rates. If a solute participates in multiple kinetically-controlled reactions, this rate is the summation of several reaction rates with different rate laws. For example, for a geogenic solute Ca, it can dissolve out of a silicate mineral such as anorthite adding Ca to the water phase but also can precipitate as calcite which is a sink term removing Ca out of the water phase. The Rm then is the net rate of addition and removal following the TST rate law in Equation (13). For DOC or nutrient species coming out of the decomposition of Soil Organic Matter (SOM), the rate follows Equation (12) with parameters specific to reactions in Table 2. Note that DOC is produced by SOC decomposition but can also be consumed in other microbe-mediated reactions as the electron donor in processes such as denitrification. The Rm term in this case would be the net rate of the two microbe-mediated redox reactions (SOC decomposition and denitrification) following their specific rate laws with corresponding rate parameters.

The module Flux-PIHM solves for water storage and fluxes, which are used to drive reactive transport calculation for solute concentrations in the BioRT module. The code is written in a way that users can define solutes, reaction type, and rate laws of interests in the input file, enabling the flexibility based on the need of the users. The code outputs spatial and temporal distributions of solute concentrations and reaction rates of interest. Below I show examples to illustrate how hydrological conditions influence water chemistry at the watershed scale.

Examples of hydological and biogeochemical coupling

How does ET influence water chemistry?

In the Shale Hills example discussed earlier (Fig. 4), although precipitation is relatively similar across the year, summer is dry due to higher ET. These transient hydrological conditions have profound impacts on water chemistry. The non-reactive tracer, chloride (Cl), originates from precipitation, whereas Mg leaches out of clay weathering and participates in ion exchange. In the dry summer, only a small proportion of the catchment, mostly swales and the very vicinity of the stream, is connected to the stream (Fig. 6A). Most Cl is ensnared in isolated pockets of soil water, leading to escalated concentrations that are more than one order of magnitude higher than rainfall concentrations (~2 µmol/L) and those in other locations connected to the stream (Fig. 6B). In wet spring and winter, large hillslope areas are connected to the stream, allowing Cl to export rapidly such that chloride concentration ([Cl]) is relatively low and is spatially homogeneous. In other words, [Cl] varies spatially and temporally because of the water dynamics. When ET is low, [Cl] is relatively low and homogeneously distributed across the entire catchment owing to a well-connected, wet watershed that rapidly flushes out Cl. When ET is high, flow pathways close up, leading to “islands” of elevated [Cl] and a much more heterogeneous concentration field. Across the catchment, [Cl] is high in “old” water but they are almost irrelevant to the stream water because of the disconnection between stream and hillslope.

The Cl mass in the watershed is determined by the balance between rainfall input and discharge output. Low discharge and connectivity in the summer result in negligible Cl export and therefore Cl accumulation within the watershed. Intense rainfall and snowmelt later mobilize trapped Cl, plummeting to ~50% of the total Cl mass, as shown in Figure 7E in Li et al. (2017a). The increase of Cl concentrations with higher ET and lower recharge is well documented and serves as the basis for the Cl mass balance method for estimating groundwater recharge (Rice and Hornberger 1998; Semenov and Zimnik 2015).

Water storage influences Mg concentrations ([Mg]) by modulating the wetted surface area and clay dissolution rates (notably chlorite). Chlorite dissolves faster in swales and valley floor with ample water and connectivity to the stream (Fig. 6C). Averaged dissolution rates at the watershed scale were estimated to drop by about half in the dry summer compared to spring. Notably, [Mg] in soil water does not increase as much as [Cl] in the summer; it also does not vary spatially as much as [Cl], primarily because cation exchange effectively acts as a buffer that dampens concentration fluctuations (Clow and Mast 2010; Herndon et al. 2015). The [Mg] on exchange sites is highest in the valley floor (Fig. 6D) because convergent water flow continuously brings Mg mass fluxes from the upslope.

The DOC derives from SOC (or SOM) decomposition, the rate of which depends on water content, temperature, and depth of water table (Eqn. 13). Its spatial distribution mostly mirrors that of the soil moisture in Figure 6A with relatively high rates in swales, which is similar to chlorite dissolution and suggests the dominant control of water content. They however differ from those of chlorite dissolution patterns because of their dependence on water table (f(Zw)). The shallower water table in swales leads to high f(Zw) values and high SOC rates, amplifying the difference from the planar hillslopes. The DOC concentrations are particularly high in streams in all seasons, because streams receive the produced DOC from the land.

This example illustrates the intimate linkage between water storage, reactions, and solute export. For Cl originating from wet and dry deposition, concentrations in soil water are higher as the watershed dries up at high ET. For the reactive Mg, although clay dissolution is driven by the wetness of mineral surfaces that largely depend on ET, ion exchange reactions buffer its concentration fluctuations, leading to much lower extent of seasonal and spatial variations. For DOC, although it is produced in land that is similar to Mg, its not-as-strong sorption on soils leads to higher spatial and temporal variations in concentrations. In general, however, the variations in solute concentrations are much smaller than the order of magnitude change in discharge.

How do extreme hydrological events influence water chemistry?

Snowmelt and flooding.

Seasonal snow covers 30% of the Earth's land surface, 98% of which is in the Northern Hemisphere, specifically in North America and Eurasia (Robinson et al. 1993; Cheng et al. 2018). In western USA, mountain Snow packs are estimated to feed about three fourth of the freshwater, earning them the name of “water tower” (Cayan 1996). The hydrology of these mountains is snow-dominated, with a distinct snow melting season in late spring to early summer contributing to > 70% of the annual water budget. The significance of snowmelt however does not merely lie in its sheer large water volumes. Snow melting also plays an outsize role in flushing out solutes and contaminants. To illustrate this, Figure 7 shows as an example of the Coal Creek watershed, a representative high elevation watershed (2,700 to 3,700 m, 53 km2) in the central Rocky Mountains of Colorado, with an annual mean temperature around 0.9 °C and average annual rainfall and snowfall of ~ 600 mm and 550 mm, respectively. The watershed has been mined for ~ 100 years in the past with lingering metal-leaching mine tailings at the site, which brings out the eminent water quality issue as Coal Creek supplies drinking water for the skiing town Crested Butte (Zhi et al. 2019).

Figure 7 shows a pronounced discharge peak during snowmelt from early May to late June, 2016. DOC remains at low concentrations under base flow and elevates to a maximum of ~5 mg/L at the discharge peak. In contrast, geogenic cations (Na, Ca, and Mg) exhibit high concentrations under base flow and low concentrations at high flow, demonstrating a dilution trend that is opposite of the DOC pattern. Trace metals (Zn, Mn, and Cd) derived from abandoned mine tailings and natural deposits approach their maxima at early flushes before the discharge peak. Total solute export (mass/d) estimated using the USGS Load Estimator (LOADEST) show that discharge from early May to early July accounts for nearly 80% of the annual discharge (Zhi et al. 2019). In this period, 90% of annual DOC export occurs, compared to 70% and 75% for geogenic solutes and trace metals, respectively. In essence, snow melt is the hot moment of the year not only for water quantity but also for water quality: solute and contaminants have outsized export in this period.

These temporal impacts are in fact similar to flooding, storm water, and melting glaciers (Saberi et al. 2019). Excessive water in large hydrological events has been shown to flush out disproportionally large pulses of “stored” DOC and contaminants that deteriorate water quality and ecosystem health (Raymond and Saiers 2010; Huntington et al. 2016). Rapid changes in discharge during storms and snowmelt can result in shifts in source waters and pronounced variations in DOC concentration, occurring concomitantly with a shift in the optimal properties of DOM that indicates a transition towards higher molecular weight, more aromatic DOM composition (O'Donnell et al. 2010; Voss et al. 2015).


Under changing climate conditions, precipitation experiences shifts from snow to rainfall, accompanied by shrinking snow packs and dwindling summer flow in mountainous watersheds (Berghuijs et al. 2014; Rhoades et al. 2018). Such hydrological shift, coupled with increasing temperature, can fundamentally alter biogeochemical reaction rates and water quality (Smith et al. 2005; Laudon et al. 2012; Milner et al. 2017). The effects of droughts on stream water composition however have remained poorly understood (Larned et al. 2010). Harjung et al. (2018) recently examined the impacts of droughts on DOM using a myriad of drought conditions in stream-side flumes. Their work showed production and release of DOM from biofilms and leaf up by 50% during droughts, accompanied by changes in absorbance and fluorescence properties suggesting more labile DOM for microbial degradation. These observations are consistent with measured DOC concentrations and DOM properties in arid areas with prolonged droughts (Jones et al. 1996; Vazquez et al. 2011). Droughts generate disconnected areas with distinct hot spots, amplifying spatial heterogeneities as observed in Figure 6.


Drivers of chemical weathering in natural systems

As discussed earlier in the section Fundamentals of soil biogeochemistry, chemical weathering and organic matter decomposition follow different rate laws. These rate laws are derived from theories derived considering only chemistry with with parameters measured typically in well-mixed systems without mass transport limitation. Natural systems are almost never well-mixed (Dentz et al. 2011; Loschko et al. 2018). In particular, at the watershed scale, hydrological processes, often shaped by both external forcing and internal structure characteristics, influence reactant concentrations, distance to equilibrium, and redox state. Although large datasets of concentrations and discharges exist and allude to rates at the watershed scale (Bluth and Kump 1994; White and Blum 1995), mechanistic understanding of controls of reaction rates at the watershed scale is generally lacking.

The role of reactive gases and roots.

In shallow subsurface, chemical weathering is influenced not only by the factors outlined in the TST rate law (Eqn. 13), but also by the transient nature of soil moisture and characteristics. The gas of O2 can drive the oxidative weathering of pyrite-containing shale, creating fractures and porosity that allow water infiltration and weathering of other minerals (Moore et al. 2019). Oxidation of pyrite also produces sulfuric acid, which can enhance mineral weathering and create permeable pathways for water and solutes (Brantley et al. 2013a). The presence and abundance of O2 therefore regulate how fast reaction fronts progress as O2-filled water infiltrate through pyrite-rich subsurface (Bolton et al. 2006; Heidari et al. 2017).

Soil CO2 is typically 1–2 orders magnitude higher in vegetated soils compared to the atmospheric CO2 level because of soil and root respiration, therefore driving soil water acidity to a level of pH around 4.0 and accelerating weathering of silicates and carbonates (Hinsinger et al. 2003; Hasenmueller et al. 2015). Organic acids from root exudates also speed up dissolution by increasing soil water acidity and by forming complexity with cations, effectively elevating mineral solubility (Leake et al. 2008; Lawrence et al. 2014). These “biogeochemistry” aspects of root effects have been extensively discussed (Taylor et al. 2009; Brantley et al. 2011; Thorley et al. 2015). Roots of different plants (e.g., grass-, shrub-, and wood-lands) can additionally leave distinct imprints on the relative distribution of soil carbon with depth. For example, the relatively deep root distributions of shrubs can lead to soil carbon profiles deeper than those in grassland (Jackson et al. 1996), which may help maintain high soil CO2 (Oh et al. 2005) and low pH level at depth that is in contact with unweathered, more reactive primary minerals.

Plant roots may also affect soil structure, hydrological redistribution, and different patterns of water flow and soil moisture exchange (Le and Kumar 2017). Macropores originated from root channels have been estimated to account for about 70% of the total macropores (Noguchi et al. 1997). Their structure highly depend on plant species (Beven and Germann 1982, 2013). In grasslands, the lateral, dense spread of roots near surface promotes the development of horizontal macropores and granular or “sandy” texture soil aggregates that facilitate near surface lateral flow (Oades 1993; Cheng et al. 2011). In contrast, forests with deep and thick roots tend to have high proportion of macropores and high connectivity in deep soil (Canadell et al. 1996; Nardini et al. 2016), leading to smaller water-holding capacities and deeper water penetration. In other words, deeper roots in forests can lead to branching, or partitioning of water more to the depth, which can pump CO2-rich water to deeper zones and enhance the contact between CO2 and primary minerals at depth. As an example, woody encroachment into grassland and associated rooting differences have been attributed to elevated chemical weathering fluxes in carbonate terrains (Sullivan et al. 2018). Much remains to be learned about feedbacks among weathering, transport, and biologically-mediated processes with realistic subsurface structures. Watershed RTMs can be useful to distinguish biogeochemical versus hydrological impacts of root structure on chemical weathering.

Dynamic reactive surface area.

The surface area A in Equations (12) and (13) is essential in determining rates and is arguably the most uncertain parameter in the rate laws. For SOC decomposition under aerobic conditions, this surface area quantifies the contact between reactants, i.e., microbe, water, SOC, and O2. For mineral dissolution, this A quantifies the effective mineral surface that is actually reacting. Five to six orders of magnitude discrepancy between rates of mineral reactions measured in laboratory versus field context (Blum and Stillings 1995; Drever and Clow 1995; Sverdrup et al. 1995; White 1995; White and Brantley 1995, 2003; Richards and Kump 2003). This rate discrepancy has been attributed to two major categories of factors. One category includes processes and conditions that directly affect the “reactivity” of mineral surfaces. Examples include mineral surface roughness or fractalness (Navarre-Sitchler and Brantley 2007), passivation or armoring layers as a result of clay coating or secondary mineral precipitation (Helgeson 1971; Nugent et al. 1998; Maher et al. 2009b), leached layers that act as diffusion barriers (Luce et al. 1972), and water films between grains that vary under different temperature and pressure conditions (Renard et al. 1997). The second category includes factors that regulate fluid compositions. For example, water flow controls reaction rates by regulating the extent of deviation from equilibrium (Maher 2010; Salehikhoo et al. 2013). The amount of wetted surface under unsaturated conditions, as well as concentrations of catalyzing or inhibiting aqueous species, also play a significant role (Oelkers et al. 1994; Lawrence et al. 2014).

In addition, porous medium characteristics—the details of pore structure and distribution of porous medium properties (e.g., permeability and mineral abundance)—can lead to large spatial variations in flow rates and deviations from equilibrium, therefore modifying the overall “effective” reaction rates (Kang et al. 2006; Li et al. 2007; Molins et al. 2012; Al-Khulaifi et al. 2017; Jung and Navarre-Sitchler 2018) or even change the reaction direction (Brunet et al. 2016). For example, a mineral that does not have any coating can be bathed in a fluid at equilibrium in the low permeability zone of a heterogeneous medium and therefore is not effectively dissolving. In fact, any factors that influence reaction kinetics (surface area, concentrations of catalyzing or inhibiting species, and distance from equilibrium) can contribute to the rate discrepancy. In natural systems, all these factors operate simultaneously resulting in much smaller effectively-dissolving surface areas than the measured Brunauer–Emmett–Teller (BET) or imaged surface area, although the dominance of different factors may vary with conditions; e.g., (White and Peterson 1990; Nugent et al. 1998; Zhu 2005; Jin et al. 2011a).

At the watershed scale, the “effective” surface area is influenced not only by the subsurface structure, not also by how much water in different parts of the watershed is connected to the stream. As shown in Figure 6, although minerals are dissolving everywhere, not all areas are connected to the stream and have active flow that drives to disequilibrium. The hydrological regime has a tremendous impact on weathering rates (Torres et al. 2015, 2017), which are controlled via effective surface area at the watershed scale. How are they relate to reaction rates at the watershed scale? Watershed RTM can be used to answer these questions.

Watershed response to hydrological changes: Concentration discharge relationships

Although the integration of reactive transport and watershed hydrology in process-based models is nascent, field measurements have in fact documented decades of stream flow and water chemistry data worldwide. Many rivers in the USA, for example, the Mississippi River and many of its tributaries, have water quality data dating back to 1950s (https://waterdata.usgs.gov/nwis?). A direct measure of biogeochemical responses to hydrological changes is the concentration (C) and discharge (Q) relationships measured at river and stream mouths. CQ relationships have been extensively used to understand watershed structure and response to external forcing (Winterdahl et al. 2016) and to quantify fluxes of chemical weathering (White and Blum 1995) and nutrient export (Zhang et al. 2016a).

Studies of CQ date back more than five decades on event-based CQ hysteresis (Johnson et al. 1969; Hooper et al. 1990; Evans and Davies 1998). The power law relationship C = aQb remains the extensively used form (Godsey et al. 2009) to characterize CQ patterns, although recent work has proposed more complex characterization (Thompson et al. 2011; Moatar et al. 2017; Zhang and Ball 2017). Chemostasis is defined as relatively small variations in concentrations compared to discharge with b values within −0.1 and 0.1 (Musolff et al. 2015). Chemodynamic patterns are characterized as having high absolute b values (> 0.1), with positive b values demonstrating flushing or enrichment behavior whereas negative b illustrating dilution.

A big puzzle in this area is the contrasting CQ relationships that have been observed for different solutes: similar CQ patterns have been observed for different solutes; contrasting CQ relationships have been observed for the same solute in different watersheds (Herndon et al. 2015; Musolff et al. 2017; Abbott et al. 2018; Zarnetske et al. 2018). Geogenic species (e.g., Na, Mg, Si) derived from chemical weathering have commonly demonstrated chemostasis (Godsey et al. 2009; Sullivan et al. 2018). Flushing behaviors are most commonly observed for DOC (Boyer et al. 1997; Zarnetske et al. 2018). For nutrients, chemostatic or biogeochemical stationary behaviors have been commonly observed in agricultural lands due to their large legacy store from decades of fertilizer applications (Basu et al. 2010; Van Meter et al. 2017) whereas dilution behavior occurs under conditions where groundwater is the major source of nutrient export (Miller et al. 2016). An example is shown in Figure 8 for the Coal Creek watershed discussed earlier. Here biogenic solutes (DOC, TN, and TON) exhibit flushing behavior, geogenic species manifest dilution behavior, and all metals show hysteresis loops.

Several schools of thoughts in the literature explain contrasting CQ patterns. Chemostasis has been attributed to hydrologically responsive dissolution rates (via wet surface area) as discharge increases (Li et al. 2017a), to buffering capacity of ion exchange reactions (Clow and Mast 2010) (Clow and Mast 2010), and to the approach to equilibrium as residence time increases (Maher 2011; Ameli et al. 2017). Flushing behaviors have been interpreted as due to rising water table tapping shallow soil zones with enriched organic carbon and nutrients, increased hydrological connectivity between uplands and stream under high flow conditions, and biotic control of vegetation and temperature (Seibert et al. 2009; Andrews et al. 2011; Mei et al. 2012). In contrast, dilution patterns have been attributed to high solute concentrations from groundwater diluted by rain water at high discharge (Miller et al. 2016; Li et al. 2017a), flow paths bypassing soil pore water with high concentration, and rapid depletion of solute store (Herndon et al. 2015b; Hoagland et al. 2017). A unifying framework has yet to emerge to interpret the diverse observations.

Using bioRT-flux-PIHM and Monte-Carlo simulations using a wide range of conditions, Zhi et al. (2019) illustrates that contrasting CQ patterns can be explained by the switch of source waters under dry and wet conditions, and the source water chemistry contrasts arising from spatial heterogeneity of source materials in the subsurface. In particular, stream water mixes different types of source water, primarily shallow soil water (e.g., perched water table), deep (relatively) groundwater, or water from riparian zones and hillslope areas. The relative contribution of different source waters vary drastically under dry and wet conditions. As shown Figure 6, under wet conditions, stream water is dominated by shallow soil water in largely connected watersheds; under dry conditions, groundwater or water in the very vicinity of streams predominate. The chemistry of these waters differ dramatically. Shallow soils are characterized by rich life and high SOC and DOC content, whereas deeper subsurface groundwater brings high concentrations of geogenic species at depths.

Solutes that are enriched in shallow soils (e.g., DOC) therefore have flushing patterns because their concentrations increase as water table rises under wet conditions and taps shallow soils; in contrast, those abundant at depths show up at high concentrations in base flow during dry times, demonstrating dilution patters. These insights have been encapsulated in a general relationship between CQ slope b and ratios of soil versus ground water concentrations as the two end-member source waters (Cratio = Csw/Cgw): b=δbCratioCratio,1/2+Cratio+bmin (Zhi et al. 2019), where bmin is the minimum b value, δb is the difference between maximum and minimum b values, and Cratio,1/2 is the concentration ratio when the b values equal the average of minimum and maximum b values. Disparate CQ patterns of 11 solutes (DOC, DP, NO3, K, Si, Ca, Mg, Na, Al, Mn, and Fe) from three watersheds (i.e., Coal Creek, Shale Hills, and Plynlimon) with different climate and geology conditions follow this general relationship, indicating its promise of quantifying CQ patterns for a wide range of solutes in watersheds of diverse characteristics. This underscores the essential role of shallow–deep biogeochemistry contrasts in shaping CQ patterns. This comes along similar line to another study that contributes distinct CQ patterns as rising from spatial structure heterogeneity of source waters and travel time that allows transformation reactions to occur within a watershed (Musolff et al. 2017). Similarly, spatial distributions in the horizontal direction or the proximity of solute source zones controls their connectivity to the stream under different flow regimes and therefore also regulates stream chemistry (Tunaley et al. 2017). The general relationship will need to be tested in diverse watersheds with additional data on soil and groundwater chemistry and subsurface structure, and with watershed RTMs that can distinguish the effects of individual processes and properites.

Reaction rates at the watershed scale: linking travel time, age, and reactions

In homogeneous groundwater systems, chemical weathering rates and geogenic concentrations have been shown to hinge on dimensionless numbers such as the Damkohler number (Maher 2011; Salehikhoo et al. 2013; Salehikhoo and Li 2015), defined as the relative time scales of reactions (to reach equilibrium, τeq) versus advective transport (i.e., residence time τr). If the time it takes to reach equilibrium is longer than the time that water stays in a system, i.e., τeq > τr, solute concentrations leaving the system depend on flow rates; if τeq < τr, solutes reach and remain at equilibrium concentrations at the system outlet and do not vary with flow rates. In heterogeneous systems where fast-dissolving minerals (e.g., carbonate) reside in low-permeability zones and slow-dissolving minerals (e.g., quartz) in high permeability zones, solute concentrations and weathering rates have been shown to depend similarly on Damkohler number but need to additionally consider the relative magnitude (ratio) of the contact time between water and reacting minerals (τad,r) versus the overall residence time in the entire domain (τa) (Wen and Li 2018). This is because only a fraction of water flow, not the total flow, channels through the low-permeability reactive zones and dissolves reacting minerals (Wen and Li 2017). Hence the dissolution rates not only depend on dissolution kinetics and surface area, but also on overall flow rates and characteristics of spatial distribution and permeability contrasts between reactive and non-reactive zones. These studies accentuate the key linkage between water contact time and dissolution rates.

Similar framework has yet to be developed at the watershed scale. How do biogeochemical rates at the watershed scale relate to different time scales? How do they differ from groundwater systems with steady-state flow? The catchment hydrology community has been grappling with the question “how old is stream water” using water isotopes and non-reactive tracers (e.g., Cl) for decades and have gleaned ample knowledge on travel time at the watershed scale. Transit (or travel) time, defined as the difference in the time rain water enters a catchment to the time it exits, has been studied extensively as “a fundamental catchment descriptor that reveals information about storage, flow pathways and source of water in a single characteristic.” (McGuire et al. 2005; McGuire and McDonnell 2006). Measurements and modeling of stable water isotopes and tracer (e.g., Cl) have estimated mean and distribution of transit times to understand where stream water comes from and how old it is under diverse climate, topography, soil property conditions (McGuire et al. 2005; Soulsby et al. 2006; Tetzlaff et al. 2009; Godsey et al. 2010; Birkel et al. 2011). The field has also evolved from assuming temporal stationary (McGuire and McDonnell 2006), to the recent development of time-variant transit time distribution (TTD) theory with StorAge Selection function (SAS), recognizing the impacts of non-stationary of (or, temporally variable) rainfall events and water storage (Botter et al. 2011; van der Velde et al. 2012; Harman 2015; Rinaldo et al. 2015). There has also been considerable advocate on using the young water fraction to characterize catchment travel time (Jasechko et al. 2016; Kirchner 2016a). In addition to travel time, other concepts have emerged as key measures of watershed hydrology. In particular, hydrological connectivity, defined as the extent of connection from hillslope uplands to streams (McGlynn and McDonnell 2003a, b), has been shown to govern water flow and travel times (Western et al. 2001; Stieglitz et al. 2003; McGuire and McDonnell 2006; Jencso et al. 2009, 2010; Gooseff et al. 2017), and can determine the extent and rates of reactions.

These concepts and quantifications have yet to be linked to biogeochemical measures at the watershed scale. As discussed earlier, concentration discharge measures abound, from which the solute export rates and potentially biogeochemical reaction rates can be inferred. These rates and concentrations have started to be linked to flow and travel time at the watershed scale (Benettin et al. 2015b; Ameli et al. 2017; Musolff et al. 2017). As an example, the use of general age theory for reactive solutes (N species: ammonia and nitrate) and comparison to a non-reactive solute in agricultural lands have found that the age of N is in fact higher in top soil and in low depression area because of the ion exchange reaction that retards N and the tile drains that mobilize N from immobile zones (Woo and Kumar 2016, 2019). In addition, biogeochemistry reactive interfaces (Krause et al. 2017; Li et al. 2017b), defined as where water flows and converge from contrasting geological, ecohydrological or biogeochemical compartments, are “hot spots” (McClain et al. 2003) where rates are often orders of magnitude higher than the rest of a domain (Krause et al. 2013; Harvey and Gooseff 2015; Salehikhoo and Li 2015). Example include surface water–ground water mixing zones, hyporheic zones at river–land-sediment interfaces, spring and seep emergence zones, subsurface high-low permeability boundaries, and riparian zones at land-river boundaries (Öquist et al. 2009; Stonedahl et al. 2010; Hefting et al. 2013; Gooseff et al. 2016; Dwivedi et al. 2018). With the capability of explicitly calculating these quantities, watershed RTMs can help unravel key linkages between these concepts and control of biogeochemical transformation at the watershed scale.


This chapter focuses on the fundamentals of watershed hydrology and biogeochemistry and the recent development of watershed RTMs. As has been shown here, biogeochemical reactive transport is water-driven and depends strongly on the dynamics of water fluxes including ET, and different source waters flowing to streams under different hydrological regimes. Unlike groundwater systems where flow is spatially heterogeneous and at steady state, watershed hydrological processes are spatially heterogeneous and temporally dynamic. This adds additional layers of complexity in understanding process coupling and feedbacks at the watershed scale.

Model development in watershed hydrology and biogeochemistry have leaped forward without much interactions until very recently. Observations that record water and water chemistry however have accumulated for decades owing to the long-term water monitoring programs worldwide. In fact, we are now at an exciting time when a deluge of data exist. The past decades have witnessed rapid advances in technology and unprecedented generation of earth surface data from remote sensing from satellites (Hrachowitz et al. 2013; McCabe et al. 2017). Within the U.S. alone, such programs include the Critical Zone Observatories (CZOs), the Long-Term Ecological Research (LTER), the Great Lake Ecological Observatory Network (GLEON), the National Ecological Observatory Network (NEON), and the testbed research sites supported by the Department of Energy Subsurface Biogeochemical Research program (e.g., East River Watershed (Hubbard et al. 2018)). Internationally a myriad of observatories exist (Bogena et al. 2018; Gaillardet et al. 2018). Water and water chemistry data are collected regularly through research and observation networks, presenting unprecedented opportunities for integrated understanding (Baatz et al. 2018). There is also a growing recognition that more integration is needed as we move forward (Grimm et al. 2003; Hrachowitz et al. 2016).

These rich data and process-based modeling tool will enable us to answer questions about linkages between travel time, water age, and biogeochemical rates, and interrogate their functional relationships and underlying mechanisms. Such efforts can guide field campaigns: instead of measuring everything everywhere, we may focus on end-members in distinct zones (surface water, soil water, and groundwater) that dominate different hydrological regime; we may need to focus on hot spots (e.g., riparian zones) and hot moments (e.g., storm events, snow melt, droughts) where biogeochemical transformation rates are disproportionally high. It is also important to go beyond individual watersheds, as each watershed has its own hydroclimatic conditions and idiosyncrasies (e.g., land cover, soil type, lithology, and topography) but general principles often underlie the apparent contrasts and offer common threads. Synthesis studies that explore a large number of watersheds are needed to answers questions including: how the linkages and functional relationships developed in individual watersheds vary under diverse climate, land cover, and geology conditions? What are the first order control and watershed characteristics that can be used for prediction? Answers to these questions will ultimately mount to general patterns, principles, and classification. Watershed RTMs can help guide the field campaigns but also set up digital watersheds and interrogate the role of external forcings and internal watershed idiosyncrasies.

Cross-site understanding can also facilitate understanding cross-scale behavior, i.e., the connection between small-scale physics and large-scale behavior (scaling). The issue of scaling is a major challenge and has been extensively discussed in hydrology and other disciplines (Beven 1989; Levin 1992; Sivapalan et al. 2003; Hrachowitz et al. 2013). Reactive transport models at the pore scale have run at the resolution of microns using synchrotron X-ray tomographic images and direct numerical simulation (Wen et al. 2016; Molins et al. 2019). At the watershed scale, models have been run at a spatial resolution as high as one meter using Light Detection and Ranging (LIDAR) images (Le and Kumar 2017; Woo and Kumar 2017). These models illuminate the important role of small-scale spatial characteristics such as microtopography, depression area, and ponding, in regulating chemical weathering and the age of nitrogen in agriculture soils. The issue of upscaling however almost always arises when prediction is needed at larger spatial scales. The resolution of the National Water Model is approximately ~ 1 km2; the resolution of Earth system Models is often considerably coarser. An individual watershed may be one homogeneous grid block or even smaller than one grid block in large scale models. To represent watershed reactive transport in these models, we need to go beyond spatial heterogeneities and process complexity that obscure effects of individual processes. Such upscaled understanding can not only reduce computational cost and facilitate explicit incorporation of feedback schemes, but also hold the key to unlock a unifying hydro-biogeochemical theory. Ultimately, going up and down scales are complementary and depend on the needs of answering target research questions.


The author acknowledges the support from the NSF EAR – 1331726 (S. Brantley) for the Susquehanna Shale Hills Critical Zone Observatory (SSHCZO) and the Department of Energy Subsurface Biogeochemical Research program DE-SC0016221 (L. Li) for the development of bioRT-flux-PIHM and the initiation of watershed reactive transport modeling. The author acknowledges collaborators and students heavily involved in the code development, including Chris Duffy, Chen Bao, Wei Zhi, and Yuning Shi. The author appreciates the hospitality and stimulation of Andrea Rinaldo, Paolo Benettin, and other ECHO laboratory members at the Ecole Polytechnique Fédérale de Lausanne (EPFL) while developing this chapter during sabbatical. Comments from the Editors Jennifer Druhan and Christopher Tournassat and an anonymous reviewer have improved the chapter.

Open access: Article available to all readers online.