Supplementary Methods
Watershed reactive transport modeling. BioRT-Flux-PIHM (BFP)
augments the Flux-PIHM family code and is a recently developed reactive
transport module BioRT for modeling watershed hydrological and
biogeochemical processes. The model couples three modules: a
multicomponent reactive transport module BioRT, the land‐surface
interaction module Flux for processes such as solar radiation and
evapotranspiration (ET), and the surface hydrology module PIHM for
hydrological processes (e.g., precipitation, infiltration, recharge,
surface runoff, and subsurface flow). The model simulates both various
types of abiotic and biotic reactions with different reaction kinetics
and thermodynamics, including mineral dissolution and precipitation, and
microbe-mediated redox reactions, ion exchange, surface complexation,
and aqueous complexation. Detailed features are documented in other
publications (Bao et al., 2017; Li, 2019; Li et al., 2017; Zhi et al.,
2019).
Conewago Creek watershed, a tributary to the Susquehanna River in the
Chesapeake Bay, was used as a model Watershed. It drains an area of 136
km2 and is low in elevation (85 – 350 m). Within a
total of 47% agricultural land, many of the main stem and tributary
floodplains are actively pastured (33%) or cultivated for crop
production (14%). The rest of the watershed is mostly forested (43%)
and developed (17%). Mean annual precipitation and temperature are
1,100 mm and 10.6 °C, respectively. The USGS has monitored discharge and
water-quality since June 2011 near Falmouth, Pennsylvania (station ID:
01573710). Discharge record is available at the daily frequency and
water quality has been measured bi-weekly.
Model setup and calibration (base case). The model was set up
in a spatially implicit manner with watershed characteristics including
initial water and geochemistry conditions, average soil properties,
average land cover, and climate forcing (time series of temperature,
wind speed, solar radiation, precipitation rate and chemistry) based on
data from National Elevation Dataset, National Land Cover Database,
Moderate Resolution Imaging Spectroradiometer (Leaf Area Index), and
Soil Survey Geographic Database. The leaching of nitrate from various
sources was parsimoniously represented as leaching from a generic soil
solid N (possible fertilizers and Soil organic Matter) with the
following reaction and rate law:\(\text{Soil\ N\ }\left(e.g.\ fertilizer,\ SOM\right)\mathbf{\rightarrow}\ \text{NO}_{3}^{-}\ ,\ \)whereSOM is soil organic matter (e.g., plant residuals, biomass).
Its rate was represented as\(\text{kAf}\left(\text{Zw}\right)f\left(T\right)f\left(\text{Sw}\right),\ \)where r is the reaction rate (mol/m3/s); \(k\)is the rate constant (10-9,
mol/m2/s); and A is the effective surface area for
nitrate leaching (m2), which is calculated based
specific surface area (SSA, m2/g), volume fraction
(0.02 m3/m3), and solid density
(urea, 1.34 g/cm3). The depth function\(f\left(Z_{w}\right)=\exp\left(-\frac{Z_{w}}{b_{m}}\right)\)accounts for decreasing soil N concentration along with depth (Van Meter
et al., 2016), where \(Z_{w}\) is the water table depth (m) and\(b_{m}\) (0.5 m) is a declining coefficient (Hagedorn et al., 2001;
Weiler & McDonnell, 2006). The \(f\left(T\right)\) and \(f(S_{w})\)describe the rate dependence on soil temperature and moisture,
respectively (Zhi et al., 2019). The \(f(T)\) takes the widely used
Q10-based form\(f\left(T\right)=\ Q_{10}^{\left|T-20\right|/10},\ \) where\(T\) is the soil temperature (°C) and \(Q_{10}=2\) is the relative
increase in reaction rates when temperature increases by 10 °C (Hararuk
et al., 2015; Zhou et al., 2009). The \(f(S_{w})\) takes the form of\(f\left(S_{w}\right)=\ (S_{w})^{n},\) where \(n=2\) is the
saturation exponent reflecting effects of soil water content on reaction
(Hamamoto et al., 2010; Yan et al., 2016).
The model was calibrated based on the best-fit criteria for the Nash -
Sutcliffe efficiency and by reproducing time series of discharge and C-Q
pattern of stream nitrate. The calibrated hydrology parameters are in
Table S1 (with most sensitive ones in bold). The porosity \(\theta\)defines subsurface water storage with a larger water storage resulting
in a larger baseflow and wide but low discharge peaks. The saturated
hydraulic conductivity \(K_{\text{satV}}\) and \(K_{\text{satH}}\) are
key parameters for recharge of infiltrated water and lateral flow to
stream, influencing the magnitude and timing of discharge peaks. The van
Genuchten parameters α and n primarily control soil water
retention and therefore have strong impacts on discharge peaks (Shi et
al., 2014). The biogeochemical calibration adjusts SSA for soil nitrate
leaching reaction. The calibrated SSA is 2.2 m2/g that
is similar to the reported range of 1.7 to 3.5 m2/g
(Eghbali Babadi et al., 2019). Groundwater flow was estimated based on
iterative calibration between discharge and stream chemistry. The model
and data comparison is in Figure S6, where the calibrated model
generally reproduced the dynamics of baseflow and discharge peaks and
stream nitrate.
Monte-Carlo Simulations. Monte-Carlo simulation (500 cases)
were carried out to cast the base model results to broader conditions
representing the four land uses with varying subsurface N-containing
material abundance. The conditions for the 500 simulations were randomly
sampled across the uncertainty range of soil N fraction (0.01% to 1%)
(Jurgensen et al., 2017; Shepherd et al., 2015) and groundwater nitrate
concentration (Cdw, 0.1 to 10 mg/L, Figure 3a and Figure
S4) while keeping all other hydrological and N leaching kinetic
parameters the same as in the base case. As the volume fraction of soil
N increases, the soil water concentration (Csw) also
increase because more NO3- is leached
out. These ranges covered the different land use types represented with
their corresponding soil and groundwater concentration ranges (Figure
S8). For each run, the C-Q slope \(b\) was calculated based on the daily
model output of concentration and discharge data. In parallel, annual
averages of soil water and groundwater concentrations were used to
calculate Cratio used to plot the b versus
Cratio in Figure 6.
In addition, four hypothetical cases for Agriculture, Mixed,
Undeveloped, and Urban were simulated based on available
Csw and Cgw to illustrate typical
nitrate export patterns under different land use conditions (Extended
Data Figure 7). Similar to the conceptual diagram (Figure 1), results
show that agriculture-dominated land uses (i.e., Agriculture and Mixed)
elevated stream nitrate levels and exhibited flushing C-Q patterns while
the Urban land showed a dilution pattern. Undeveloped, pristine land was
lowest in nitrate concentration yet exhibited a flushing pattern due to
higher Cratio (= Csw /
Cdw).