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