Where the past meets the present: connecting nitrogen from watersheds to streams through groundwater flowpaths

Groundwater discharge to streams is a nonpoint source of nitrogen (N) that confounds N mitigation efforts and represents a significant portion of the annual N loading to watersheds. However, we lack an understanding of where and how much groundwater N enters streams and watersheds. Nitrogen concentrations at the end of groundwater flowpaths are the culmination of biogeochemical and physical processes from the contributing land area where groundwater recharges, within the aquifer system, and in the near-stream riparian area where groundwater discharges to streams. Our research objectives were to quantify the spatial distribution of N concentrations at groundwater discharges throughout a mixed land-use watershed and to evaluate how relationships among contributing and riparian land cover, modeled aquifer characteristics, and groundwater discharge biogeochemistry explain the spatial variation in groundwater discharge N concentrations. We accomplished this by integrating high-resolution thermal infrared surveys to locate groundwater discharge, biogeochemical sampling of groundwater, and a particle tracking model that links groundwater discharge locations to their contributing area land cover. Groundwater N loading from groundwater discharges within the watershed varied substantially between and within streambank groundwater discharge features. Groundwater nitrate concentrations were spatially heterogeneous ranging from below 0.03–11.45 mg-N/L, varying up to 20-fold within meters. When combined with the particle tracking model results and land cover metrics, we found that groundwater discharge nitrate concentrations were best predicted by a linear mixed-effect model that explained over 60% of the variation in nitrate concentrations, including aquifer chemistry (dissolved oxygen, Cl−, SO4 2−), riparian area forested land cover, and modeled physical aquifer characteristics (discharge, Euclidean distance). Our work highlights the significant spatial variability in groundwater discharge nitrate concentrations within mixed land-use watersheds and the need to understand groundwater N processing across the many spatiotemporal scales within groundwater cycling.


Introduction
Nitrogen (N) applied to land surfaces infiltrates with groundwater recharge and has substantially increased N concentrations in aquifers over the last century (Galloway and Cowling 2002, Puckett et al 2011, Houlton and Morford 2014).Groundwater N discharge to streams represents a widespread nonpoint source of N that can be significantly delayed from terrestrial applications based on varying groundwater transit times (Tesoriero et al 2013, Stets et al 2020) causing groundwater legacies to confound surface water quality management (Rosenberry et al 2016, Liu et al 2017, Vero et al 2018).Although groundwater N is known to be a substantial component of N loading to streams (Johnson and Stets 2020, Wherry et al 2021), we lack the ability to predict the spatial distribution of groundwater N loading to streams, hindering mitigation efforts (Sanford andPope 2013, Van Meter et al 2018).
Groundwater N loading to streams depends on complex biogeochemical interactions across space and time (Rivett et al 2008).Nitrogen deposition on land surfaces is often heterogeneous, and N may undergo multiple transformations before entering surface water (Kaushal et al 2011, Kolbe et al 2019).Initially, some N is retained by vegetation and soil through biological and physical processes (Fowler et al 2013).Subsequently, the fate of N percolating to aquifers is controlled by groundwater residence times, redox gradients, and electron donor availability along groundwater flowpaths (Ocampo et al 2006, Kolbe et al 2019, Tesoriero et al 2021, Henri and Harter 2022).As groundwater interfaces with surface water at the end of flowpaths, organic carbon within the river corridor may promote further N processing through denitrification, the microbially mediated reduction of nitrate to nitrous oxide or nitrogen gas (Lutz et al 2020).Riparian buffers and reducing N application mitigate loading to surface waters (Ranalli and Macalady 2010).However, the effects of best management practices are not always realized downstream (Chang et al 2021, Martin et al 2021) and in some cases, riparian buffers can create preferential flowpaths that allow N to bypass management practices (Hester and Fox 2020).
Groundwater N is transported from the land surface along groundwater flowpaths to surface waters.Groundwater flowpaths form in response to recharge patterns, surface and bedrock topography, and hydraulic conductivity (Winter et al 1998).However, identifying where flowpaths interface with surface water throughout a stream network is challenging.Field mapping of groundwater discharge at broad scales is difficult due to access to continuous reaches and time-intensive field methods (Harvey andWagener 2000, Rosenberry et al 2021).Predicting groundwater discharge locations based on subsurface characteristics is limited by the spatial resolution and accuracy of existing datasets (e.g.SSURGO, Nauman et al 2014), and is further complicated by land use (Scanlon et al 2005), which can disconnect groundwater from surface water.
Groundwater models can be used to predict where groundwater discharges, but field validation is challenging because of the significant heterogeneity in how groundwater is expressed along streams in terms of discharge, spatial distribution, and lateral extent (Briggs et al 2021, Barclay et al 2022).Thermal infrared (TIR) cameras can be used to map groundwater discharge from the point-to-reach scale (Dugdale et al 2015, Hare et al 2017, Sullivan et al 2021) and have recently been used to evaluate modular finite difference flow model (MODFLOW)predicted discharge locations (Barclay et al 2022).In this study, we integrate high-resolution groundwater discharge TIR surveys and biogeochemical sampling with a particle tracking model to spatially connect groundwater discharge to its contributing land area and physical aquifer characteristics.Reducing the uncertainty of groundwater N inputs represents a critical knowledge gap for building more resilient riverscapes, particularly as streams face further human-induced changes (Wherry et al 2021).
We hypothesized that spatial variation in groundwater NO 3 − concentrations is a function of (1) land cover contributing to groundwater flowpaths, (2) physical characteristics that describe water movement through the aquifer (3) groundwater chemistry representing reactivity along flowpaths, and (4) riparian characteristics controlling biogeochemical processing along and at the end of groundwater flowpaths (figure 1).The goal of this study is to quantify the spatial patterns of N in groundwater discharges throughout a mixed land-use stream network.Our objectives were to (1) map the spatial variation of groundwater discharge N concentrations, and (2) evaluate the relationships between groundwater discharge N concentrations and groundwater chemistry, modeled aquifer characteristics, and contributing area and riparian land cover.

Study area
The Farmington River watershed (1572 km 2 ) is a fifth-order tributary of the Connecticut River located in Connecticut and Massachusetts, USA.It is predominantly forested (67%) with development (18%), agriculture (3%), and grassland (1%) primarily along the mainstem and forested wetlands (8%) and open water (3%) throughout the watershed (Dewitz and U.S. Geological Survey 2021, figure 2).Bedrock geology consists of New England crystalline rock, Mesozoic sandstone, and Newark Supergroup basalt, which is overlain by coarse-grained, stratified glacial deposits in the lower portion of the watershed while the upper portion consists of fine-grained, unstratified glacial till (Olcott 1995).Bedrock depth along the mainstem is 23.2 m ± 5.5 m (Jackson et al 2023) and mean saturated hydraulic conductivity is 9.3 m d −1 (Briggs et al 2021).The daily mean discharge is 30.9 m 3 s −1 (USGS NWIS, gage 01189995, 2010-2022).

Field surveys of groundwater discharge
We used handheld TIR cameras (FLIR, Wilsonville, OR) to identify areas of groundwater discharge along streambanks (Deitchman and Loheide 2009,  Briggs and Hare 2018, Barclay et al 2022).Following Barclay et al (2022) and Briggs et al (2021), we surveyed streams via wading and canoeing to map groundwater discharges at the meter scale along 36 stream reaches totaling over 60 km from 2017 and 2019-2020.Surveys were completed during baseflow (July-November, <30.9 m 3 s −1 ) when surface water temperatures were above 15 • C and warmer than discharging groundwater (figure 2).Groundwater discharges were characterized using GPS coordinates, TIR, and direct temperature and lateral extent measurements.We identified over 300 unique groundwater discharges ranging from 1 m to over 200 m in lateral extent (Moore et al 2020, Barclay et al 2022).

Sample collection
We quantified chemical and spatial variability of groundwater chemistry by sampling mapped groundwater discharge locations (2017,(2019)(2020)(2021).Because groundwater flux rates within discharge features can vary significantly (Haynes et al 2023), we collected our samples from the coldest point assuming this was the most representative sample of groundwater flow through the streambank.We used pushpoint samplers (MHE Products, East Tawas, MI) inserted 20 cm, or until resistance, into the sediment (n = 363).We sampled previously mapped groundwater discharges along 25 km of the 5thorder mainstem and 19 small streams (figure 2).Along the mainstem, we randomly chose locations to sample (n = 107) and sampled as many groundwater discharge locations as possible in smaller streams (n = 123).Groundwater was pulled through the push-point sampler and purged until it ran clear before sample collection.We measured specific conductance, temperature, and dissolved oxygen (DO) in situ using a calibrated YSI-6000 (YSI Inc., Yellow Springs, OH).Groundwater was collected to analyze nitrate (NO 3 − ), ammonia (NH 4 + ), total dissolved nitrogen (TDN), chloride (Cl − ), sulfate (SO 4 2− ), and dissolved organic carbon (DOC) in acid-washed and field-rinsed HDPE bottles.We quantified dissolved nitrogen gas (N 2 ), a proxy for denitrification (Böhlke et al 2002, Kennedy et al 2009) following Lamberti and Hauer (2017), and nitrous oxide (N 2 O) following Helton et al (2014).Coincident surface water samples were collected (1-2 per sampling day).Additional details on sample treatment and analysis can be found in Moore et al (2023) and supplemental information (SI).

Particle tracking model
Groundwater discharge locations throughout the network were estimated using a MODPATH particle tracking model (Pollock 2012), which tracks hypothetical water particles through an aquifer based on the outputs of a MODFLOW groundwater model (Niswonger et al 2011).We used an existing model (Barclay et al 2020b, model name: RivK_BedK_drn) with 250,000 particles across a 300 m 2 grid.Particles were distributed to starting cells (n = 18,113) based on recharge rate (Reitz et al 2017) and were tracked from the model start cell on the land surface to model end cells (n = 2998) along the network with unique flowpath lines.Modeled aquifer characteristics were generated for each model end cell discharging along the network, including median residence time (years), median Euclidean distance (m), maximum flowpath depth (m), and median groundwater discharge per reach length (m 3 m −1 ).Additional details on model implementation can be found in SI and Barclay et al (2020a), (2020b).

Riparian and contributing area land cover
We evaluated watershed land cover from 2001 to 2019 finding minimal change (forest −1.47%, developed +1.30%, agriculture −0.38%, wetland +0.38%, grassland +0.61%, water −0.43%); therefore, we opted to use the 2019 National Land Cover Dataset (Dewitz and U.S. Geological Survey 2021).We combined NLCD classes into forest (gridcodes 41-43),  developed (21-31), agriculture (81-82), wetland (90-95), grassland (51-74), and water (11-12) (figure 2).Because fertilizer application rates are only available at the county-scale (Byrnes et al 2020), we used percent land cover as a proxy for anthropogenic nitrogen sources.We also accounted for anthropogenic N sources by calculating the density of people on septic.We created a spatial data layer representing where people are most likely to dwell in homes with septic systems, or 'settled land' , by using a 150 m buffer of developed (low density, medium density, and open space) classes of the 1996 C-CAP land cover (NOAA 2022).We paired 'settled land' with the number of people on septic systems reported in the 1990 census (U.S. Census Bureau 1990) by dividing the number of people on septic in each census block by the area of 'settled land' in each census block.The result is a septic density value for 'settled land' in each census block.We estimated percent land cover and septic density for the riparian area using a 300 m buffer radius around each sampling location.
We also estimated percent land cover and septic density for the land area contributing groundwater to each groundwater discharge sampling location (i.e.contributing area).We discretized the area of each 2019 land cover class and septic density for each MODPATH model start cell.For each model end cell along the stream network, we traced all particles transported to that end cell back to their start cells.The group of start cells connected to a model end cell make up that end cell's contributing area.We then calculated the weighted average land cover based on the number of particles connecting each starting cell to the end cell.Thus, we connected sampling locations to the weighted average contributing area land cover supplying recharge.The contributing area to a groundwater discharge location includes all land surface cells connected by groundwater flow paths, regardless of their hydrologic connection via surface topography.Spatial analyses were performed using ArcMap (ESRI 2021 Inc., Redlands, CA).

−
We conducted five linear mixed-effects model selections based on four sets of predictor variables that represent the major drivers of discharging NO 3 − concentrations: (1) contributing area weighted percent forest, wetland, and agriculture, and density of people on septic, (2) simulated residence time, Euclidean distance traveled, flowpath depth, and groundwater discharge rate, (3) measured DO, SO 4 2− , Cl − , and DOC concentrations in groundwater (representative of aquifer biogeochemistry), and (4) riparian area percent forest, wetland, and agriculture, and density of people on septic.We used sampling year as a random effect to account for differences in the annual flow regime because 2021 had consistently higher river discharge due to three hurricanes.Percent forest and developed land cover were highly correlated (Pearson's r = −0.92for contributing area and −0.53 for riparian area) therefore we opted to exclude developed land cover.Model selection was performed through an exhaustive search of all model variables.Best fit models were selected as those with the lowest Bayesian information criterion and highest model weight (table S1).We also conducted model selection for the combined categories of variables including only those variables within the best-fit models of the first four model selections.In summary, we conducted four individual sets of models (contributing area, aquifer characteristics, groundwater biogeochemistry, and riparian buffer) and one combined model using the best fit model from each individual model set.Variables were log 10 transformed to meet assumptions of normality (Shapiro-Wilks, p < 0.05).Statistical analyses were performed using R 4.1.2(R Core Team 2021).

Spatial patterns of N at groundwater discharges
Groundwater N was heterogeneous throughout the watershed (figure 3).Groundwater collected only meters apart often varied substantially; for example, at river km 22, NO 3 − decreased from 7.42 to 0.37 mg-N/L within 30 m. TDN and NO 3 -were highly correlated across the watershed (rho = 0.937, p < 0.05, figure S1), and NO 3 − constituted over 90% of TDN, indicating that most of the N sampled was NO 3 − .We generally observed low NH 4 + along the 25 km reach, yet there were rare locations with NH 4 + up to 2 mg l −1 , and NH 4 + was negatively correlated with NO 3 − (rho = −0.28,p < 0.05).Small stream concentrations and ranges of NO 3 − and NH 4 + in streams were lower than the mainstem (table 1), reflecting the shift from predominantly forest to mixed land use.
Denitrification end-products also varied considerably, indicating heterogeneity in N processing.Mainstem N 2 O was positively correlated with NO 3 − (rho = 0.58, p < 0.05) whereas excess N 2 , our surrogate measure of denitrification, was not (rho = −0.03,p > 0.05).Small stream N 2 O and N 2 groundwater concentrations were significantly lower than mainstem concentrations and N 2 O was positively correlated with NO 3 − (rho = 0.66, p < 0.05), whereas N 2 was not (rho = 0.02, p > 0.05).Only 60% of N 2 samples contained excess N 2 from complete denitrification.DO varied considerably across the 25 km reach, indicating significant redox gradient heterogeneity within the aquifer (Tesoriero et al 2021).A laterally extensive groundwater face near river km 11 contained DO ranging from 1.35 to 6.07 mg l −1 within 150 m (figure 3).Mainstem and small stream DO were not significantly different but were positively correlated with NO 3 − (mainstem rho = 0.49, p < 0.05; small stream rho = 0.28, p < 0.05).
We also collected other geochemical variables that may influence groundwater N cycling.DOC concentrations were generally low (median: 0.85 mg l −1 ) and were not correlated with NO 3 − along the mainstem (rho = −0.19,p > 0.05) or small streams (rho = −0.02,p > 0.05).Mainstem DOC was significantly higher than headwaters (table 1).Sulfate concentrations were significantly higher along the mainstem than small streams, and SO 4 2− was positively correlated with NO 3 − along the mainstem (rho = 0.32, p < 0.05) and small streams (rho = 0.30, p < 0.05).Chloride concentrations were also heterogeneous ranging from 2.04 to 2337.84 mg l −1 along the mainstem and 0.03-123.31mg l −1 in small streams.Sediment grab samples show that over 90% of the sediment was sand or larger grain size (Moore et al 2023).Organic matter ranged from 0.2 to 13.34% along the mainstem and 0.2%-73.91% in small streams.

Particle tracking model
Particle tracking applied to a calibrated groundwater flow model from Barclay et al (2020b) spatially connected our samples to contributing area land cover.Across the entire watershed, modeled groundwater flowpath characteristics were generally consistent with model predictions at sampling locations (SI).MODPATH particles terminated at 242 of the 363 sampling locations.The remaining 121 samples occurred along reaches that were not predicted to receive particles, particularly in 1st and 2nd order streams (figure 4).Modeled flowpath characteristics increased with stream order (figure S2).At sampling locations, groundwater residence times ranged from 0.1 to 406.7 years (median 9.2 years) and the maximum flowpath depth ranged from 0.7 to 183.6 m (median = 96.0m).Euclidean distance to sampling locations ranged from 150.0-3400.8m (median = 1006.1 m) and median groundwater discharge rates ranged from 3 × 10 −3 to 18.8 m 3 m −1 (median = 3.8 m 3 m −1 ).

−
The model that best explained groundwater discharge NO 3 − concentrations included variables that describe riparian land cover, aquifer biogeochemistry, and modeled groundwater discharge and flowpath distance (table 2).Aquifer biogeochemistry and riparian land cover variables were the strongest individual group of predictors of NO 3 − (R 2 = 0.41 and 0.32 respectively, table 2).Concentrations of DO, SO 4 2− , and Cl − were positively correlated with NO 3 − .Riparian and contributing area forest cover were both negatively related to NO 3 − , though contributing area forest cover was not included in the bestperforming model.Interestingly, riparian agricultural land cover did not correlate with groundwater NO 3 − and was not selected through model subsetting (table S1).However, the density of people on septic systems is an influential variable within the riparian area.Aquifer characteristics derived from particle tracking predicted groundwater NO 3 − to increase with increased discharge and decrease with distance traveled but explained relatively little variation (R 2 = 0.15).Contributing area land cover was not a significant parameter within the full model selection indicating that flowpath characteristics and biogeochemistry convolute connections between discharging groundwater NO 3 − patterns and contributing area land cover.
We accounted for flow year as a random effect resulting in increased predictive power in all models except the riparian area model indicating that sampling year flow condition was an influential variable in predicting groundwater NO 3 − concentrations (table 2).Groundwater NO 3 − and NH 4 + were elevated and more heterogeneous during the 2019 low-flow year than the 2021 high-flow year and indicators of denitrification (N 2 O and N 2 ) were also more variable during the 2019 low-flow year (figure S3).

Discussion
Our research shows that groundwater discharge N concentrations are highly heterogeneous both within and among streambank discharge features, and that NO 3 − concentrations are best predicted by (1) aquifer biogeochemistry, (2) riparian area forest cover and density of people on septic, and (3) flowpath distance and groundwater discharge rate.Contrary to our expectation, we found that contributing area land cover was not a major predictor of groundwater NO 3 − .This contrasts prior studies that have directly linked surface land cover to groundwater well N concentrations (Lockhart et al 2013), baseflow N concentrations in streams (Johnson and Stets 2020), and groundwater discharge N concentrations (Cole et al 2006, Shabaga andHill 2010).These studies were conducted along comparatively short riparian flowpaths, where groundwater N concentrations are currently elevated due to contemporary N application.Within our study watershed attempts to link contributing area land cover to groundwater discharge N concentrations were likely complicated by heterogeneous N inputs from mixed land uses and differential N processing along flowpaths.

Farmington River mainstem
Small streams NO3 Throughout the watershed, groundwater well NO 3 -concentrations ranged from non-detect to 9.8 mg-N/L (median = 0.83 mg-N/L, n = 90, 1953-2017, USGS NWIS) and are consistent with our groundwater discharge NO 3 − data (table 1).Groundwater well NO 3 − concentrations show little change over time (figure S4) indicating that while our study watershed contains historic agriculture (Jenkins 1925) groundwater N is not consistently elevated.Spatial variation in N application, including lawn and agricultural fertilizers and septic drainage, may drive areas of high groundwater NO 3 − .For example, high groundwater NO 3 − concentrations were found along the mainstem surrounded by golf courses and agricultural fields.
The extensive chemical heterogeneity in discharging groundwater (figure 3) suggests high uncertainty in predicting groundwater N loading to streams across the watershed.For example, we located approximately 4000 m of streambank groundwater discharge along the mainstem that extended approximately 0.5 m vertically up the streambank (Haynes et al 2023).When we use the 5th (0.078 mg-N/L) and 95th (5.12 mg-N/L) percentile of groundwater discharge NO 3 − concentrations and estimate groundwater flux to be 0.5 m d −1 (Haynes et al 2023), groundwater NO 3 − load ranges from 6740 to 442,360 kg d −1 .Estimating surface water NO 3 − load from the mainstem using our 5th (0.78 mg-N/L) and 95th (1.92 mg-N/L) percentile surface water NO 3 − concentrations and mean daily river discharge (30.9 cms) results in a daily surface water load of 208-5,125 kg d −1 of NO 3 − .These estimates indicate that groundwater is the primary source of N loading  Groundwater discharge N also varied significantly between sampling years where NO 3 − concentrations were lower during the high versus low (2019, 2020) flow years, suggesting that annual variability in groundwater discharge N must also be considered when budgeting N loading from groundwater (figure S3).
Connecting individual groundwater flowpaths from their point of discharge back to surface land use is impossible with empirical methods; therefore, we used a novel integration of field-collected samples and flowpath modeling to connect discharging groundwater chemistry back to contributing land areas.Like Barclay et al (2022) we were able to accurately predict relative groundwater discharge patterns along 3rd and 5th-order streams but struggled to predict discharge in smaller streams (figure 4).Of the 363 groundwater samples collect only 66% of them were predicted to receive particles by the MODPATH model limiting our understanding of how contributing area land cover affects groundwater N concentrations at the river network scale.Modeled flowpath characteristics varied by stream order with 3rd and 5th receiving higher discharge from longer and deeper flowpaths (figure S2).We related predicted locations of groundwater discharge back to their land cover source.However, it did not allow us to generate concrete conclusions about the fate of N from the contributing area, particularly in small streams.This may be due to model resolution; future modeling efforts may consider finer, or even flexible, spatial resolutions to better predict groundwater discharge in small streams.Equally plausible is that reactivity along the groundwater flowpaths could be obscuring the relationship between discharge chemistry and contributing land area (Rivett et al 2008).
Rates of N removal along groundwater flowpaths are typically elevated upon infiltration and at the terminal ends of flowpaths where organic carbon availability is often high (Gorski et al 2022).N removal rates can also be high within aquifers if conditions are appropriate for denitrification (Rivett et al 2008, Kolbe et al 2019, Henri and Harter 2022).Our groundwater discharge samples were low in DOC due to a lack of calcareous bedrock (Olcott 1995), and we found low sediment organic matter at sampling locations, limiting carbon-based denitrification.We found elevated SO 4 2− within discharging groundwater suggesting potential for sulfur-fueled denitrification within the aquifer (Megonigal et al 2003, Ben Maamar et al 2015).We also found evidence of denitrification along groundwater flowpaths (figure 3 − and N 2 O in discharging groundwater.Most of our groundwater DO samples were above the 2 mg l −1 anoxic threshold for denitrification (Rivett et al 2008).High groundwater DO concentrations support our finding that only 60% of samples contained excess N 2 indicating conditions favorable for incomplete rather than complete denitrification.
Low DOC, sub-oxic groundwater, and short residence times within the watershed likely created heterogeneous patterns of N removal and constrained reduction of NO 3 − .However, our physical methods could have biased this interpretation; we identified preferential discharges, which by definition have relatively high flux rates that reduce N attenuation potential (Ocampo et al 2006, Briggs andHare 2018).We did not sample more diffuse groundwater discharge, which may have lower groundwater flux rates but longer residence times and higher rates of biogeochemical reactivity.Even so, our results are consistent with Tesoriero et al (2013) who suggested that watersheds in New England are less vulnerable to legacy groundwater N contamination than watersheds in the agricultural Midwest and urbanized east coast.Yet, heterogenous N concentrations in groundwater and apparent N attenuation within the aquifer suggest localized areas of groundwater N loading to streams are important locations for management (Mullaney and Schwarz 2013).Overall, the Farmington River watershed does not indicate watershed-wide groundwater N pollution and elevated N concentrations within the watershed are likely sourced from local and contemporary anthropogenic sources.Perhaps in New England watersheds, site (and even subsite) level understanding of groundwater N loading is particularly important when making management decisions.Our results indicate that nonpoint source N pollution from groundwater is driven by aquifer biogeochemistry and riparian forests within the river corridor which suggests that riparian buffers are helping to reduce, not only surface water N loading, but also N loading from groundwater sources (Mayer et al 2007).However, we also found elevated N concentrations in groundwater coinciding with little to no N attenuation which also suggests that some groundwater flowpaths are not effectively reducing N within riparian buffers (Hester and Fox 2020).Groundwater N loading is lagged in time behind surface water N loading and will continue to confound surface water quality management from months to decades (Van Meter and Basu 2017).We suggest that watershed managers continue to mitigate N loading through wide riparian buffers and to be cautious when extrapolating nonpoint source groundwater N loading due to significant heterogeneity in where, and how much, groundwater N enters river networks.

Conclusion
Nitrogen concentrations in groundwater discharging to streams-along individual expressions of groundwater discharge, single stream reaches, and across the watershed-were highly heterogeneous, varying 20-fold along the mainstem.Although in this mixed land use watershed, we expected N application at the recharging land area to be a dominant predictor of the variation in N concentration in discharging groundwater, we found instead that the combined effects of groundwater chemistry, flowpath distance and flux rate, and riparian land cover best predicted N concentrations.This suggests that the predictability of N loading from land cover is obscured by reactivity along groundwater flow paths; interestingly, relatively high and variable DO and SO 4 2− paired with low DOC concentrations indicated NO 3 − production through nitrification, and removal through sulfur-based denitrification and incomplete denitrification.Our results show that in mixed land use watersheds groundwater is not a well-mixed system and that N loading from groundwater is highly heterogenous across space further complicating our understanding of how land cover influences groundwater biogeochemistry.The consequences of heterogeneous N delivery from groundwater discharge are particularly important for extrapolating and predicting how groundwater-delivered N loads contribute to watershed export.Improved understanding of the fate of nonpoint source loading to watersheds will benefit efforts to mitigate N loading downstream.

Figure 1 .
Figure 1.Conceptual model linking zones along groundwater flowpaths where controls of nitrogen removal are expected to vary-upon recharge in contributing areas, within the aquifer, and riparian area and near-stream sediments-with particle tracking model structure.Gray squares with black points represent particle tracking model cells with groundwater flowpath lines colored by representative contributing area land cover.Blue squares near the stream represent locations of predicted groundwater discharge.Groundwater discharge sampling locations are shown with a white star.

Figure 2 .
Figure 2. The Farmington River watershed.Sampled groundwater discharges are shown as black circles (n = 363).The mainstem 25 km reach is highlighted with a blue box.The right panel shows two groundwater discharge faces with coincident thermal and real-color images where regions of red and white are near air temperature and blue regions are cold groundwater discharging along the river.Images by Eric Moore 8/2022.

Figure 3 .
Figure 3. Longitudinal groundwater discharge sampling along the mainstem Farmington River, CT, USA.The blue lines on the bottom panel (g) represent the longitudinal extent of mapped groundwater discharges.Gray lines in panel (f) represent the subset of groundwater discharge locations sampled for chemistry.Red boxes highlight examples of areas where solute concentrations vary dramatically over short distances.Violin plots of small stream groundwater discharge concentrations are shown on the right for comparison.

Figure 4 .
Figure 4. Particle tracking model flowpath lines along the mainstem of the Farmington River, CT overlain on 2019 National Land Cover Dataset.Gray to black lines represent predicted particle residence time from recharge to discharge.Blue boxes show modeled groundwater discharge rates along the river network.Black points represent not sampled streambank discharge locations, while colored points indicate sampled groundwater discharge NO3 − concentrations.

Table 2 .
Coefficients for best fit linear mixed-effects models for groundwater discharge nitrate (NO3 − ) concentrations in the Farmington River, CT, USA.Dependent variables are shown on the left with coefficients and significance ( * indicates p < 0.05) of each best fit model.Detailed model selection and regression output is included in SI (tableS1).Mulholland et al  2008)and the importance of in-stream N retention and removal in surface waters as shown by significantly less N flux in surface waters versus groundwaters(Wollheim et al 2006, Mullholand et al 2008).
Buss S R, Morgan P, Smith J W N and Bemment C D 2008 Nitrate attenuation in groundwater: a review of biogeochemical controlling processes Water Res.42 4215-32 Rosenberry D O, Briggs M A, Delin G and Hare D K 2016 Combined use of thermal methods and seepage meters to efficiently locate, quantify, and monitor focused groundwater discharge to a sand-bed stream Water Resour.Res.52 4486-503 Rosenberry D O, Engesgaard P and Hatch C 2021 Hydraulic conductivity can no longer be considered a fixed property when quantifying flow between groundwater and surface water Hydrol.Process.35 e14226 Sanford W E and Pope J P 2013 Quantifying groundwater's role in delaying improvements to Chesapeake Bay water quality Environ.Sci.Technol.47 13330-8 Scanlon B R, Reedy R C, Stonestrom D A, Prudic D E and Dennehy K F 2005 Impact of land use and land cover change on groundwater recharge and quality in the southwestern US Glob.Change Biol.11 1577-93 Shabaga J A and Hill A R 2010 Groundwater-fed surface flow path hydrodynamics and nitrate removal in three riparian zones in southern Ontario, Canada J. Hydrol.388 52-64 Stets E G, Sprague L A, Oelsner G P, Johnson H M, Murphy J C, Ryberg K and Riskin M L 2020 Landscape drivers of dynamic change in water quality of U.S. rivers Environ.Sci.Technol.54 4336-43 Sullivan C J, Vokoun J C, Helton A M and Briggs M A 2021 An ecohydrological typology for thermal refuges in streams and rivers Ecohydrology 14 e2295 Tesoriero A J, Duff J H, Saad D A, Spahr N E and Wolock D M 2013 Vulnerability of streams to legacy nitrate sources Environ.Sci.Technol.47 3623-9 Tesoriero A J, Liebscher H and Cox S E 2000 Mechanism and rate of denitrification in an agricultural watershed: electron and mass balance along groundwater flow paths Water Resour.Res.36 1545-59 Tesoriero A J, Stratton L E and Miller M P 2021 Influence of redox gradients on nitrate transport from the landscape to groundwater and streams Sci.Total Environ.800 150200 U. S. Census Bureau 1990 TIGER/line data (available at: www.census.gov/data/data-tools.html) U. S. Geological Survey 2023 National Water Information System (available at: http://waterdata.usgs.gov/nwis/)(Accessed 1 February 2023) Van Meter K J and Basu N B 2017 Time lags in watershed-scale nutrient transport: an exploration of dominant controls Environ.Res.Lett.12 084017 Van Meter K J, Van Cappellen P and Basu N B 2018 Legacy nitrogen may prevent achievement of water quality goals in the Gulf of Mexico Science 360 427-30 Vero S E, Basu N B, Van Meter K, Richards K G, Mellander P E, Healy M G and Fenton O 2018 Revue: L'état environnemental et les implications du décalage temporel du nitrate en Europe et Amérique du Nord Hydrogeol.J. 26 7-22 Ward B B 2008 Nitrification Encyclopedia of Ecology, Five-Volume Set (Academic Press) pp 2511-8 Wherry S A, Tesoriero A J and Terziotti S 2021 Factors affecting nitrate concentrations in stream base flow Environ.Sci.Technol.55 902-11 Winter T C, Harvey J W, Franke O L and Alley W M 1998 Ground Water and Surface Water-a Single Resourcex (U.S. Geological Survey) Circular 1139 Wollheim W M, Vörösmarty C J, Peterson B J, Seitzinger S P and Hopkinson C S 2006 Relationship between river size and nutrient removal Geophys.Res.Lett.33