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 © 2023 The Author(s). Published by IOP Publishing Ltd 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.

Methods

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), andwater (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).

Linear mixed-effects models of groundwater NO 3

-

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

Results

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.31 mg 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.0 m). Euclidean distance to sampling locations ranged from 150.0-3400.8 m (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 ).

Predicting spatial patterns of groundwater NO 3

-

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.

Table 1. Groundwater discharge solute concentrations along the mainstem of the Farmington River, CT, and small streams. Bolding indicates significant differences of medians between mainstem and small stream groundwater chemical samples (Wilcoxon, p < 0.05). Counts (n) represent any chemical data missed during sample collection. Concentrations below detection are listed as n.d. (non-detect). Nitrate (NO3 -), ammonia (NH4 + ), total dissolved nitrogen (TDN), chloride (Cl -), sulfate (SO4 2-), dissolved oxygen (DO), dissolved organic carbon (DOC), dissolved nitrogen gas (N2), nitrous oxide (N2O).

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 varied significantly between sampling years where NO 3 -concentrations were lower during the high (2021) 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.