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

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 and Pope 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 and Wagener 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₃ ⁻ 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.

Zoom In Zoom Out Reset image size

2.1. Study area

The Farmington River watershed (1572 km²) 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⁻¹ (Briggs et al 2021). The daily mean discharge is 30.9 m³ s⁻¹ (USGS NWIS, gage 01189995, 2010–2022).

Zoom In Zoom Out Reset image size

2.2. Field surveys of groundwater discharge

2.3. Sample collection

2.4. Particle tracking model

2.5. Riparian and contributing area land cover

2.6. Linear mixed-effects models of groundwater NO₃ ⁻

3.1. 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₃ ⁻ decreased from 7.42 to 0.37 mg-N/L within 30 m. TDN and NO₃ ^- were highly correlated across the watershed (rho = 0.937, p < 0.05, figure S1), and NO₃ ⁻ constituted over 90% of TDN, indicating that most of the N sampled was NO₃ ⁻. We generally observed low NH₄ ⁺ along the 25 km reach, yet there were rare locations with NH₄ ⁺ up to 2 mg l⁻¹, and NH₄ ⁺ was negatively correlated with NO₃ ⁻ (rho = −0.28, p < 0.05). Small stream concentrations and ranges of NO₃ ⁻ and NH₄ ⁺ in streams were lower than the mainstem (table 1), reflecting the shift from predominantly forest to mixed land use.

Zoom In Zoom Out Reset image size

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 (NO₃ ⁻), ammonia (NH₄ ⁺), total dissolved nitrogen (TDN), chloride (Cl⁻), sulfate (SO₄ ²⁻), dissolved oxygen (DO), dissolved organic carbon (DOC), dissolved nitrogen gas (N₂), nitrous oxide (N₂O).

		Farmington River mainstem	Small streams
NO3 − (mg-N/L)	Median:	1.1	0.49
	Range:	n.d.—11.65	n.d.—4.34
	Std. error:	0.25	0.09
	n:	107/107	123/123
NH4 + (mg-N/L)	Median:	0.05	0.03
	Range:	0.01–2.05	0.02–0.74
	Std. error:	0.03	0.01
	n:	107/107	123/123
N2O (ug-N/L)	Median:	3.83	3.78
	Range:	n.d.—146.41	n.d.—45.32
	Std. error:	2.02	0.59
	n:	97/107	123/123
N2 (mg l−1)	Median:	0.34	0.15
	Range:	n.d.—4.42	n.d.—1.17
	Std. error:	0.08	0.04
	n:	106/107	123/123
DO (mg l−1)	Median:	5.16	5.68
	Range:	0.1–10.55	0.49–10.56
	Std. error:	0.28	0.22
	n:	101/107	123/123
DOC (mg l−1)	Median:	0.84	0.73
	Range:	n.d.—3.63	0.31–2.40
	Std. error:	0.07	0.03
	n:	107/107	123/123
SO4 2− (mg-S/L)	Median:	3.55	2.89
	Range:	0.12–52.57	0.50–12.88
	Std. error:	0.76	0.20
	n:	107/107	123/123
Cl− (mg l−1)	Median:	27.5	19.63
	Range:	n.d.—2337.84	n.d.—123.31
	Std. error:	23.62	2.33
	n:	107/107	123/123

Denitrification end-products also varied considerably, indicating heterogeneity in N processing. Mainstem N₂O was positively correlated with NO₃ ⁻ (rho = 0.58, p < 0.05) whereas excess N₂, our surrogate measure of denitrification, was not (rho = −0.03, p > 0.05). Small stream N₂O and N₂ groundwater concentrations were significantly lower than mainstem concentrations and N₂O was positively correlated with NO₃ ⁻ (rho = 0.66, p < 0.05), whereas N₂ was not (rho = 0.02, p > 0.05). Only 60% of N₂ samples contained excess N₂ 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⁻¹ within 150 m (figure 3). Mainstem and small stream DO were not significantly different but were positively correlated with NO₃ ⁻ (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⁻¹) and were not correlated with NO₃ ⁻ 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₄ ²⁻ was positively correlated with NO₃ ⁻ 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⁻¹ along the mainstem and 0.03–123.31 mg l⁻¹ 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.

3.2. 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⁻³ to 18.8 m³ m⁻¹ (median = 3.8 m³ m⁻¹).

Zoom In Zoom Out Reset image size

3.3. Predicting spatial patterns of groundwater NO₃ ⁻

Our research shows that groundwater discharge N concentrations are highly heterogeneous both within and among streambank discharge features, and that NO₃ ⁻ 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₃ ⁻. 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 and Hill 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.

Throughout the watershed, groundwater well NO₃ ^-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₃ ⁻ data (table 1). Groundwater well NO₃ ⁻ 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₃ ⁻. For example, high groundwater NO₃ ⁻ 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₃ ⁻ concentrations and estimate groundwater flux to be 0.5 m d⁻¹ (Haynes et al 2023), groundwater NO₃ ⁻ load ranges from 6740 to 442,360 kg d⁻¹. Estimating surface water NO₃ ⁻ load from the mainstem using our 5th (0.78 mg-N/L) and 95th (1.92 mg-N/L) percentile surface water NO₃ ⁻ concentrations and mean daily river discharge (30.9 cms) results in a daily surface water load of 208–5,125 kg d⁻¹ of NO₃ ⁻. These estimates indicate that groundwater is the primary source of N loading within our study watershed, highlighting the challenge of empirically quantifying groundwater N contribution within entire watersheds (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). Groundwater discharge N also varied significantly between sampling years where NO₃ ⁻ 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₄ ²⁻ 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: N₂O and N₂) consistent with other groundwater studies (Tesoriero et al 2000, Hinshaw et al 2020, McAleer et al 2017, Chen et al 2014), however only N₂O was strongly related to NO₃ ⁻ concentrations (R² = 0.51, p < 0.05). Low NH₄ ⁺ coinciding with elevated DO and SO₄ ²⁻ suggest that nitrification (Ward 2008) and sulfur denitrification (Hayakawa et al 2021) could lead to elevated NO₃ ⁻ and N₂O in discharging groundwater. Most of our groundwater DO samples were above the 2 mg l⁻¹ anoxic threshold for denitrification (Rivett et al 2008). High groundwater DO concentrations support our finding that only 60% of samples contained excess N₂ 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₃ ⁻. 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 and Hare 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 sub-site) 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.

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₄ ²⁻ paired with low DOC concentrations indicated NO₃ ⁻ 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.

We acknowledge review comments from two anonymous reviewers and one USGS reviewer that substantially improved the manuscript. We would like to thank Fiona Liu for significant field and laboratory help as an undergraduate researcher. Any use of trade, firm, or product names is for descriptive purposes only and does not imply endorsement by the U.S. Government.

The data that support the findings of this study are openly available at the following https://doi.org/10.5066/P941XKST (Moore et al 2023).

No humans subjects were used in this research.

This research was funded by NSF-EAR Award Number 1824820 and USDA HATCH Project CONS00938.

The authors declare no conflict of interest