\documentclass[10pt]{article}

\usepackage{fullpage}
\usepackage{setspace}
\usepackage{parskip}
\usepackage{titlesec}
\usepackage[section]{placeins}
\usepackage{xcolor}
\usepackage{breakcites}
\usepackage{lineno}
\usepackage{hyphenat}





\PassOptionsToPackage{hyphens}{url}
\usepackage[colorlinks = true,
            linkcolor = blue,
            urlcolor  = blue,
            citecolor = blue,
            anchorcolor = blue]{hyperref}
\usepackage{etoolbox}
\makeatletter
\patchcmd\@combinedblfloats{\box\@outputbox}{\unvbox\@outputbox}{}{%
  \errmessage{\noexpand\@combinedblfloats could not be patched}%
}%
\makeatother


\usepackage{natbib}




\renewenvironment{abstract}
  {{\bfseries\noindent{\abstractname}\par\nobreak}\footnotesize}
  {\bigskip}

\titlespacing{\section}{0pt}{*3}{*1}
\titlespacing{\subsection}{0pt}{*2}{*0.5}
\titlespacing{\subsubsection}{0pt}{*1.5}{0pt}


\usepackage{authblk}


\usepackage{graphicx}
\usepackage[space]{grffile}
\usepackage{latexsym}
\usepackage{textcomp}
\usepackage{longtable}
\usepackage{tabulary}
\usepackage{booktabs,array,multirow}
\usepackage{amsfonts,amsmath,amssymb}
\providecommand\citet{\cite}
\providecommand\citep{\cite}
\providecommand\citealt{\cite}
% You can conditionalize code for latexml or normal latex using this.
\newif\iflatexml\latexmlfalse
\providecommand{\tightlist}{\setlength{\itemsep}{0pt}\setlength{\parskip}{0pt}}%

\AtBeginDocument{\DeclareGraphicsExtensions{.pdf,.PDF,.eps,.EPS,.png,.PNG,.tif,.TIF,.jpg,.JPG,.jpeg,.JPEG}}

\usepackage[utf8]{inputenc}
\usepackage[ngerman,greek,english]{babel}



\usepackage{float}








\begin{document}

\title{The shallow and deep hypothesis: vertical chemical contrasts shape
nitrate export patterns from different land uses}



\author[1]{Wei Zhi}%
\author[1]{Li Li}%
\affil[1]{Pennsylvania State University}%


\vspace{-1em}



  
  \date{May 07, 2024}


\begingroup
\let\center\flushleft
\let\endcenter\endflushleft
\maketitle
\endgroup





\selectlanguage{english}
\begin{abstract}
Excessive nutrients have penetrated water bodies worldwide yet our
forecasting capabilities remain elusive. This work tests the shallow and
deep hypothesis: chemical contrasts in the subsurface shape nitrate
export patterns. We use data synthesis for 228 U.S. watersheds and
reactive transport modeling (500 simulations) spanning broad climate,
geology, and land conditions. Data synthesis showed that human
perturbation has amplified chemical contrasts in shallow versus deep
waters, inducing primarily flushing patterns (concentration increase
with streamflow) in agriculture lands and dilution (concentration
decrease with streamflow) patterns in urban watersheds. Data and model
reveal a general relationship between export patterns and shallow versus
deep nitrate contrasts. This underscores the often-overlooked role of N
distribution over depth. The results challenge commonly-held perception
that legacy stores in agricultural lands induce chemostasis with
negligible concentration variations with discharge. They suggest
nutrient export will exacerbate as extreme events such as flooding
intensify in the future climate.%
\end{abstract}%



\sloppy


\textbf{INTRODUCTION}

Eutrophication in rivers, lakes, and coastal areas worldwide has
persisted for decades due to intensive agriculture and urbanization
(Diaz \& Rosenberg, 2008; Galloway et al., 2004; Le Moal et al., 2019;
Royer et al., 2006; Van Cappellen \& Maavara, 2016; Van Meter et al.,
2017). In the U.S. alone, over 60\% of U.S. estuaries and coastal water
bodies have been degraded by excessive nutrient inputs (Howarth et al.,
2002). The concentration levels and export patterns of nitrate, a major
dissolved nitrogen (N) species, are essential to understand and quantify
N export from terrestrial to aquatic ecosystems. Export patterns are
often quantified using concentration - discharge (C-Q) relationships
with a power law of C = \emph{a} Q\emph{\textsuperscript{b}} . The C-Q
relationships determine loads of solute export (load = C x Q) and
reflect the response of earth systems to hydrological changes. The value
of \emph{b} indicates different export patterns: a positive, high value
means a flushing pattern where concentrations increase with streamflow,
whereas a negative \emph{b} indicates decreasing concentrations with
streamflow (Godsey et al., 2009; Moatar et al., 2017; Musolff et al.,
2015). Chemostatic patterns (with absolute b values close to zero) occur
when concentrations vary little compared to hydrological variations.
These export patterns have important implications on how and how much
nitrate exports across different hydrological regimes. High
\emph{b}values with pronounced flushing patterns indicate highly
sensitive, escalated export during extreme hydrological events such as
flooding.

Agriculture and urbanization have significantly modified subsurface
physical and chemical structures, including the spatial distribution of
N (Figure 1). Agricultural lands are often characterized by shallow flow
via tile drainage and heavy fertilization that enriches N in shallow
soils (Van Meter et al., 2016; Woo \& Kumar, 2019). Urban watersheds
represent a different type of human perturbation, often characterized by
impervious surfaces that facilitate surface runoff during storms and
sewer pipes that enhance rapid shallow subsurface flow (Grimm et al.,
2008). Nitrogen in urban watersheds can come from a variety of surface
and underground sources including atmospheric deposition, lawn
fertilizers, and domestic and industrial wastewater (Baker et al., 2001;
Divers et al., 2014). Buried leaky sewer systems often elevate the N
input to urban groundwater (Groffman et al., 2003; Pickett et al.,
2011). These manipulations do not occur in undeveloped, pristine
watersheds, where N content is often low and tightly cycled (Weitzman \&
Kaye, 2018). These differences, in particular their vertical
distribution of N content, can lead to distinct C-Q patterns across
different land use conditions.

Given the tremendous human influence, an open question that is key to
understand nutrient export is:\emph{How do subsurface physical and
chemical structures from different land uses in governing nitrate export
patterns?} Contrasting observations in literature do not converge to a
unifying framework that answers these questions. In agricultural lands,
both dilution and flushing have been observed (Jiang et al., 2010;
Miller et al., 2017). The common perception is that chemostasis or
biogeochemical stationarity prevails in agriculture lands, primarily due
to the large legacy store of nitrogen (N) that induces transport
limitation and buffers concentration variability (Basu et al., 2010;
Basu et al., 2011; Thompson et al., 2011). In undeveloped, forest
watersheds, nitrate concentrations have been shown to often peak during
spring floods (Creed \& Band, 1998; Creed et al., 1996; Pellerin et al.,
2012; Sickman et al., 2003), whereas chemostasis has not been as
commonly observed (Huang et al., 2012; Jacobs et al., 2018). Although
much less examined, existing urban studies have observed both dilution
(Duncan et al., 2017) and flushing patterns (Poor \& McDonnell, 2007;
Wollheim et al., 2005).

Growing literature has recognized the importance of subsurface
structure, particularly the hydrological and biogeochemical contrasts in
shallow and deeper zones in governing nitrate levels, denitrification
potentials (Bishop et al., 2004; Creed et al., 1996; Kolbe et al.,
2019), and solute export in general (Musolff et al., 2017; Zhi et al.,
2019). Seibert et al. (2009) developed a flow-concentration integration
model that quantitatively links variations of stream chemistry to flow
rates and vertical profile of shallow soil water chemistry. Zhi et al.
(2019) showed that vertical chemical contrasts between shallow soil
water and deeper groundwater can explain diverse export patterns in
multiple watersheds under a gradient of climate, geology, and land cover
conditions. Here we propose the \emph{shallow versus deep hypothesis:
the N concentration contrasts in shallow and deep zones determine
nitrate export patterns under different land use conditions} . In
particular, we hypothesize that nitrate concentration contrasts in
shallow water (C\textsubscript{sw}) and deep water (C\textsubscript{dw})
differ in different land uses, and that the magnitude of concentration
contrasts governs export patterns. Here ``shallow water'' and ``deep
water'' are broadly defined as waters from different depths contributing
to streams (Figure 1). The shallow water can come from the shallow
surface / subsurface zone that dominates streamflow at high flow; the
deep water may derive from deeper subsurface zones and predominate at
low flow. These shallow and deep waters are loosely defined. In pristine
sites, the shallow water may be the shallow soil water; in agriculture
sites, it can be the combination of surface runoff, tile drainage, and
soil water; in urban lands, it may be the runoff from impervious
surfaces and shallow subsurface pipes. The deep water is the groundwater
that interacts with streams and carries chemical signatures that are
distinctive from shallow water.

If this hypothesis is true, we expect higher concentrations in shallow
water in agricultural lands, leaning more toward a flushing pattern. In
contrast, in urban watersheds, concentrated nutrients accumulated in
leaky pipes in deeper subsurface are often higher than shallow, rapid
runoff on impervious surfaces, leading more toward a dilution pattern
(Duncan et al., 2017). Nitrogen in forests and pristine sites can come
from decomposition from organic matter in shallow soils and leaching
from N-containing rocks (Mayer et al., 2002; Montross et al., 2013). In
addition, they are often tightly cycled. These characteristics can lead
to varied patterns. We test this hypothesis combining a national-scale
synthesis for 228 sites in the United States data synthesis and
mechanism-based reactive transport modeling. Our goal is to propose a
general, simple conceptual model that can explain and quantify export
patterns under diverse conditions.

\textbf{APPROACHES AND METHODS}

\textbf{Water data source.} Nitrate data from 228 sites across the
contiguous United States were analyzed. This includes 147 sites from the
National Water-Quality Assessment (NAWQA) project, a network established
to monitor stream water quality since 1991 (Spahr et al., 2010). The
NAWQA has daily discharge records and monthly water quality data across
four land use categories: Agriculture, Mixed, Undeveloped, and Urban.
The NAWQA sites have median and mean drainage areas of 150 and 294
km\textsuperscript{2}, respectively. The USGS Water Quality Watch (WQW)
database (https://waterwatch.usgs.gov/wqwatch/) includes 71 sites and
continuous water quality data (5 minutes to hourly) since 2014. The WQW
sites are primarily located in the Midwestern ``Corn Belt'' (e.g., Iowa,
Illinois, Indiana) with sizes ranging from 1 to \textgreater{} 29,000
km\textsuperscript{2}. The USGS discharge and nitrate data were
retrieved from the National Water Information System
(https://waterdata.usgs.gov/nwis).

We also examined intensively-measured data in 6 undeveloped and 4 urban
sites from research network sites, including Critical Zone Observatories
(CZOs, http://criticalzone.org/national/), Long Term Ecological Research
Network (LTER, https://lternet.edu/), and Watershed Function Scientific
Focus Area (SFA, https://watershed.lbl.gov/). The 6 undeveloped sites
include Shale Hills watershed (CZO, Pennsylvania) (Brantley et al.,
2018), East River and Coal Creek watersheds (Colorado) (Zhi et al.,
2019), Sleepers River Research Watershed (Vermont) (Armfield et al.,
2018; Sebestyen et al., 2008), Hubbard Brook Experimental Forest W1 and
W6 watershed (LTER, New Hampshire) (Dittman et al., 2007; LaBaugh et
al., 2013). The four urban sites are Dead Run, Glyndon, Gwynnbrook, and
Carroll Park (Duncan et al., 2017; Kaushal et al., 2008), all located
within the Gwynns Fall Watershed (Baltimore Ecosystem Study LTER,
Maryland).

\textbf{Land use classification} . The NAWQA sites were classified based
on the percentage of major land use according to the USGS criteria
(Spahr et al., 2010). Agricultural sites have \textgreater{} 50\%
agricultural land and \textless{}= 5\% urban land; urban sites have
\textgreater{} 25\% urban and \textless{}= 25\% agricultural land;
undeveloped sites have \textless{}= 5\% urban and \textless{}= 25\%
agricultural land; all other combinations of urban, agricultural, and
undeveloped land are classified as mixed. Land use in the WQW sites was
compiled from multiple sources including EPA water quality reports, USGS
watershed management plan, USDA watershed assessment report, and the
National Land Cover Database (Hyer et al., 2016; Machesky et al., 2010;
Tedesco et al., 2005). Out of the 228 sites, the numbers of Agriculture,
Mixed, Undeveloped, and Urban sites are 61, 78, 49, and 40,
respectively.

\textbf{Concentration-discharge (C-Q) data analysis and estimation of
shallow and deep water concentrations.} A variety of metrics have been
proposed for characterizing export regimes, including load-discharge
relationship, temporal inequality of load and discharge, and various
forms of concentration and discharge relationships (Basu et al., 2010;
Jawitz \& Mitchell, 2011; Thompson et al., 2011; Zhang, 2018). Here we
used the power law form of \(C=aQ^{b}\), a common form for
cross-site comparison (Godsey et al., 2009; Herndon et al., 2015). Here
\(C\) is the nitrate concentration, \(Q\) is
discharge, and \(a\) and \(b\) are fitting
parameters with \emph{b} quantifying the C-Q slopes. For the NAWQA sites
with monthly data, C-Q slopes were calculated using all available paired
daily data over multiple years. Duplicate and incomplete records and
concentrations below the Method Detection Limit~(MDL) of 0.01 mg/L were
excluded for C-Q fitting (Nolan and Stoner, 2000). Flow-weighted stream
concentrations (C\textsubscript{stream} in Figures 2 and 3) were
calculated by dividing total nitrate load by total discharge for the
entire data period.

The WQW sites have data up to 15 min frequency for multiple years (2-4
years typically). High-frequency data were aggregated to daily means of
discharge and nitrate concentrations. They were then used to calculate
the slope \(b\) for each year so that the \emph{b} values
were at consistent time scale as those from NAWQA. The mean and standard
deviation of slope \emph{b} are reported in Figure 5 and 6. In addition,
these high-frequency stream data were used to infer end-member
concentrations in shallow and deep waters. Shallow soil water often
dominates at high flow conditions, whereas deeper base flow predominates
in low flow conditions. Based on this, deep water concentrations
(C\textsubscript{dw}) were estimated as the average of stream
concentrations during low flows at the 5\textsuperscript{th} percentile
of annual discharge. Shallow water concentrations (C\textsubscript{sw})
were estimated as the average of stream concentrations during high flows
(95\textsuperscript{th} percentile of annual discharge). The
concentration ratio (C\textsubscript{ratio}) was calculated as
C\textsubscript{ratio} = C\textsubscript{sw} / C\textsubscript{dw}.

\textbf{Validation of shallow and deep water chemistry estimation.}For
shallow water concentrations, we used a recent study that compiled
measured stream water and soil water concentration from 94 sites with
diverse land cover (e.g., agriculture, forest, disturbed forest,
grassland, urban) across the globe (Sudduth et al., 2013). To be
consistent in land use classification, we grouped ``forest'',
``disturbed forest'', and ``grassland'' in Sudduth et al. (2013) as
``undeveloped'' in this work. With limited measurements of
C\textsubscript{sw}, we cannot directly compare estimated and measured
C\textsubscript{sw} in individual sites. Instead, we examined whether
the relationship between estimated C\textsubscript{sw} and
C\textsubscript{stream} from this study is close to the derived
relationship between measured C\textsubscript{sw} and
C\textsubscript{stream} from Sudduth et al. (2013). The relationships
would be similar if the end-member estimation could provide a realistic
estimation of C\textsubscript{sw} (Figure 4a).

For deep water chemistry, we used the National Ground-Water Monitoring
Network (NGWMN), a compiled database for groundwater level and water
quality (https://cida.usgs.gov/ngwmn/). Based on well location and data
availability, groundwater nitrate data were retrieved from the Midwest
(i.e., Minnesota, Wisconsin, Iowa, Nebraska, Illinois), New Jersey, and
Florida for the maximum coverage. Based on the proximity of sites, a
number of sub-regions were selected for comparing estimated vs. measured
C\textsubscript{dw} (Figure S5). Within each sub-region, estimated and
measured C\textsubscript{dw} were averaged and compared (Figure 4b).

\textbf{Watershed reactive transport modeling.} Watershed modeling was
carried out using a recently developed code BioRT-Flux-PIHM, to explore
the effects of subsurface vertical profile for N (Zhi et al., 2019). A
base model was set up at the Conewago Creek watershed, a tributary to
the Susquehanna River in the Chesapeake Bay. It drains an area of 136
km\textsuperscript{2} and is dominated by agricultural (47\%), forested
(43\%), and developed (17\%) lands. The base model reproduced the
temporal discharge and the slightly flushing C-Q pattern as observed at
the site (Figure S6 and S7). Monte-Carlo simulations of 500 cases were
carried out to cast the base model results to broader conditions
spanning a range of N vertical distribution and reaction kinetics
representing different land use conditions (Figure 6). We then derived a
general relationship for the dependence of C-Q \emph{b} slope on the
ratio of shallow versus deep water concentration (C\textsubscript{sw} /
C\textsubscript{dw}). This simple relationship was further verified
by\emph{b} values estimated based on stream water chemistry from 81
sites, and a few intensively monitored sites with detailed soil and
groundwater measurements. Other details of model setup and calibration,
and Monte-Carlo simulation are in the Supporting Information (SI).

\textbf{RESULTS}

\textbf{Nitrate concentrations under different land uses.} Nitrate data
from 228 sites across the contiguous United States were analyzed (See
Detailed Methods). This includes 147 sites from the National
Water-Quality Assessment (NAWQA) project (Spahr et al., 2010) and 71
sites from the Water Quality Watch (WQW)
(https://waterwatch.usgs.gov/wqwatch/) within the United States
Geological Survey (USGS), and 10 sites from intensively monitored sites
via research networks in the U.S.. Higher stream nitrate concentrations
generally cluster in the Midwestern ``Corn Belt'' that is dominated by
agricultural and mixed lands, whereas undeveloped and urban sites
consistently have lower nitrate concentrations (Figure 2). The soil
water and groundwater nitrate concentrations (Figure S1, 81 sites)
estimated from stream concentrations generally follow the spatial
pattern of stream nitrate concentration, indicating the penetration of
nitrate into deeper groundwater.

Statistical boxplot (Figure 3) confirmed that stream nitrate
concentrations (median and mean in mg/L) decrease from Agriculture (3.2,
4.0), to Mixed (1.4, 2.5), to Urban (0.58, 0.94), and to Undeveloped
(0.16, 0.26). The deep water concentrations (mg/L) have a sequence
similar to stream water concentrations, decreasing from Agriculture
(2.9, 3.0) to Urban (1.2, 1.7) to Undeveloped land (0.13, 0.25). These
values are close to the national median groundwater concentration (mg/L)
of agricultural land (3.1), urban land (1.4), and major aquifer wells
(0.56) (Burow et al., 2010). Compared to pristine sites, agriculture and
urbanization have elevated nitrate level by 3 to 20 times. Agriculture
and Mixed lands also show much higher shallow water concentrations
(C\textsubscript{sw}) compared to deep water concentrations
(C\textsubscript{dw}), whereas Undeveloped lands exhibit the opposite
trend. The Urban sites are the only land type with
C\textsubscript{dw}\textgreater{} C\textsubscript{sw}. Stream and soil
water concentrations correlate strongly with the percentage of
agricultural land with the correlation C\textsubscript{stream} =
0.28e\selectlanguage{ngerman}\textsuperscript{0.031 × \%Ag} (R\textsuperscript{2} = 0.54) and
C\textsubscript{sw} = 0.35e\textsuperscript{0.034 × \%Ag}
(R\textsuperscript{2} = 0.54), respectively. Deep nitrate concentrations
do not correlate well with agricultural area fraction
(R\textsuperscript{2} = 0.25), indicating that other factors might be
important in determining nitrate in deep water (Figure S3).

\textbf{Estimated versus measured concentrations in shallow and deep
waters.} Co-located measurements of stream, soil water, and groundwater
concentrations are rare, except in a few intensively measured sites.
End-member estimation such as hydrograph separation and chemical mixing
analysis are often used to quantify baseflow and quick flow contribution
(Raffensperger et al., 2017) and solute concentrations in end-members
(Miller et al., 2017). To validate the estimation of C\textsubscript{sw}
and C\textsubscript{dw}using stream concentration at the
95\textsuperscript{th} and 5\textsuperscript{th} percentiles, we compare
estimation from this work (Figure 3a and Figure S1) against literature
data and the National Ground-Water Monitoring Network (NGWMN) database
(see Detailed Methods). Results show that the linear relationship
between C\textsubscript{stream} and estimated C\textsubscript{sw} from
this work (slope = 0.9, R\textsuperscript{2} = 0.94) is close to the
linear relationship estimated for 94 sites of diverse land uses covering
4 orders of magnitude (Sudduth et al., 2013) (Figure 4a), suggesting a
strong correlation between C\textsubscript{sw}and
C\textsubscript{stream}. The agricultural sites dominate the top right
corner with high C\textsubscript{sw} and C\textsubscript{stream}.
Comparison of estimated and measured groundwater concentrations as a
proxy for C\textsubscript{dw} exhibits a consistent spatial pattern from
closely located sites in the NGWMN (Figure S5). Most estimated and
measured C\textsubscript{dw} falls on the 1:1 line in Figure 4b,
indicating a close match (Figure 4b). The relatively large error bars of
measured C\textsubscript{dw} may be due to varying sampling depths in
different wells.

\textbf{Nitrate export patterns under different land uses.} Figure 5
(and Figure S2) showed that although different export patterns (i.e.,
multiannual scale C-Q relationship) occur in all land use conditions,
more agricultural sites have high \emph{b} values and more urban sites
have low \emph{b} values. In other words, flushing patterns dominate in
agriculture and mixed lands. This challenges the existing perception
that agricultural lands typically have chemostatic or biogeochemical
stationary patterns (Basu et al., 2010; Basu et al., 2011; Thompson et
al., 2011). Diamond and Cohen (2018) also found that a greater
agricultural land cover was not associated with chemostasis patterns.
For urban watersheds, chemostasis and dilution patterns are most
commonly observed but flushing also occurs. In pristine sites, both
chemostasis and flushing patterns are common.

\textbf{Concentration contrasts in shallow and deep waters drive C-Q
patterns.} To explain nitrate export patterns under different land use
conditions, we resort to a recently-developed, process-based watershed
reactive transport model BioRT-Flux-PIHM (Zhi et al., 2019) to explore
the drivers of C-Q patterns. The model was set up first to reproduce
hydrological and nitrate data in Conewago Creek, a watershed with
\textasciitilde{} 47\% agricultural land in the Chesapeake Bay. The
model has three major water components contributing to the stream:
surface runoff, soil water, and deep water. The shallow water combines
surface runoff and shallow soil interflow water, and the deep water is
the groundwater that interacts with the stream. Annually the shallow and
deep water comprise 80.7\% and 19.3\% of the total discharge,
respectively, although at daily scale the deep water can vary from 1\%
under dry conditions to 99\% during large storm events. The hydrology
aspect of the model captured the general trend that shallow water
dominates stream discharge at high flow and deep water predominates at
low flow. The model reproduced the daily discharge dynamics that was
highly responsive to precipitation, as well as daily concentration
range, temporal trends, and C-Q patterns (Figure S6). The slope
\emph{b}from the C-Q data is 0.13, and the annual concentrations of
shallow water (C\textsubscript{sw}) and deep water (C\textsubscript{dw})
were 3.8 and 2.1 mg/L, respectively, yielding a C\textsubscript{ratio}
of 1.8.

To cast the model to broader conditions, 500 Monte-Carlo simulations
were run using the hydrology and reaction kinetic conditions from the
base case and variable soil N concentrations that yield the
representative range of shallow and deep waters in the four land use
types in Figure 1. The model reproduced the range of soil water and
groundwater concentrations and C\textsubscript{ratio} in Figure 2 (see
Methods and Figure S7 for C-Q relationships in four cases representing
four different land uses). In each simulation case, the model outputs of
concentrations and discharge were used to calculate \emph{b} values and
the annual average C\textsubscript{sw}, C\textsubscript{dw}, and
C\textsubscript{ratio} (Figure S8). Model results from the 500 cases
(overlapping gray circles, Figure 6) collapse to an ``S'' shaped curve
(\(b=\frac{\delta_{b}\ C_{\text{ratio}}}{C_{\text{ratio},\ 1/2}\ +\text{\ C}_{\text{ratio}}\ }+b_{\min}\)) that illustrates the consistent dependence of
\emph{b} values on C\textsubscript{ratio}, where \(\delta_{b}=1.66\) is
the difference between maximum \emph{b} (\(b_{\max}=0.73\)) and minimum
\emph{b}(\(b_{\min}=-0.93\)) and \(C_{\text{ratio},\ 1/2}=0.80\) is the concentration
ratio when \emph{b} = ½ (\(b_{\max}+b_{\min}\)). This modeled ``S'' curve
explained the \emph{b} values from 81 sites with available
C\textsubscript{ratio} data. The curve shows that most agricultural and
mixed lands have high \emph{b} values whereas urban and undeveloped
lands have lower \emph{b} values.

\textbf{DISCUSSION}

\textbf{The significance of N distribution and processes in the shallow
and deep subsurface.} Figure 6 suggests that the subsurface N
concentration contrast (or C\textsubscript{ratio}) arising from
different land uses predominantly control export patterns. The
representatives of shallow and deep waters, such as soil and
groundwater, have been shown before to correlate well with stream
nitrate concentration across different land uses (Sudduth et al., 2013).
The export pattern dependence on concentration contrast echoes the
decades-long End Member Mixing Approach (EMMAS) that uses end-member
concentrations to infer stream water chemistry (Hooper et al., 1990).
The spatial data synthesis across the U.S. continent from this work
enables the generalization of this idea across diverse climate, geology,
and land use conditions. The physics-based watershed reactive transport
modeling facilitates the derivation of the general equation that can be
used to estimate export patterns based on streamflow concentrations or
measured shallow and deep water concentrations.

These shallow versus deep physical contrasts originate from chemical and
physical weathering in pristine sites (Brantley et al., 2017). In urban
and agriculture sites, these contrasts arise from human engineering
efforts with tile drains, impervious surfaces, and water pipes that
modify environments (Grimm et al., 2008). The concentration differences
in shallow and deep waters arise not only from different subsurface
distribution of N source but also from different processes. In
agricultural lands, nitrate concentrations in shallow waters are high
not only because abundant N sources leach nitrate; it is also because
denitrification is limited with the presence of abundant
O\textsubscript{2}. In deeper groundwater, with limited
O\textsubscript{2}, denitrification can occur and reduce nitrate (Kolbe
et al., 2019).

Under broad conditions, the relative magnitude of shallow versus deep
water concentrations may hinge on soil properties and geologic
structure. These subsurface structures determine shallow and deep
connectivity, recharge and water table depth (Brantley et al., 2017), as
well as local biogeochemical conditions (e.g., anoxic condition, organic
carbon availability) that control the extent of denitrification and
nitrate removal (Kolbe et al., 2019; Tesoriero et al., 2015). Miller et
al. (2017) showed that nitrate export exhibits a dilution pattern in
Tomorrow River, a site with permeable sand and gravel, but a flushing
pattern in Duck Creek with low-permeability clayey soil clayey soils
developed from glacial tills. The overwhelming convergence toward
high\emph{b} values in agriculture lands, however, indicates that the
concentration contrasts in shallow versus deep waters are the
predominant control of export patterns.

In urban watersheds, we observe varied export patterns and a large
number of sites have higher concentrations in deep water compared to
shallow water. This is possibly caused by underground leaky sewage and
pipes that can contaminate groundwater and surface water (Divers et al.,
2013; Lerner et al., 1999; Pennino et al., 2016). A large number of
urban watersheds also exhibit chemostatic patterns, potentially due to
the co-occurrence of both shallow N sources (e.g., lawn fertilizer, pet
waste, atmospheric deposition, automobile emission) and deeper
underground input from buried sewage and septic systems. Kaushal et al.
(2014) and Newcomer Johnson et al. (2014) suggest that urban N export
can be influenced by the degree of hydrologic connectivity associated
with impervious surface, stormwater infrastructure and sewage pipe, and
stream restoration. Point source discharge from wastewater treatment
plant (WWTP) can also increase nutrient loading (Luthy et al., 2015),
potentially contributing to varied export patterns. In arid and semiarid
regions with small urban streams, high nitrate concentration from WWTP
can dominate the base flow at the dry time and become diluted under high
flow conditions, resulting in dilution pattern (Marti et al., 2010). In
fact, urban watersheds are complicated as the groundwater-soil-surface
water interactions are modified by the level of urbanization and
management, and non-point sources from leaky infrastructure and chronic
groundwater contamination (Kaushal \& Belt, 2012).

\textbf{Limitations, simplifications, and uncertainties} . The
conceptual model in Figure 1 and the general b equation in Figure 6
emphasize two end members in shallow and deep zones and are meant to
build a simple relationship between export patterns and major components
of water contributing to the stream. Here we lump the waters into
shallow and deep waters components to illustrate the first-order
control. Such simplification is often necessary in practice. In urban
watersheds, for example, the impervious surface often contributes to
large surface runoff in storms. Surface runoff however typically
dominates as short pulses in early stages of storm events often
resulting in an overall limited contribution of surface runoff to annual
discharge (e.g., 11\% by Pellerin et al. (2008). These temporarily large
contributions are followed by rapid subsurface flow (via underground
stormwater pipes) with much longer duration. Field studies typically
have incomplete information (e.g., hydrograph separation, isotopic
signature) of contributing flow paths. It is therefore often necessary
to lump different water sources into major compartments (Barnes \&
Raymond, 2010), in order to capture average behaviors.

With its simplicity, the model does not take into account specifics of
individual sites, which can lead to deviations from the b curve in
Figure 6. For example, in some places, N is distributed and processed in
a way that demands more than two end members (Cowie et al., 2017; Miller
et al., 2017). The model also does not explicitly account for N removal
in streams. Instream removal depends on a wide variety of parameters
including local climate and seasonality, landscape structure (e.g.,
topography, hyporheic zone), and biogeochemical conditions (e.g.,
nitrate legacy, stream carbon source) (Dodds et al., 2002; Hill, 1996;
Mulholland et al., 2008; Vidon \& Hill, 2004). Significant instream N
removal can lead to underestimation of C\textsubscript{sw} and
C\textsubscript{dw} but to a different extent. As C\textsubscript{sw} is
estimated under high flow conditions where instream N removal is often
not efficient, it may not be heavily influenced by instream processes
(Dodds et al., 2002). The estimation of deep water concentrations
however is based on low flow conditions where instream N removal can be
highly effective such that C\textsubscript{dw} may be underestimated.
This possible different extent of underestimation can shift
C\textsubscript{ratio} (= C\textsubscript{sw} / C\textsubscript{dw}) and
b values to higher values, moving data points toward the top right end
of the d b curve. This may be the reason that some data points from Ag
and Mixed fall in the right-hand side of the b curve in Figure 6.

\textbf{Inferring shallow and deep water chemistry from stream
chemistry.} This study indicates that we can infer shallow and deep
water chemistry from streamflow chemistry under different flow regimes.
This is important as detailed subsurface characterization and
concentration measurements are often limited to only a few long-term
monitoring sites in developed countries (Brantley et al., 2018; Gran et
al., 2019), although we often claim an era of ``big data''. As shown in
Figure 6, few intensively measured watersheds have both soil water and
groundwater chemistry measurements. The extrapolation of shallow and
deep water chemistry from high flow and low flow stream chemistry
therefore enables estimation of water chemistry (without digging holes
in the ground) to infer possible export patterns and loads. This
approach can be applied not only for nitrate, but also for other solutes
in general.

\textbf{N export under changing climate.} The agriculture sites have
increased nitrate concentrations in shallow and deep waters by 15-20
times and 9-16 times, respectively, compared to undeveloped sites
(Figure 3a). In fact, extensive tile drainage networks in agricultural
lands (e.g., 80\% of the landscape in the Midwest) shortcut the shallow
water directly to stream, lowers water table (Blann et al., 2009; King
et al., 2015). In effect, these draining tiles decrease vertical
connectivity to deeper aquifers (Association, 2018) and therefore
increase concentration contrasts between shallow water (e.g., soil water
and tile drainage water) and deep groundwater, inducing pronounced
flushing patterns (Figure 6). The flushing pattern indicates that export
is sensitive to large hydrological events such as flooding, which has
been predicted to intensify as the pace of climate change accelerates
(Prein et al., 2017).

On the other hand, if high nitrate in shallow water is redirected more
to the deeper subsurface via higher vertical connectivity, the water
will reroute through longer flow paths, enhancing transformation via
denitrification and nitrate removal (Kolbe et al., 2019). The extent of
such transformation will depend on local conditions. The longer and
slower groundwater flow paths however will also lead to decade-to
century-scale time lags between anthropogenic N inputs and riverine
outputs (Sebilo et al., 2013; Van Meter \& Basu, 2015). This can present
significant challenges for balancing nitrogen budget (e.g., the
`missing' N in mass-balance) (Boyer et al., 2002) and management
effectiveness (e.g., land-use management, N-loading reduction) (Sebilo
et al., 2013; Van Meter et al., 2017). In addition, this downward N flux
also raises tantalizing questions about long-term subsurface structure
and functioning. How and how much do these man-made, high nutrient
levels elevate carbon effluxes (Zamanian et al., 2018), accelerate
weathering processes, and alter watershed functioning (Kaushal et al.,
2013; Perrin et al., 2008)? These broader earth system responses can
have far-reaching impacts on carbon and nutrient cycles at the global
scale.

\textbf{Acknowledgments}

This study was supported by the U.S. Department of Energy (DOE)
Subsurface Biogeochemical Research (SBR) program (DE\selectlanguage{english}-SC0016221) and the
Penn State Institute of CyberScience.

\textbf{Figures in the main text}\selectlanguage{english}
\begin{figure}[H]
\begin{center}
\includegraphics[width=0.70\columnwidth]{figures/image1/image1}

\end{center}
\end{figure}

\textbf{Figure 1} . A conceptual figure and hypothesis about nitrate
export patterns, exemplified in Concentration-Discharge (C-Q)
relationship under different land use conditions (i.e., Ag =
Agriculture, Undev = Undeveloped, and Urban). Triangle wedges at the
left of each hillslope describe the depth profile of N abundance with
dark and light red represents high and low N abundance, respectively.
Here we define two ``loose'' end-members: the ``shallow'' water (from,
for example, surface runoff, shallow soils, or shallow pipes) and the
``deep'' water (e.g., from deeper zones such as groundwater). Stream
water is often dominated by shallow water under high flow conditions and
by deep water under low flow conditions. \emph{We propose the shallow
versus deep hypothesis: nitrate concentration contrasts in shallow water
(C\textsubscript{sw}) and deep water (C\textsubscript{dw}) drive the
export patterns.} If this hypothesis is true, we expect higher
concentrations in shallow water in agricultural lands, leaning toward a
flushing pattern. In contrast, in urban watersheds, concentrated
nutrients accumulated in leaky pipes in deeper subsurface are often
higher than shallow, rapid runoff on impervious surfaces, possibly
leaning toward a dilution pattern. Nitrogen in forests and pristine
sites can come from the decomposition of organic matters in shallow
soils and leaching from N-containing rocks in deeper zones. They are
often tightly cycled with very low concentrations. These characteristics
can lead to diverse patterns\selectlanguage{english}
\begin{figure}[H]
\begin{center}
\includegraphics[width=0.70\columnwidth]{figures/image2/image2}

\end{center}
\end{figure}

\textbf{Figure 2} . The spatial pattern of stream nitrate
(NO\textsubscript{3}-N) concentration under Ag (agriculture, 61 sites),
Mixed (Ag + other land use, 78 sites), Undev (Undeveloped, pristine, 49
sites), and Urban (40 sites) conditions. Please refer to Detailed
Methods for the classification of land uses. Compared to undeveloped
sites, agricultural lands generally have elevated stream water
concentrations by more than one order of magnitude.\selectlanguage{english}
\begin{figure}[H]
\begin{center}
\includegraphics[width=0.70\columnwidth]{figures/image3/image3}

\end{center}
\end{figure}

\textbf{Figure 3.} Nitrate (NO\textsubscript{3}-N) concentrations: (a)
Statistics of stream (C\textsubscript{stream}, left box), shallow water
(C\textsubscript{sw}, middle box), and deep water (often groundwater,
C\textsubscript{dw}, right box) under different land use conditions. The
C\textsubscript{stream} was the direct field measurement from 228 sites
while C\textsubscript{sw} and C\textsubscript{dw} were estimated and/or
measured end-member concentrations from 81 sites. The central line in
each box indicates the median, and the bottom and top edges are the
25\textsuperscript{th} and 75\textsuperscript{th} percentiles,
respectively. The whiskers extend to the most extreme data points not
considered outliers (not shown). Outliers are defined as data points
more than~1.5 interquartile ranges~(IQR, Q3-Q1) below the first~quartile
(Q1)~or above the third~quartile (Q3), i.e., \textless{} Q1-1.5\selectlanguage{ngerman}×IQR or
\textgreater{} Q3 + 1.5×IQR. (b) C\textsubscript{stream}as a function of
agricultural land cover. The solid gray line is the fit equation
C\textsubscript{stream} = 0.28e\textsuperscript{0.031 ×
\%Ag}(R\textsuperscript{2} = 0.54), and the dashed gray lines are
confidence interval at 99.9\%, with upper and lower bound equations of
C\textsubscript{stream} = 0.39e\textsuperscript{0.038 × \%Ag} and
C\textsubscript{stream} = 0.20e\textsuperscript{0.025 × \%Ag},
respectively.\selectlanguage{english}
\begin{figure}[H]
\begin{center}
\includegraphics[width=0.70\columnwidth]{figures/image4/image4}

\end{center}
\end{figure}

\textbf{Figure 4.} Validation of estimated C\textsubscript{sw} and
C\textsubscript{dw} (this study). (a) A linear relationship (red line
with slope = 0.95 and R\textsuperscript{2} = 0.94) derived from
estimations in this work (colored markers), compared to the linear
relationship (black line, slope = 0.83, R\textsuperscript{2} = 0.70)
from measurements (gray markers) compiled in Sudduth et al. (2013). The
dotted gray circle at the top left corner is an outliner and was
excluded from data fitting. The similar slope and correlation indicate
that C\textsubscript{sw}estimated from stream chemistry approximates
measured C\textsubscript{sw}. The C\textsubscript{stream} from this
study are flow-weighted concentrations calculated by dividing annual
nitrate load by annual discharge. (b) Estimated C\textsubscript{dw} from
this work vs. measured groundwater concentration (as a proxy for
C\textsubscript{dw}) from the National Ground Water Monitoring Network.
The measured groundwater concentrations were averaged for close by
regions (Figure S5). Most symbols fall close to the 1:1 line, indicating
estimated C\textsubscript{dw} can be used to approximate groundwater
concentration. Error bars are for one standard deviation in nearby
regions.\selectlanguage{english}
\begin{figure}[H]
\begin{center}
\includegraphics[width=0.70\columnwidth]{figures/image5/image5}

\end{center}
\end{figure}

\textbf{Figure 5.} Boxplot of slope \emph{b} under different land use
conditions. The slope \emph{b} was calculated from multiple years of
low-frequency concentration and discharge data and high-frequency data
averaged to daily values so the focus here is not event-based C-Q (see
Detailed Methods). The gray dashed lines are thresholds separating
flushing (b\textgreater{}0.1) and dilution patterns (b\textless{}-0.1)
from the chemostatic pattern, using b = 0.1 and =0.1 as thresholds
(Herndon et al., 2015; Hoagland et al., 2017). The gray numbers are the
fractions of flushing, chemostatic, and dilution patterns. Two-sample
\emph{t} -tests were performed for statistical analysis: the single star
marker indicates statistical significance at the level of 0.05, and the
double star maker indicates statistical significance at the level of
0.01.\selectlanguage{english}
\begin{figure}[H]
\begin{center}
\includegraphics[width=0.70\columnwidth]{figures/image6/image6}

\end{center}
\end{figure}

\textbf{Figure 6} . The dependence of export pattern, or C-Q
slope\emph{b} , on the ratio of shallow water versus deep water
concentrations (C\textsubscript{ratio} = C\textsubscript{sw} /
C\textsubscript{dw}). The gray line is an ensemble of 500 gray circles
from 500 Monte-Carlo simulations under different land use conditions.
Each gray circle represents one simulation with a calculated slope
\emph{b} from simulated C and Q values and C\textsubscript{ratio} from
an annual average of simulated C\textsubscript{sw} and
C\textsubscript{dw} (see Figure S8). The triangles and stars are from
individual watersheds: some are based on measured stream discharge,
chemistry, and soil and groundwater concentrations
C\textsubscript{ratio} (triangles), and others are based on
high-frequency stream data (stars). The error bars are for one standard
deviation. The black line represents the fitting
equation\(b=\frac{\delta_{b}\ C_{\text{ratio}}}{C_{\text{ratio},\ 1/2}\ +\text{\ C}_{\text{ratio}}\ }+b_{\min}\), which describes the black line and the
dependence of \emph{b} on C\textsubscript{ratio}. Here
\(\delta_{b}\)is the difference between maximum\emph{b}
(\(b_{\max}\)) and minimum \emph{b} (\(b_{\min}\)),
and\(C_{\text{ratio},\ 1/2}\) is the concentration ratio when \emph{b} = \selectlanguage{ngerman}½
(\(b_{\max}+b_{\min}\)). Here \(b_{\max}=0.73\) and\(b_{\min}=-0.93\),
\(\delta_{b}=1.66\), and\(C_{\text{ratio},\ 1/2}=0.80\). The agricultural and mixed
lands cluster toward high b and C\textsubscript{ratio} values,
indicating flushing as the predominant export pattern. The urban
watersheds lean toward chemostasis (near zero \emph{b} ) and dilution
(negative \emph{b} ) patterns. The undeveloped watersheds cover a wide
range with more flushing than dilution patterns.

\textbf{Supporting Information}

\begin{itemize}
\tightlist
\item
  \textbf{Supplementary Methods}
\item
  \textbf{Supplementary Figures: 8}
\item
  \textbf{Supplementary Tables: 1}
\end{itemize}

\textbf{Supplementary Methods}

\textbf{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\selectlanguage{english}-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
km\textsuperscript{2} 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 degC, 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.

\textbf{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 \emph{r} is the reaction rate
(mol/m\textsuperscript{3}/s); \(k\)is the rate constant
(10\textsuperscript{-9}, mol/m\textsuperscript{2}/s); and A is the
effective surface area for nitrate leaching (m\textsuperscript{2}),
which is calculated based specific surface area (SSA,
m\textsuperscript{2}/g), volume fraction (0.02
m\textsuperscript{3}/m\textsuperscript{3}), and solid density (urea,
1.34 g/cm\textsuperscript{3}). 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 Q\textsubscript{10}-based
form\(f\left(T\right)=\ Q_{10}^{\left|T-20\right|/10},\ \) where\(T\) is the soil temperature
(degC) and \(Q_{10}=2\) is the relative increase in reaction rates
when temperature increases by 10 degC (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 \selectlanguage{greek}\emph{α} \selectlanguage{english}and \emph{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 m\textsuperscript{2}/g that is similar to the reported range of 1.7
to 3.5 m\textsuperscript{2}/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.

\textbf{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 (C\textsubscript{dw}, 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 (C\textsubscript{sw}) also
increase because more NO\textsubscript{3}\textsuperscript{-} 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 C\textsubscript{ratio} used to plot the \emph{b}
versus C\textsubscript{ratio} in Figure 6.

In addition, four hypothetical cases for Agriculture, Mixed,
Undeveloped, and Urban were simulated based on available
C\textsubscript{sw} and C\textsubscript{gw} 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 C\textsubscript{ratio} (= C\textsubscript{sw} /
C\textsubscript{dw}).

\textbf{Supplementary Figures (Figure S1 -- S8)}\selectlanguage{english}
\begin{figure}[H]
\begin{center}
\includegraphics[width=0.70\columnwidth]{figures/image7/image7}

\end{center}
\end{figure}

\textbf{Figure S1.} The spatial pattern of shallow water
(C\textsubscript{sw}) and deep water (C\textsubscript{dw}) nitrate
concentration under different land use conditions (81 sites). Note that
filled circles are estimated end-member concentrations from
high-frequency stream nitrate (C\textsubscript{stream}, Figure 2) while
star markers are direct field measurements of C\textsubscript{sw} and
C\textsubscript{dw} from intensively studied sites.\selectlanguage{english}
\begin{figure}[H]
\begin{center}
\includegraphics[width=0.70\columnwidth]{figures/image8/image8}

\end{center}
\end{figure}

\textbf{Figure S2.} The spatial pattern of C-Q slope \emph{b} (228
sites), which was calculated by fitting the power law
\(C=aQ^{b}\)using all available censored data of nitrate
concentration (C) and discharge (Q).\selectlanguage{english}
\begin{figure}[H]
\begin{center}
\includegraphics[width=0.70\columnwidth]{figures/image9/image9}

\end{center}
\end{figure}

\textbf{Figure S3.} Shallow water (C\textsubscript{sw}) and deep water
(C\textsubscript{dw}) nitrate concentration as a function of percent of
agricultural land (81 sites). The solid gray line is the exponential fit
while dashed lines are upper and lower bounds at the confidence interval
of 95\%. Note that filled circles are estimated end-member concentration
from high-frequency stream nitrate (C\textsubscript{stream}, see
Methods) while stars are direct field measurements of
C\textsubscript{sw} and C\textsubscript{dw} from intensively studied
sites.\selectlanguage{english}
\begin{figure}[H]
\begin{center}
\includegraphics[width=0.70\columnwidth]{figures/image10/image10}

\end{center}
\end{figure}

\textbf{Figure S4.} Shallow water (C\textsubscript{sw}) vs. deep water
(C\textsubscript{dw}) nitrate concentration (Figure S1\textbf{,} 81
sites). Filled circles are estimated end-member concentrations from
high-frequency stream nitrate (C\textsubscript{stream}) while stars are
direct field measurements of C\textsubscript{sw} and
C\textsubscript{dw}from intensively studied sites.\selectlanguage{english}
\begin{figure}[H]
\begin{center}
\includegraphics[width=0.70\columnwidth]{figures/image11/image11}

\end{center}
\end{figure}

\textbf{Figure S5.}Validation of estimated deep groundwater nitrate
concentration (C\textsubscript{dw}) against field measurements from the
National Ground-Water Monitoring Network. Within each red circle,
estimated and measured C\textsubscript{dw} are first averaged and then
compared (results see Figure 4b). Red circles are determined by the
proximity of location and data availability (i.e., at least one
estimated C\textsubscript{gw} point for comparison with measurements).
Iowa and Florida have more direct comparison because there are more
estimated C\textsubscript{gw} points (triangle markers) located in these
regions and more measured C\textsubscript{gw} points (circle markers) in
proximity.\selectlanguage{english}
\begin{figure}[H]
\begin{center}
\includegraphics[width=0.70\columnwidth]{figures/image12/image12}

\end{center}
\end{figure}

\textbf{Figure S6.} Water fluxes and stream nitrate in 2017: discharge
(data vs. model, left-Y axis), precipitation (data,
right-1\textsuperscript{st} Y-axis), stream nitrate (data vs. model,
right-2\textsuperscript{nd} Y-axis) at Conewago Creek watershed. Annual
precipitation and discharge of 2017 was 1.04 m and 0.47 m, respectively,
resulting in a runoff ratio of 0.45.\selectlanguage{english}
\begin{figure}[H]
\begin{center}
\includegraphics[width=0.70\columnwidth]{figures/image13/image13}

\end{center}
\end{figure}

\textbf{Figure S7.} Modeled C-Q relationships from 4 hypothetical cases
reproduced typical nitrate export patterns under Agriculture, Mixed,
Undeveloped, and Urban land uses. The 4 hypothetical cases used
representative soil N distribution to reproduce measured shallow water
(C\textsubscript{sw}) and deep water (C\textsubscript{dw}) nitrate under
each land use condition (Figure 3). All colored dots are from models.
Open pink circles are measured stream nitrate data measured from
Conewago Creek watershed. Dashed lines are C-Q fits for each model case.\selectlanguage{english}
\begin{figure}[H]
\begin{center}
\includegraphics[width=0.70\columnwidth]{figures/image14/image14}

\end{center}
\end{figure}

\textbf{Figure S8.} The annual average of shallow water
(C\textsubscript{sw}) versus deep water (C\textsubscript{dw}) nitrate
concentrations from 500 Monte-Carlo simulation cases. Each open circle
represents concentrations from one of the 500 runs. The two gray dashed
lines are the thresholds separating flushing (red,
\emph{b}\textgreater{} 0.1) and dilution (blue, \emph{b} \textless{}
0.1) patterns from chemostatic (gray, -0.1\textless{}= \emph{b}
\textless{}= 0.1) patterns. This figure shows that as long as
C\textsubscript{sw} / C\textsubscript{dw} \textgreater{} 1.7, nitrate
exhibits flushing patterns (red); when C\textsubscript{sw} /
C\textsubscript{dw}\textless{} 0.8, nitrate export exhibits dilution
patterns (blue).

\textbf{Supplementary Tables}

\textbf{Table S1.} Hydrology calibration: key land-surface and soil
parameters\selectlanguage{english}
\begin{longtable}[]{@{}lll@{}}
\toprule
Parameter & Calibrated value & Parameter note\tabularnewline
\midrule
\endhead
\(C_{\text{zil}}\) & 7.13 & Zilitinkevich coefficient\tabularnewline
\(D_{\text{root}}\) & 0.34 & Rooting depth (m)\tabularnewline
RGL & 0.80 & Reference visible solar radiation
(W/m\textsuperscript{2})\tabularnewline
RS & 0.89 & Minimum stomatal resistance
(s/m)\tabularnewline
HS & 1.42 & Water vapor exchange
coefficient\tabularnewline
\(\mathbf{K}_{\mathbf{\text{satH}}}\) & 1.32 & Horizontal saturated hydraulic conductivity
(m/d)\tabularnewline
\(\mathbf{K}_{\mathbf{\text{satV}}}\) & 26.7 & Vertical saturated hydraulic conductivity
(m/d)\tabularnewline
\(K_{\inf}\) & 0.76 & Infiltration hydraulic conductivity
(m/d)\tabularnewline
\(D_{\text{mac}}\) & 0.31 & Macropore depth (m)\tabularnewline
\(\mathbf{\theta}\) & 0.41 & Porosity
(m\textsuperscript{3}/m\textsuperscript{3})\tabularnewline
\(\mathbf{\alpha}\) & 4.67 & van Genuchten alpha, inversely proportional
to pore diameter (m\textsuperscript{-1})\tabularnewline
\(\mathbf{n}\) & 1.88 & van Genuchten n, inversely proportional to
water retention\tabularnewline
\(Z_{\text{soil}}\) & 1.0 & soil thickness (m)\tabularnewline
\bottomrule
\end{longtable}

Note that soil properties are initially estimated using pedotransfer
function (W\selectlanguage{ngerman}östen et al., 2001) based on soil composition derived from
the SSURGO database. Bolded parameters are the most sensitive parameters
that reproduced stream discharge.

\textbf{Reference:}

Armfield, J., Perdrial, J., Gagnon, A., Ehrenkranz, J., Perdrial, N.,
Cincotta, M., Ross, D., Shanley, J., Underwood, K., Ryan, P. 2018. Does
stream water composition at Sleepers River in Vermont reflect dynamic
changes in soils during recovery from acidification? \emph{Frontiers in
Earth Science} , \textbf{6} , 246.

Association, M.G. 2018. Drain tiles and groundwater resources:
Understanding the relations. Minnesota Ground Water Association.

Baker, L.A., Hope, D., Xu, Y., Edmonds, J., Lauver, L. 2001. Nitrogen
balance for the Central Arizona--Phoenix (CAP)
ecosystem.\emph{Ecosystems} , \textbf{4} (6), 582-602.

Bao, C., Li, L., Shi, Y., Duffy, C. 2017. Understanding watershed
hydrogeochemistry: 1. Development of RT\selectlanguage{english}-Flux\selectlanguage{english}-PIHM. \emph{Water Resour.
Res.} , \textbf{53} (3), 2328-2345.

Barnes, R.T., Raymond, P.A. 2010. Land\selectlanguage{english}-use controls on sources and
processing of nitrate in small watersheds: insights from dual isotopic
analysis. \emph{Ecological Applications} , \textbf{20} (7), 1961-1978.

Basu, N.B., Destouni, G., Jawitz, J.W., Thompson, S.E., Loukinova, N.V.,
Darracq, A., Zanardo, S., Yaeger, M., Sivapalan, M., Rinaldo, A. 2010.
Nutrient loads exported from managed catchments reveal emergent
biogeochemical stationarity. \emph{Geophysical Research Letters}
,\textbf{37} (23).

Basu, N.B., Thompson, S.E., Rao, P.S.C. 2011. Hydrologic and
biogeochemical functioning of intensively managed catchments: A
synthesis of top\selectlanguage{english}-down analyses. \emph{Water Resources Research}
,\textbf{47} (10).

Bishop, K., Seibert, J., Kohler, S., Laudon, H. 2004. Resolving the
Double Paradox of rapidly mobilized old water with highly variable
responses in runoff chemistry. \emph{Hydrological Processes}
,\textbf{18} (1), 185-189.

Blann, K.L., Anderson, J.L., Sands, G.R., Vondracek, B. 2009. Effects of
Agricultural Drainage on Aquatic Ecosystems: A Review. \emph{Critical
Reviews in Environmental Science and Technology} , \textbf{39} (11),
909-1001.

Boyer, E.W., Goodale, C.L., Jaworski, N.A., Howarth, R.W. 2002.
Anthropogenic nitrogen sources and relationships to riverine nitrogen
export in the northeastern USA. \emph{Biogeochemistry} , \textbf{57}
(1), 137-169.

Brantley, S.L., Lebedeva, M.I., Balashov, V.N., Singha, K., Sullivan,
P.L., Stinchcomb, G. 2017. Toward a conceptual model relating chemical
reaction fronts to water flow paths in hills. \emph{Geomorphology}
,\textbf{277} , 100-117.

Brantley, S.L., White, T., West, N., Williams, J.Z., Forsythe, B.,
Shapich, D., Kaye, J., Lin, H., Shi, Y., Kaye, M. 2018. Susquehanna
Shale Hills Critical Zone Observatory: Shale Hills in the context of
Shaver's Creek watershed. \emph{Vadose Zone Journal} , \textbf{17} (1).

Burow, K.R., Nolan, B.T., Rupert, M.G., Dubrovsky, N.M. 2010. Nitrate in
groundwater of the United States, 1991\selectlanguage{english}- 2003. \emph{Environmental
Science \& Technology} , \textbf{44} (13), 4988-4997.

Cowie, R.M., Knowles, J.F., Dailey, K.R., Williams, M.W., Mills, T.J.,
Molotch, N.P. 2017. Sources of streamflow along a headwater catchment
elevational gradient. \emph{Journal of Hydrology} , \textbf{549} ,
163-178.

Creed, I., Band, L. 1998. Export of nitrogen from catchments within a
temperate forest: evidence for a unifying mechanism regulated by
variable source area dynamics. \emph{Water Resources Research}
,\textbf{34} (11), 3105-3120.

Creed, I., Band, L., Foster, N., Morrison, I., Nicolson, J., Semkin, R.,
Jeffries, D. 1996. Regulation of nitrate\selectlanguage{english}-N release from temperate
forests: A test of the N flushing hypothesis. \emph{Water Resources
Research} , \textbf{32} (11), 3337-3354.

Diamond, J.S., Cohen, M.J. 2018. Complex patterns of catchment
solute-discharge relationships for coastal plain
rivers.\emph{Hydrological Processes} , \textbf{32} (3), 388-401.

Diaz, R.J., Rosenberg, R. 2008. Spreading dead zones and consequences
for marine ecosystems. \emph{science} , \textbf{321} (5891), 926-929.

Dittman, J.A., Driscoll, C.T., Groffman, P.M., Fahey, T.J. 2007.
Dynamics of nitrogen and dissolved organic carbon at the Hubbard Brook
Experimental Forest. \emph{Ecology} , \textbf{88} (5), 1153-1166.

Divers, M.T., Elliott, E.M., Bain, D.J. 2013. Constraining Nitrogen
Inputs to Urban Streams from Leaking Sewers Using Inverse Modeling:
Implications for Dissolved Inorganic Nitrogen (DIN) Retention in Urban
Environments. \emph{Environmental Science \& Technology} ,\textbf{47}
(4), 1816-1823.

Divers, M.T., Elliott, E.M., Bain, D.J. 2014. Quantification of nitrate
sources to an urban stream using dual nitrate
isotopes.\emph{Environmental science \& technology} , \textbf{48} (18),
10580-10587.

Dodds, W.K., Lopez, A.J., Bowden, W.B., Gregory, S., Grimm, N.B.,
Hamilton, S.K., Hershey, A.E., Marti, E., McDowell, W.H., Meyer, J.L.,
Morrall, D., Mulholland, P.J., Peterson, B.J., Tank, J.L., Valett, H.M.,
Webster, J.R., Wollheim, W. 2002. N uptake as a function of
concentration in streams. \emph{Journal of the North American
Benthological Society} , \textbf{21} (2), 206-220.

Duncan, J.M., Welty, C., Kemper, J.T., Groffman, P.M., Band, L.E. 2017.
Dynamics of nitrate concentration\selectlanguage{english}-discharge patterns in an urban
watershed. \emph{Water Resources Research} , \textbf{53} (8), 7349-7365.

Eghbali Babadi, F., Yunus, R., Abbasi, A., Masoudi Soltani, S. 2019.
Response Surface Method in the Optimization of a Rotary Pan-Equipped
process for Increased Efficiency of Slow-Release Coated
Urea.\emph{Processes} , \textbf{7} (3), 125.

Galloway, J.N., Dentener, F.J., Capone, D.G., Boyer, E.W., Howarth,
R.W., Seitzinger, S.P., Asner, G.P., Cleveland, C.C., Green, P.,
Holland, E.A. 2004. Nitrogen cycles: past, present, and
future.\emph{Biogeochemistry} , \textbf{70} (2), 153-226.

Godsey, S.E., Kirchner, J.W., Clow, D.W. 2009. Concentration-discharge
relationships reflect chemostatic characteristics of US
catchments.\emph{Hydrological Processes} , \textbf{23} (13), 1844-1864.

Gran, K.B., Dolph, C., Baker, A., Bevis, M., Cho, S.J., Czuba, J.A.,
Dalzell, B., Danesh-Yazdi, M., Hansen, A.T., Kelly, S., Lang, Z.,
Schwenk, J., Belmont, P., Finlay, J.C., Kumar, P., Rabotyagov, S.,
Roehrig, G., Wilcock, P., Foufoula-Georgiou, E. 2019. The Power of
Environmental Observatories for Advancing Multidisciplinary Research,
Outreach, and Decision support: The Case of the Minnesota River
Basin.\emph{Water Resources Research} , \textbf{55} (4), 3576-3592.

Grimm, N.B., Faeth, S.H., Golubiewski, N.E., Redman, C.L., Wu, J., Bai,
X., Briggs, J.M. 2008. Global change and the ecology of
cities.\emph{science} , \textbf{319} (5864), 756-760.

Groffman, P.M., Bain, D.J., Band, L.E., Belt, K.T., Brush, G.S., Grove,
J.M., Pouyat, R.V., Yesilonis, I.C., Zipperer, W.C. 2003. Down by the
riverside: urban riparian ecology. \emph{Frontiers in Ecology and the
Environment} , \textbf{1} (6), 315-321.

Hagedorn, F., Bucher, J.B., Schleppi, P. 2001. Contrasting dynamics of
dissolved inorganic and organic nitrogen in soil and surface waters of
forested catchments with Gleysols. \emph{Geoderma} , \textbf{100} (1),
173-192.

Hamamoto, S., Moldrup, P., Kawamoto, K., Komatsu, T. 2010.
Excluded-volume expansion of Archie's law for gas and solute
diffusivities and electrical and thermal conductivities in variably
saturated porous media. \emph{Water Resources Research} , \textbf{46}
(6).

Hararuk, O., Smith, M.J., Luo, Y. 2015. Microbial models with
data-driven parameters predict stronger soil carbon responses to climate
change. \emph{Glob. Chang. Biol.} , \textbf{21} (6), 2439-53.

Herndon, E.M., Dere, A.L., Sullivan, P.L., Norris, D., Reynolds, B.,
Brantley, S.L. 2015. Landscape heterogeneity drives contrasting
concentration--discharge relationships in shale headwater
catchments.\emph{Hydrology and earth system sciences} , \textbf{19} (8),
3333-3347.

Hill, A.R. 1996. Nitrate Removal in Stream Riparian Zones. \emph{Journal
of Environmental Quality} , \textbf{25} , 743-755.

Hoagland, B., Russo, T.A., Gu, X., Hill, L., Kaye, J., Forsythe, B.,
Brantley, S.L. 2017. Hyporheic zone influences on
concentration\selectlanguage{english}-discharge relationships in a headwater sandstone
stream.\emph{Water Resources Research} , \textbf{53} (6), 4643-4667.

Hooper, R.P., Christophersen, N., Peters, N.E. 1990. Modelling
streamwater chemistry as a mixture of soilwater end-members --- An
application to the Panola Mountain catchment, Georgia,
U.S.A.\emph{Journal of Hydrology} , \textbf{116} (1--4), 321-343.

Howarth, R.W., Sharpley, A., Walker, D. 2002. Sources of nutrient
pollution to coastal waters in the United States: Implications for
achieving coastal water quality goals. \emph{Estuaries} , \textbf{25}
(4), 656-676.

Huang, J.-C., Lee, T.-Y., Kao, S.-J., Hsu, S.-C., Lin, H.-J., Peng,
T.-R. 2012. Land use effect and hydrological control on nitrate yield in
subtropical mountainous watersheds. \emph{Hydrology and Earth System
Sciences} , \textbf{16} (3), 699-714.

Hyer, K.E., Denver, J.M., Langland, M.J., Webber, J.S., Bohlke, J.,
Hively, W.D., Clune, J.W. 2016. Spatial and temporal variation of stream
chemistry associated with contrasting geology and land-use patterns in
the Chesapeake Bay watershed---Summary of results from Smith Creek,
Virginia; Upper Chester River, Maryland; Conewago Creek, Pennsylvania;
and Difficult Run, Virginia, 2010--2013. US Geological Survey.
2328-0328.

Jacobs, S.R., Weeser, B., Guzha, A.C., Rufino, M.C., Butterbach\selectlanguage{english}-Bahl,
K., Windhorst, D., Breuer, L. 2018. Using High\selectlanguage{english}-Resolution Data to Assess
Land Use Impact on Nitrate Dynamics in East African Tropical Montane
Catchments. \emph{Water Resources Research} , \textbf{54} (3),
1812-1830.

Jawitz, J.W., Mitchell, J. 2011. Temporal inequality in catchment
discharge and solute export. \emph{Water Resources Research}
,\textbf{47} (10).

Jiang, R., Woli, K.P., Kuramochi, K., Hayakawa, A., Shimizu, M., Hatano,
R. 2010. Hydrological process controls on nitrogen export during storm
events in an agricultural watershed. \emph{Soil Science \& Plant
Nutrition} , \textbf{56} (1), 72-85.

Jurgensen, M.F., Page-Dumroese, D.S., Brown, R.E., Tirocke, J.M.,
Miller, C.A., Pickens, J.B., Wang, M. 2017. Estimating carbon and
nitrogen pools in a forest soil: influence of soil bulk density methods
and rock content. \emph{Soil Science Society of America Journal}
,\textbf{81} (6), 1689-1696.

Kaushal, S.S., Belt, K.T. 2012. The urban watershed continuum: evolving
spatial and temporal dimensions. \emph{Urban Ecosystems} ,\textbf{15}
(2), 409-435.

Kaushal, S.S., Delaney-Newcomb, K., Findlay, S.E.G., Newcomer, T.A.,
Duan, S.W., Pennino, M.J., Sivirichi, G.M., Sides-Raley, A.M.,
Walbridge, M.R., Belt, K.T. 2014. Longitudinal patterns in carbon and
nitrogen fluxes and stream metabolism along an urban watershed
continuum. \emph{Biogeochemistry} , \textbf{121} (1), 23-44.

Kaushal, S.S., Groffman, P.M., Band, L.E., Shields, C.A., Morgan, R.P.,
Palmer, M.A., Belt, K.T., Swan, C.M., Findlay, S.E., Fisher, G.T. 2008.
Interaction between urbanization and climate variability amplifies
watershed nitrate export in Maryland. \emph{Environmental Science \&
Technology} , \textbf{42} (16), 5872-5878.

Kaushal, S.S., Likens, G.E., Utz, R.M., Pace, M.L., Grese, M., Yepsen,
M. 2013. Increased river alkalinization in the Eastern
US.\emph{Environmental science \& technology} , \textbf{47} (18),
10302-10311.

King, K.W., Williams, M.R., Macrae, M.L., Fausey, N.R., Frankenberger,
J., Smith, D.R., Kleinman, P.J.A., Brown, L.C. 2015. Phosphorus
Transport in Agricultural Subsurface Drainage: A Review. \emph{Journal
of Environmental Quality} , \textbf{44} (2), 467-485.

Kolbe, T., De Dreuzy, J.-R., Abbott, B.W., Aquilina, L., Babey, T.,
Green, C.T., Fleckenstein, J.H., Labasque, T., Laverman, A.M., Marcais,
J. 2019. Stratification of reactivity determines nitrate removal in
groundwater. \emph{Proceedings of the National Academy of Sciences}
,\textbf{116} (7), 2494-2499.

LaBaugh, J.W., Harte, P.T., Shapiro, A.M., Hsieh, P.A., Johnson, C.D.,
Goode, D.J., Wood, W.W., Buso, D.C., Likens, G.E., Winter, T.C. 2013.
Physical, chemical, and isotopic data from groundwater in the watershed
of Mirror Lake, and in the vicinity of Hubbard Brook, near West
Thornton, New Hampshire, 1983 to 1997. US Geological Survey. 2331-1258.

Le Moal, M., Gascuel-Odoux, C., Menesguen, A., Souchon, Y., Etrillard,
C., Levain, A., Moatar, F., Pannard, A., Souchu, P., Lefebvre, A. 2019.
Eutrophication: a new wine in an old bottle? \emph{Science of the Total
Environment} , \textbf{651} , 1-11.

Lerner, D.N., Yang, Y.S., Barrett, M.H., Tellam, J.H. 1999. Loadings of
non-agricultural nitrogen in urban groundwater. \emph{Impacts of Urban
Growth on Surface Water and Groundwater Quality} (259), 117-123.

Li, L. 2019. Watershed Reactive Transport. in: \emph{Reactive transport
modeling in natural and engineered system} , Vol. 85, Society of
Mineralogy and Geochemistry. Reviews in Mineralogy \& Geochemistry.

Li, L., Bao, C., Sullivan, P.L., Brantley, S., Shi, Y., Duffy, C. 2017.
Understanding watershed hydrogeochemistry: 2. Synchronized hydrological
and geochemical processes drive stream chemostatic behavior. \emph{Water
Resources Research} , \textbf{53} (3), 2346-2367.

Luthy, R.G., Sedlak, D.L., Plumlee, M.H., Austin, D., Resh, V.H. 2015.
Wastewater-effluent-dominated streams as ecosystem-management tools in a
drier climate. \emph{Frontiers in Ecology and the Environment}
,\textbf{13} (9), 477-485.

Machesky, M.L., Holm, T.R., Slowikowski, J.A. 2010. Phosphorus
speciation in stream bed sediments from an agricultural watershed:
Solid-phase associations and sorption behavior. \emph{Aquatic
geochemistry} , \textbf{16} (4), 639-662.

Marti, E., Riera, J.L., Sabater, F. 2010. Effects of Wastewater
Treatment Plants on Stream Nutrient Dynamics Under Water Scarcity
Conditions. \emph{Water Scarcity in the Mediterranean: Perspectives
under Global Change} , \textbf{8} , 173-195.

Mayer, B., Boyer, E.W., Goodale, C., Jaworski, N.A., Van Breemen, N.,
Howarth, R.W., Seitzinger, S., Billen, G., Lajtha, K., Nadelhoffer, K.
2002. Sources of nitrate in rivers draining sixteen watersheds in the
northeastern US: Isotopic constraints. \emph{Biogeochemistry}
,\textbf{57} (1), 171-197.

Miller, M.P., Tesoriero, A.J., Hood, K., Terziotti, S., Wolock, D.M.
2017. Estimating Discharge and Nonpoint Source Nitrate Loading to
Streams From Three End-Member Pathways Using High-Frequency Water
Quality Data. \emph{Water Resources Research} , \textbf{53} (12),
10201-10216.

Moatar, F., Abbott, B.W., Minaudo, C., Curie, F., Pinay, G. 2017.
Elemental properties, hydrology, and biology interact to shape
concentration\selectlanguage{english}-discharge curves for carbon, nutrients, sediment, and
major ions. \emph{Water Resources Research} , \textbf{53} (2),
1270-1287.

Montross, G.G., McGlynn, B.L., Montross, S.N., Gardner, K.K. 2013.
Nitrogen production from geochemical weathering of rocks in southwest
Montana, USA. \emph{Journal of Geophysical Research: Biogeosciences}
,\textbf{118} (3), 1068-1078.

Mulholland, P.J., Helton, A.M., Poole, G.C., Hall, R.O., Hamilton, S.K.,
Peterson, B.J., Tank, J.L., Ashkenas, L.R., Cooper, L.W., Dahm, C.N.,
Dodds, W.K., Findlay, S.E.G., Gregory, S.V., Grimm, N.B., Johnson, S.L.,
McDowell, W.H., Meyer, J.L., Valett, H.M., Webster, J.R., Arango, C.P.,
Beaulieu, J.J., Bernot, M.J., Burgin, A.J., Crenshaw, C.L., Johnson,
L.T., Niederlehner, B.R., O'Brien, J.M., Potter, J.D., Sheibley, R.W.,
Sobota, D.J., Thomas, S.M. 2008. Stream denitrification across biomes
and its response to anthropogenic nitrate loading. \emph{Nature}
,\textbf{452} (7184), 202-U46.

Musolff, A., Fleckenstein, J.H., Rao, P.S.C., Jawitz, J.W. 2017.
Emergent archetype patterns of coupled hydrologic and biogeochemical
responses in catchments. \emph{Geophysical Research Letters}
,\textbf{44} (9), 4143-4151.

Musolff, A., Schmidt, C., Selle, B., Fleckenstein, J.H. 2015. Catchment
controls on solute export. \emph{Advances in Water Resources}
,\textbf{86} , 133-146.

Newcomer Johnson, T.A., Kaushal, S.S., Mayer, P.M., Grese, M.M. 2014.
Effects of stormwater management and stream restoration on watershed
nitrogen retention. \emph{Biogeochemistry} , \textbf{121} (1), 81-106.

Pellerin, B.A., Saraceno, J.F., Shanley, J.B., Sebestyen, S.D., Aiken,
G.R., Wollheim, W.M., Bergamaschi, B.A. 2012. Taking the pulse of
snowmelt: in situ sensors reveal seasonal, event and diurnal patterns of
nitrate and dissolved organic matter variability in an upland forest
stream. \emph{Biogeochemistry} , \textbf{108} (1-3), 183-198.

Pellerin, B.A., Wollheim, W.M., Feng, X., Vorosmarty, C.J. 2008. The
application of electrical conductivity as a tracer for hydrograph
separation in urban catchments. \emph{Hydrological Processes: An
International Journal} , \textbf{22} (12), 1810-1818.

Pennino, M.J., Kaushal, S.S., Mayer, P.M., Utz, R.M., Cooper, C.A. 2016.
Stream restoration and sewers impact sources and fluxes of water,
carbon, and nutrients in urban watersheds. \emph{Hydrology and Earth
System Sciences} , \textbf{20} (8), 3419-3439.

Perrin, A.-S., Probst, A., Probst, J.-L. 2008. Impact of nitrogenous
fertilizers on carbonate dissolution in small agricultural catchments:
Implications for weathering CO2 uptake at regional and global
scales.\emph{Geochimica et Cosmochimica Acta} , \textbf{72} (13),
3105-3123.

Pickett, S.T., Cadenasso, M.L., Grove, J.M., Boone, C.G., Groffman,
P.M., Irwin, E., Kaushal, S.S., Marshall, V., McGrath, B.P., Nilon, C.H.
2011. Urban ecological systems: Scientific foundations and a decade of
progress. \emph{Journal of Environmental Management} , \textbf{92} (3),
331-362.

Poor, C.J., McDonnell, J.J. 2007. The effects of land use on stream
nitrate dynamics. \emph{Journal of Hydrology} , \textbf{332} (1-2),
54-68.

Prein, A.F., Rasmussen, R.M., Ikeda, K., Liu, C., Clark, M.P., Holland,
G.J. 2017. The future intensification of hourly precipitation
extremes.\emph{Nature Climate Change} , \textbf{7} (1), 48-52.

Raffensperger, J.P., Baker, A.C., Blomquist, J.D., Hopple, J.A.
2017.\emph{Optimal hydrograph separation using a recursive digital
filter constrained by chemical mass balance, with application to
selected Chesapeake Bay watersheds} . US Geological Survey.

Royer, T.V., David, M.B., Gentry, L.E. 2006. Timing of riverine export
of nitrate and phosphorus from agricultural watersheds in Illinois:
Implications for reducing nutrient loading to the Mississippi
River.\emph{Environmental Science \& Technology} , \textbf{40} (13),
4126-4131.

Sebestyen, S.D., Boyer, E.W., Shanley, J.B., Kendall, C., Doctor, D.H.,
Aiken, G.R., Ohte, N. 2008. Sources, transformations, and hydrological
processes that control stream nitrate and dissolved organic matter
concentrations during snowmelt in an upland forest. \emph{Water
Resources Research} , \textbf{44} (12).

Sebilo, M., Mayer, B., Nicolardot, B., Pinay, G., Mariotti, A. 2013.
Long-term fate of nitrate fertilizer in agricultural
soils.\emph{Proceedings of the National Academy of Sciences of the
United States of America} , \textbf{110} (45), 18185-18189.

Seibert, J., Grabs, T., Kohler, S., Laudon, H., Winterdahl, M., Bishop,
K. 2009. Linking soil- and stream-water chemistry based on a Riparian
Flow-Concentration Integration Model. \emph{Hydrol. Earth Syst. Sci.}
,\textbf{13} (12), 2287-2297.

Shepherd, M., Ghani, A., Rajendram, G., Carlson, B., Pirie, M. 2015.
Soil total nitrogen concentration explains variation in pasture response
to spring nitrogen fertiliser across a single farm. \emph{Nutrient
Cycling in Agroecosystems} , \textbf{101} (3), 377-390.

Shi, Y., Davis, K.J., Zhang, F., Duffy, C.J. 2014. Evaluation of the
parameter sensitivities of a coupled land surface hydrologic model at a
critical zone observatory. \emph{Journal of Hydrometeorology}
,\textbf{15} (1), 279-299.

Sickman, J.O., Leydecker, A., Chang, C.C., Kendall, C., Melack, J.M.,
Lucero, D.M., Schimel, J. 2003. Mechanisms underlying export of N from
high-elevation catchments during seasonal
transitions.\emph{Biogeochemistry} , \textbf{64} (1), 1-24.

Spahr, N.E., Dubrovsky, N.M., Gronberg, J.M., Franke, O., Wolock, D.M.
2010. Nitrate loads and concentrations in surface-water base flow and
shallow groundwater for selected basins in the United States, water
years 1990-2006. U. S. Geological Survey.

Sudduth, E.B., Perakis, S.S., Bernhardt, E.S. 2013. Nitrate in
watersheds: Straight from soils to streams? \emph{Journal of Geophysical
Research-Biogeosciences} , \textbf{118} (1), 291-302.

Tedesco, L., Pascual, D., Shrake, L., Casey, L., Hall, B., Vidon, P.,
Hernly, F., Barr, R., Ulmer, J., Pershing, D. 2005. Eagle Creek
watershed management plan: An integrated approach to improved water
quality. \emph{Eagle creek watershed alliance, CEES publication}
,\textbf{7} , 182.

Tesoriero, A.J., Terziotti, S., Abrams, D.B. 2015. Predicting Redox
Conditions in Groundwater at a Regional Scale. \emph{Environmental
Science \& Technology} , \textbf{49} (16), 9657-9664.

Thompson, S., Basu, N., Lascurain, J., Aubeneau, A., Rao, P. 2011.
Relative dominance of hydrologic versus biogeochemical factors on solute
export across impact gradients. \emph{Water Resources Research}
,\textbf{47} (10).

Van Cappellen, P., Maavara, T. 2016. Rivers in the Anthropocene: Global
scale modifications of riverine nutrient fluxes by
damming.\emph{Ecohydrology \& Hydrobiology} , \textbf{16} (2), 106-111.

Van Meter, K., Basu, N., Van Cappellen, P. 2017. Two centuries of
nitrogen dynamics: Legacy sources and sinks in the Mississippi and
Susquehanna River Basins. \emph{Global Biogeochemical Cycles}
,\textbf{31} (1), 2-23.

Van Meter, K.J., Basu, N.B. 2015. Catchment legacies and time lags: A
parsimonious watershed model to predict the effects of legacy storage on
nitrogen export. \emph{PLoS One} , \textbf{10} (5), e0125971.

Van Meter, K.J., Basu, N.B., Veenstra, J.J., Burras, C.L. 2016. The
nitrogen legacy: emerging evidence of nitrogen accumulation in
anthropogenic landscapes. \emph{Environmental Research Letters}
,\textbf{11} (3), 035014.

Vidon, P.G.F., Hill, A.R. 2004. Landscape controls on nitrate removal in
stream riparian zones. \emph{Water Resources Research} , \textbf{40}
(3).

Weiler, M., McDonnell, J.R.J. 2006. Testing nutrient flushing hypotheses
at the hillslope scale: A virtual experiment approach. \emph{J. Hydrol.}
, \textbf{319} (1-4), 339-356.

Weitzman, J.N., Kaye, J.P. 2018. Nitrogen Budget and Topographic
Controls on Nitrous Oxide in a Shale-Based Watershed. \emph{Journal of
Geophysical Research: Biogeosciences} , \textbf{123} (6), 1888-1908.

Wollheim, W.M., Pellerin, B.A., Vorosmarty, C.J., Hopkinson, C.S. 2005.
N retention in urbanizing headwater catchments. \emph{Ecosystems}
,\textbf{8} (8), 871-884.

Woo, D.K., Kumar, P. 2019. Impacts of Subsurface Tile Drainage on
Age---Concentration Dynamics of Inorganic Nitrogen in Soil. \emph{Water
Resources Research} , \textbf{55} (2), 1470-1489.

Wosten, J., Pachepsky, Y.A., Rawls, W. 2001. Pedotransfer functions:
bridging the gap between available basic soil data and missing soil
hydraulic characteristics. \emph{Journal of hydrology} ,\textbf{251}
(3-4), 123-150.

Yan, Z.F., Liu, C.X., Todd-Brown, K.E., Liu, Y.Y., Bond-Lamberty, B.,
Bailey, V.L. 2016. Pore-scale investigation on the response of
heterotrophic respiration to moisture conditions in heterogeneous
soils.\emph{Biogeochemistry} , \textbf{131} (1-2), 121-134.

Zamanian, K., Zarebanadkouki, M., Kuzyakov, Y. 2018. Nitrogen
fertilization raises CO2 efflux from inorganic carbon: A global
assessment. \emph{Global Change Biology} , \textbf{24} (7), 2810-2817.

Zhang, Q. 2018. Synthesis of nutrient and sediment export patterns in
the Chesapeake Bay watershed: Complex and non-stationary
concentration-discharge relationships. \emph{Science of the Total
Environment} , \textbf{618} , 1268-1283.

Zhi, W., Li, L., Dong, W.M., Brown, W., Kaye, J., Steefel, C., Williams,
K.H. 2019. Distinct Source Water Chemistry Shapes Contrasting
Concentration-Discharge Patterns. \emph{Water Resources Research}
,\textbf{55} (5), 4233-4251.

Zhou, T., Shi, P., Hui, D., Luo, Y. 2009. Global pattern of temperature
sensitivity of soil heterotrophic respiration (Q10) and its implications
for carbon\selectlanguage{english}-climate feedback. \emph{Journal of Geophysical Research:
Biogeosciences} , \textbf{114} (G2).

\selectlanguage{english}
\FloatBarrier
\end{document}

