Introduction

Groundwater quality concerns in the karst-dominated Driftless Area of southeastern Minnesota, USA (Fig. 1) have led to extensive water quality monitoring over the past several decades. One of the most frequently detected contaminants is nitrate-nitrogen (reported as nitrite plus nitrate-nitrogen, NO2–N + NO3–N), hereafter simply referred to as nitrate. An assessment of domestic drinking water wells in areas with a dominance of agricultural land use in the Driftless Area by the Minnesota Department of Agriculture (MDA) from 2013 to 2019 found that 13.6% of the sampled wells (8,837) exceeded the United States Environmental Protection Agency (EPA) nitrate drinking water standard of 10 mg L−1. Additionally, 18.1% of km-scale stream sections in the Minnesota Driftless Area are listed or proposed to be listed as nitrate impaired, with half of these sites added since 2022 (Minnesota Pollution Control Agency, unpublished 2024).

Fig. 1
figure 1

The Driftless Area in relation to: a North America; and b the State of Minnesota. c Location of groundwater residence time study springs and wells. Map numbering of springs (1–16) and wells (17–30) corresponds to the numbering used in figures and tables. Deeply incised valleys are represented by dark green hues across the Prairie du Chien Plateau. Gray lines in inset map (a) represent state boundaries. Legend abbreviations: PdC Prairie du Chien, SL St. Lawrence, LR Lone Rock

A report completed in 2013 estimates that 89% of the nitrogen loading to surface waters in southeast Minnesota watersheds originates from cropland, with a substantial portion derived from groundwater (57%) and tile drainage pathways (23%). Most remaining nitrogen loads are from wastewater treatment plants, septic systems, urban runoff, and atmospheric deposition (MPCA 2013). Farmers and crop advisors are encouraged to use conservation practices and nitrogen best management practices (BMPs) to improve economic efficiencies while reducing nitrate leaching from crops. Some of these practices are regulated by statute through the passage of Minnesota’s Chapter 7020 Feedlot and Chapter 1573 Groundwater Protection Rules. However, measuring the effectiveness of these efforts beyond the field scale is challenging because practices may not produce immediate results where groundwater is sampled—for instance, monitoring results from some springs and streams do not appear to reflect relatively recent changes in local agricultural practices (Watkins et al. 2013; Runkel et al. 2014; Steenberg et al. 2014; Barry et al. 2023), nor event-scale climatic conditions (Luhmann et al. 2011). Groundwater nitrate concentrations are increasing at some sites while decreasing at others with similar agricultural practices. This phenomenon has been suggested to be caused in part by the complexity of the multiaquifer groundwater system in this landscape, which results in monitoring sites that have variable mixing of local and regional sources of water with distinctly different residence times and physical and chemical properties (Runkel et al. 2014; Barry et al. 2023). Residence time and legacy nitrate are considered key factors when evaluating the impacts of land-use practices on groundwater quality and trends (Puckett and Cowdry 2002; Katz 2004; Van Meter et al. 2018; Johnson and Stets 2020; Green et al. 2021; Ascott et al. 2021), but analysis and interpretations in the Driftless Area are limited.

This study conducted a trend analysis of nearly 35,000 nitrate water samples collected from 1,174 wells, springs, and stream locations. The residence time was quantified for a subset of springs and wells (Fig. 1) using a novel method based on the historical use of the herbicide alachlor and the concentration trends of its degradate. To date, this is the most in-depth and comprehensive analysis available, as it integrates multiple data sources and is interpreted within the context of a hydrogeologic framework. Results were combined with other independent age-dating methods and supporting data, including nitrogen fertilizer input histories, land use, and precipitation. These results inform nitrate trend interpretations and allow for an improved understanding of the relationships between land use and groundwater quality, which can be applied to complex multiaquifer systems with similar geology elsewhere.

Groundwater residence time

Groundwater ‘age’ or residence time is defined as the elapsed time along a flow path from the land surface to a discharge point such as a well, spring, or stream (Kazemi et al. 2006; Bethke and Johnson 2008; Green et al. 2021). For this study, the terms residence time, groundwater age, and lag time have similar meanings and are used interchangeably. The word ‘apparent’ often prefaces these terms, since residence time is based on interpretation from measured concentrations of environmental tracers in groundwater (Plummer and Freidman 1999) and reflects an integrated mixture of young and old water with multidimensional flow paths (Bethke and Johnson 2008; Cartwright et al. 2017; Alexander and Alexander 2018; Green et al. 2021). Estimates in this report using the alachlor ethane sulfonic acid (ESA) method are assumed to reflect the modern (post-1953) fraction of groundwater mixtures. Analysis of other tracers and methods, as described in the following, provides insight into premodern (pre-1953) age components. Where available, the ‘total’ mean age, defined here as the weighted average of both the modern and premodern components is provided.

Alachlor

Alachlor is a selective chloroacetamide (formerly acetanilide) herbicide that was first registered for use in the United States in 1969 and was used extensively on row crops in the 1970s and 1980s (EPA 1990a, b). Due to regulatory requirements established after registration, and more cost-effective alternatives, alachlor use in Minnesota declined rapidly in the 1980s, with little or no use by 2002. By 2016, the EPA canceled all registrations, making it illegal to sell or distribute alachlor (OLA 2020). Trends in concentrations of certain pesticides and their degradates are related to the timing of their introduction and use. Degradates, and to a lesser extent their parent compounds, are typically detected in groundwater following their introduction (Tesoriero et al. 2007). When detected in water, a minimum age is the difference between the year a chemical was first detected and the year it was first introduced (Plummer et al. 2003).

There are many complex chemical, physical, and biological processes that influence the fate and transport of pesticides in the environment but, in general, parent pesticides tend to transform relatively rapidly in aerobic soil and much more slowly in water (Barbash and Resek 1996; Gilliom et al. 2006). Therefore, pesticide degradates are often detected at a greater frequency in groundwater when compared to their parent compounds (Kolpin et al. 1995; Scribner et al. 2005; MDA 2023). One of the more frequently detected alachlor transformation products in water samples is alachlor ESA (Lee and Strahan 2003; MDA 2023). Alachlor ESA was a preferred age tracer for this study due to its high solubility in water, high leaching potential, and conservative chemical properties (EPA 1998a, b). Additionally, the usage of the herbicide is relatively well documented, which allows distinct periods of rise, peak, or decrease of usage to be statistically correlated with corresponding trends of alachlor ESA in groundwater. Hereafter, the study method in which alachlor ESA is utilized as a tracer is known as the alachlor ESA method.

Other tracers

All water age dating techniques have limitations; therefore, results are often referenced against multiple, independent methods (USGS 1999). In addition to the alachlor ESA method, this study uses conventional atmospheric tracers, including enriched tritium (3H), chlorofluorocarbons (CFCs), sulfur hexafluoride (SF6), and tritium–helium (3H/3He) (Phillips and Castro 2003; Plummer et al. 2006; Faulkner et al. 2022). To provide insights into the extent of the premodern fraction of groundwater mixtures, carbon-14 (14C) dates for Minnesota waters were used (Alexander and Alexander 2018).

Study area, climate, and regional hydrogeologic context

The springs and wells evaluated in this study are in the Driftless Area ecoregion (Omernik 1987; White 2020) of southeast Minnesota. The Driftless Area also includes parts of Wisconsin, Iowa, and Illinois, covering an area of 62,000 km2 (24,000 mi2) in the upper Midwest (Fig. 1a,b). The Driftless Area is unique because it is dominated by a relatively high relief landscape of incised sedimentary bedrock, without the thick and extensive cover of mostly clay-rich diamicton (till) deposited during the most recent glaciations elsewhere in the region. This landscape transitions westward to feature progressively more extensive till cover (Fig. 1b), referred to as the Western Corn Belt Plains (WCBP; Omernik 1987). A subset of streams evaluated in this study are in this till-dominated landscape which, unlike the Driftless Area, is extensively tile drained.

The climate of southeast Minnesota consists of cold winters and hot, humid summers. The average annual temperature between 1991 and 2020 was 7.2 °C (45.0 °F), and precipitation was 908.7 mm (35.8 in; DNR 2023b). The geology of the Driftless Area is characterized by subhorizontal Lower Paleozoic sedimentary bedrock formations ranging from 1 to 60 m (3–200 ft) thickness and are dominated by either carbonate rock (limestone and dolostone); fine to coarse-grained quartz-rich sandstone (coarse clastic); or very fine-grained sandstone, siltstone, and shale (fine clastic, Fig. 2). Due to the thin cover of mostly sandy to silty unconsolidated sediment (generally <15 m, 50 ft), the landscape largely reflects the topography of the bedrock surface. At a regional scale, two particularly thick bedrock layers, composed mostly of karstic carbonate rock, form regionally extensive plateaus (Fig. 1c; Runkel et al. 2014). The western Upper Carbonate Plateau (UC Plateau) is composed of Devonian and Middle Ordovician-aged carbonate and shale-dominated formations. The eastern Prairie du Chien Plateau (PdC Plateau) is dominated by the Lower Ordovician dolostone of the Prairie du Chien Group and extends to the edge of the Mississippi River. Incised streams are especially pronounced across the PdC Plateau, with valleys as deep as 150 m (490 ft), penetrating deeply into Cambrian siliciclastic strata, and truncating multiple aquitards. These areas are referred to as deeply incised valleys in this study (Fig. 1c).

Fig. 2
figure 2

Hydrostratigraphic positions of springs (1–16) and open hole intervals of wells (17–30) evaluated in this study, modified from Runkel et al. (2014). Springs emanate from three general stratigraphic positions A–C as also shown in Fig. 3: a Galena-Cummingsville Formation (Galena), b Oneota Dolomite and Jordan Sandstone (PdC-Jordan), and c St. Lawrence and Lone Rock Formations (SL-LR)

Geological layers designated as aquifers easily transmit water through conduits, fracture networks, and/or porous media. Sandstone aquifers have a high matrix porosity (15–30%) and vertical hydraulic conductivity (Kv) values of up to 100 m day–1 (Runkel et al. 2003, 2006). Carbonate aquifers transmit water largely via secondary porosity features (conduits and fractures) with typical matrix porosity values of less than 5% and matrix Kv values of less than 10–5 m day–1. Layers designated as aquitards restrict flow and increase the residence times of groundwater in underlying aquifers. Aquitards are dominated by very fine-grained sandstone, siltstone and shale, or carbonate rock with poorly connected vertical fractures (Fig. 2). The aquitards are variably leaky, as exemplified by the three key aquitards in this study, with the Cummingsville-Glenwood and St. Lawrence aquitards having Kv values as low as 10–7 m day–1 and the Oneota aquitard featuring much higher Kv values ranging from 10–3 to 10–5 m day–1 (Runkel et al. 2003; Tipping et al. 2006; Runkel et al. 2018).

The interlayered aquifers and aquitards result in an anisotropic groundwater system, limiting vertical flow and promoting lateral flow that discharges as baseflow to springs and streams (Figs. 2 and 3). This results in pronounced stratification in groundwater age and anthropogenic contaminants, as documented in numerous regional groundwater reports and maps (Zhang and Kanivetsky 1996; Berg 2003; Petersen 2005; Barry 2021; Barry et al. 2018, 2023) and summarized in Runkel et al. (2014, 2018) and Steenberg et al. (2014). Although aquitards in this setting exhibit very low matrix permeability, they commonly contain high permeability, bedding-parallel partings (Fig. 2) in a horizontal direction which can yield large quantities of water to wells (Runkel et al. 2018). Many springs in Minnesota, including several in this study, discharge from such partings in aquitards (Green et al. 2012; Barry et al. 2015; Steenberg and Runkel 2018; Runkel et al. 2018). Additionally, aquitard integrity is diminished at deeply incised valley edges, where enhanced vertical fracturing results in the mixing of upgradient stratified waters discharged to springs and baseflow to streams (Runkel et al. 2018; Barry et al. 2018).

Fig. 3
figure 3

Schematic west to east cross section across the bedrock-dominated landscape of southeastern Minnesota showing the position of spring groups A, B, and C and study wells within the context of the hydrogeologic system. The stratigraphic position of the individual springs within each spring group is shown in Fig. 2. Dark green shading depicts areas with median nitrate concentrations typically above 4 mg L−1 and lighter green shading less than 1 mg L−1. The inset map shows well locations with nitrate concentrations between 4–10 mg/L and >10 mg L−1. Data sourced from the ‘well full nitrate’ dataset (dataset 1 of the electronic supplementary material, ESM 1)

Hydrogeologic context of residence time study springs

Spring discharge generally represents an integrated mixture of flow paths and contains variable proportions of younger and older waters (Runkel et al. 2014; Barry et al. 2023). Springs with high landscape positions in the unconfined karstic carbonate rocks of the UC Plateau commonly feature groundwater discharges that are strongly dominated by relatively recent recharge components near the top of the phreatic zone or top of the water table (Fig. 3). In contrast, incised valleys of the PdC Plateau commonly breach multiple bedrock aquifers and aquitards, which varies the proportion of relatively young water from the unconfined water-table aquifer and older confined groundwater discharged from springs in this setting (Runkel et al. 2014; Barry et al. 2023). The 16 springs with residence times estimated in this study (Figs. 1, 2, 3; Table 1) are classified into three groups (A, B, C) based on their hydrostratigraphic (Steenberg and Runkel 2018) and landscape positions. Group A springs emanate from the upper part of the Galena Group in the UC Plateau and are referred to as Galena springs. Group B springs, with lower landscape positions within the Oneota Dolomite or upper Jordan Sandstone of the PdC Plateau, are grouped and abbreviated as PdC-Jordan. Group C includes springs in the lowermost landscape positions that emanate from the St. Lawrence and upper Lone Rock Formations, referred to as SL-LR springs.

Table 1 Site details for study springs and wells. Sites were grouped by type and sorted from highest to lowest hydrostratigraphic position. Estimated mean apparent residence time results using the alachlor ESA method are displayed in the age modern and decadal category columns. Atmospheric tracer results reference data from Faulkner et al. (2022)

The five Group A Galena study springs are located within shallow incised valleys of the UC Plateau. One emanates from the Stewartville Formation; the other four emanate from lower in the stratigraphic section in the upper part of the combined lower Cummingsville-Glenwood Aquitard (Steenberg and Runkel 2018; Fig. 2). The presence of this thick underlying aquitard, composed largely of shale, and a downward gradient across the aquitard (Zhang and Kanivetsky 1996) indicate that Galena springs have no upwelling from deeper confined groundwater. Spring discharge is thus dominated by water from a karstic, unconfined, water-table aquifer that is commonly perched near valley edges (Carter et al. 2011). Dye tracer studies demonstrate conduit flow from sinkholes and stream sinks to many Galena springs (Alexander and Milske 1986; Ustipak et al. 2012; Kuehner et al. 2017; Barry et al. 2020; DNR 2021b). Galena springs are also characterized by abrupt and significant changes in temperature, flow, and nitrate concentrations in response to recharge (Luhmann et al. 2011; Barry et al. 2020); however, these springs maintain constant baseflow through the slow release of water in storage.

The Group B PdC Amherst Spring emanates from the upper Oneota Dolomite, a position hydrogeologically considered the lowermost part of the Shakopee Aquifer (Tipping et al. 2006). It is along or near a well-recognized karstic high-transmissivity zone (Runkel et al. 2003; Tipping et al. 2006; Alexander et al. 2013). Although positioned more deeply in an incised valley than Galena springs, and lower in the stratigraphic section, it shares several hydrologic characteristics with them. Dye tracing substantiated a rapid karst connection (DNR 2021b), and temperature monitoring showed abrupt and significant changes in response to recharge (Luhmann et al. 2011; unpublished MDA data). Furthermore, unlike the other study springs that emanate from the Oneota Dolomite (described further on), the elevation of this spring is above the potentiometric surface of the underlying confined Jordan Aquifer (MDH 2011a).

Group B PdC-Jordan springs are positioned lower in the stratigraphic section in hydrogeologically complex settings that include the potential for significant contributions of older confined regional flow. Five of the PdC-Jordan springs emanate from the Oneota Aquitard at a surface elevation lower than the hydraulic head of the underlying Jordan Aquifer, allowing for the possibility of significant upwelling. Upwelling was demonstrated at the Lanesboro State Fish Hatchery (SFH), where dye tracing from a (684602) in the Jordan aquifer demonstrated an upward gradient towards the Main Spring (Tipping et al. 2006; Barry et al. 2023). Temperature profiles of the hatchery springs (Luhmann et al. 2011; unpublished MDA data) lack event-scale (hours to days) perturbations, suggesting an isothermal contribution from a confined aquifer that is buffered from rapid recharge from the land surface.

Despite their stratigraphically lower position within the Oneota Aquitard, a component of relatively recent water from the overlying water-table aquifer is expected to contribute to PdC-Jordan springs. Dye tracing from sinkholes and stream sinks to the Lanesboro SFH springs documented a component of locally derived, recently recharged water, traveling via fast flow paths (Wheeler 1993; Barry et al. 2023). Mixing models suggest approximately half of the discharge emanating from the Lanesboro Main Spring comes from the water-table aquifer and the other half from the confined Jordan Aquifer (Barry et al. 2023). The proportions of older confined water to younger water at the other Oneota Aquitard springs have not been estimated.

Two PdC-Jordan springs, Big Spring and Hwy 76, are positioned hydrostratigraphically deeper, at or near the top of the Jordan Aquifer, and have temperature profiles similar to the four Oneota Aquitard springs (Luhmann et al. 2011; unpublished MDA data). Buffered temperature profiles coupled with their positions fully beneath the Oneota Aquitard suggest they may have a greater proportion of confined Jordan Aquifer water compared to the Oneota Aquitard springs.

The four Group C SL-LR springs are the hydrostratigraphically lowermost study springs, emerging from St. Lawrence and Lone Rock Formations beneath parts of two aquitards: the Oneota, and the upper part of the lower Jordan-St. Lawrence. Dye tracing and water chemistry show these springs can have minor contributions of recent water from nearby losing reaches of small tributaries (Green et al. 2008, 2012; Barry et al. 2015, 2018), but their hydrostratigraphic position and hydraulic head conditions dictate that flux includes confined water from upgradient of the springs and upwelling from underlying confined aquifers (Runkel et al. 2014; Barry et al. 2015, 2018; Barry 2021). Temperature profiles at each of the SL-LR springs, other than Jerry Lee Spring (JLS; not monitored for temperature), lacked event scale perturbations, which is indicative of buffering from rapid recharge (Luhmann et al. 2011; unpublished MDA data).

Hydrogeologic context of residence time study wells

Residence time was also quantified for a subset of 14 wells (Table 1). Two wells on the UC Plateau draw water from unconfined karstic aquifers. The others, from the PdC Plateau or outer edge of the UC Plateau, draw water from parts of the karstic PdC Group, the Jordan, and upper St. Lawrence aquifers. Some of these wells are completed in confined aquifers (Figs. 2 and 3; Table 1). Well construction records indicate that 10 wells have grouted casings; the remaining four have inadequate records to determine grouting status. These four wells are Galena and PdC wells with no overlying aquitards.

Methods

Time series alachlor ESA concentrations, combined with historical alachlor usage, were used to approximate residence time from 13 springs and four well sites across varied hydrostratigraphic settings. Atmospheric tracers and 14C from previous studies were also used. These residence times, combined with other supporting information, such as historical nitrogen fertilizer use, precipitation, and cropping trends, were used to interpret nitrate monitoring results from some of these same springs and wells. This information was then used to help interpret a comprehensive nitrate trend analysis of 1,174 wells, springs, and stream sites. Trends and concentrations were analyzed for 1,097 wells, 16 springs, and 22 stream sites located in northeast Iowa (RC&D 2023) and 21 stream stations located in southeast Minnesota. An additional 18 streams in the tile-drained glacial till landscape in the WCBP were also evaluated as a comparison to the Driftless Area to provide insights about the influence of precipitation patterns and land use on nitrate trends in short residence time groundwater systems. Nitrate data sources were compiled from a variety of local, state, and federal agencies. Additional details about data sources and sampling methods can be found in the electronic supplementary material (ESM 1). Additional details on the methods, data sources, and analyses are presented as tables and figures in ESM 2 and ESM 3.

Key assumptions

The alachlor ESA method uses four primary assumptions to estimate residence time:

  1. 1.

    Alachlor use in the recharge area of the monitored spring or well was consistent with state or national use patterns. The resolution of alachlor use information is insufficient to confirm this assumption, but it is deemed reasonable given that herbicide choices were limited in the 1970s and 1980s (EPA 1990a, b) and detection of alachlor degradates in Minnesota groundwater throughout agricultural areas is common (MDA 2023).

  2. 2.

    Alachlor ESA moves through soil and geologic layers at the same rate as water. Therefore, residence time quantification assumed a linear dependency between statewide alachlor sales, alachlor herbicide applications on the land surface, and alachlor ESA concentrations measured in groundwater samples.

  3. 3.

    The ESA degradate of alachlor is persistent and conservative in groundwater, and attenuation is primarily the result of dispersion rather than biotic or abiotic degradation or chemical sorption processes. There were no detections of the parent alachlor compound in the groundwater samples of this study, consistent with the findings of Tesoriero et al. (2007).

  4. 4.

    Alachlor was first introduced in 1969, so estimated ages are limited to the modern fraction of groundwater mixtures.

Data comparisons and use of other age-tracer datasets

Residence time estimates from the alachlor ESA method were compared to independent age tracers:

  1. 1.

    Groundwater age was estimated using atmospheric tracers in samples from 10 wells and four springs (Faulkner et al. 2022). The tracers included CFCs (CFC-11, CFC-12, CFC-113), SF6, and 3H/3He, along with ancillary major dissolved-gas data (N2, O2, Ar, CH4 and CO2). Additionally, a lumped parameter model called TracerLPM developed by Jurgens et al. (2012) was used to model mean residence time and groundwater age distributions.

  2. 2.

    Presence or absence of degradates associated with herbicides atrazine, metolachlor, and acetochlor was used as an indicator for residence times. For instance, frequent detections of desethylatrazine, metolachlor ESA, and acetochlor ESA were considered an approximate marker for post-1958, 1976, and 1994 groundwater residence times, respectively, with those dates corresponding to when those pesticides were introduced in Minnesota (EPA 1990a, b, 1995, 1998a, b; USDA 1994).

  3. 3.

    Tritium data collected by the DNR Groundwater Atlas Program was used to identify the extent of modern water (DNR 2023a). Tritium data were available for the following springs: Crystal, Hwy 76, Big Spring, and Lanesboro. To inform the extent of the premodern fraction of groundwater mixtures, a larger dataset of 964 wells was compiled using tritium and 14C results (Alexander and Alexander 2018).

Spring, well, and watershed supporting information

Table 1 provides supporting information for study springs and wells including plateau setting, primary groundwater source, and relative depth. Additional hydrostratigraphic position details for each spring can be found in Table S1 of ESM 2. Recharge areas (springsheds) for each spring were estimated using a combination of dye tracer results, potentiometric maps, and surface-water drainage areas. For Galena springs, row cropping acres were selected within springsheds defined through dye tracing, available through the Minnesota Groundwater Tracing Database (MGTD; DNR 2021a). Springshed areas shown in the MGTD are, in some cases, underestimates of the maximum area of the springshed, as mapping takes many dye traces across different water-table conditions to fully delineate a springshed. For springs without dye trace delineated springshed boundaries, all non-Galena springs, the apparent recharge area was defined using a combination of potentiometric contours (MDH 2011a) and surface-water drainage areas from the DNR watershed suite (DNR 2009). Study springs range from shallow, unconfined water-table springs to springs emanating from deeply incised valleys. The relative depths of springs, used for comparisons and correlations, were calculated by taking the difference between the median elevation within the springshed and the elevation of each spring. Elevations were derived from a one-meter digital elevation model. Well depth represents the distance from the land surface to the completed depth of the well. Spring and well depth provide a relatively simplistic metric of vertical travel distance from recharge to discharge, although depth is not the only factor controlling residence time. Hydraulic gradient, in particular artesian conditions at groups B and C springs, can greatly influence residence time estimates and was included in interpretations. Land-use trends from 2010 to 2022 were evaluated for spring and well recharge areas and watersheds using 30-m pixel resolution rasters from CropScape (USDA 2022).

Statistical methods

Correlation and monotonic trend analysis of nitrate and alachlor ESA concentration time series were conducted using Mann–Kendall’s correlation test (Ktau). When datasets contained one or more samples below the analytical detection limit, the Akritas-Theil-Sen (ATS) slope estimator was used; hereafter, both methods are simply abbreviated as Ktau. This nonparametric method and the resulting Theil-Sen line, represent the linear median and are less sensitive to outliers and skewed residuals (Helsel et al. 2020). Trends were evaluated for sites with a minimum of 8 monitoring years for springs and streams, and a minimum of 6 years for wells, while the average number of site years analyzed was 15, 18, and 11 years for springs, streams, and wells, respectively. The primary trend evaluation period was from 1997 to 2021, as this was the period when the data available for most monitoring sites met the criteria. Statistical significance was determined when p-values were less than 0.05, and the lower and upper 95% confidence interval for the slope prediction did not contain zero. The Mann–Whitney test was used to evaluate differences between groups. Descriptive statistics containing left-censored values with multiple detection limits were computed using robust regression on order statistics (ROS; Helsel and Cohn 1988; Helsel 2005). Analysis was conducted using open-source coding developed in R (R Core Team 2021) and commercial software (Minitab™ ver. 21) using methodologies and scripts developed by Helsel et al. (2020). R packages utilized included applied environmental statistics (AES), data analysis for censored environmental data (NADA and NADA2), and survival and environmental statistics (EnvStats) packages.

Quantifying residence time

Nonparametric linear regression models were used to generate equations for a subset of years in which statewide alachlor usage demonstrated a significant increasing or decreasing trend. Similarly, linear models were applied to corresponding alachlor ESA concentration trends measured in groundwater from springs (Fig. 4a,b) and wells. The difference in years of the alachlor usage and ESA concentration regression models was used to estimate the residence time, similar to an approach by Crawford and Lee (2015).

Fig. 4
figure 4

Examples of linear time-series alachlor ESA concentrations measured from springs: a decreasing concentrations at Spring Creek Spring (PdC-Jordan); b increasing concentrations at Hwy 76 Spring (PdC-Jordan). Examples of non-linear: c decreasing concentrations at Burr Oak Spring (Galena); d near or past peak concentrations at Big Spring (PdC-Jordan) and Peterson Spring (SL-LR)

To generate equations of the fitted lines, the following regression models and methods were used: (1) Ktau with a Microsoft® Excel® add-in Solver; (2) nonlinear regression models and intersection with the alachlor ESA detection limit; and (3) peak-to-peak method. When using the Ktau method, if there was more than one detection limit during the study period, the highest value was subtracted from the intercept to create an adjusted y-intercept.

Method 1, or the Ktau method, the x-intercept was calculated by dividing the adjusted y-intercept by the slope. The difference between the predicted x-intercept of the alachlor ESA concentration year and the normalized alachlor use year was considered the estimated residence time in years. Alachlor use was normalized by dividing the mean alachlor ESA concentrations by the mean alachlor use during the trend period. In many cases, the slopes of the fitted lines were not parallel, which can result in exaggerated differences between the two intercepts. To estimate a more centralized difference between the fitted lines, the Excel® add-in Solver was used to identify the optimal combination of values that minimized the residual sum of squares (RSS) or error of the model. This approximates the difference between the center of the fitted lines rather than one specific point at which they cross the x-axis. The year where the RSS was either at a minimum or least error, represented the best approximation for the difference between the fitted lines.

Method 2, or the nonlinear intersect method, was used for instances when alachlor ESA concentrations were asymptotic near the detection limit. The point at which the lower 95% confidence interval of the nonlinear model intersected with the lab detection limit was the estimated x-intercept year. Nonlinear models that were evaluated included quadratic, exponential, and asymptotic convex, with the latter providing the best-fit results. This nonlinear solving method was optimal for all Galena springs because alachlor ESA concentrations were consistently decreasing and plateauing near the detection limit and exhibiting long-concentration tails, as demonstrated by Burr Oak Spring in Fig. 4c.

Method 3, or the peak-to-peak approach, was used for sites like Big Spring and Peterson Spring (Fig. 4d), which emanate from the top of the Jordan and Lone Rock aquifers, respectively. The year in which alachlor ESA concentrations had presumably peaked was subtracted from the year alachlor use had peaked, which was quantified by solving for the second derivative of the quadratic equation (additional method details are provided in ESM 3).

Results

Residence time

Historical alachlor use

Estimated alachlor use in Minnesota (Fig. 5) was created using a compilation of data sources detailed in ESM 1. Alachlor use in Minnesota climbed substantially after it was first introduced in 1969 and peaked in 1980 at just over 5,000 metric tons (mt). Over 20 years, from 1980 to 2000, alachlor use declined by 99% at an average of ~5% year–1 (260 mt year–1), whereas after 2016, alachlor was no longer registered in Minnesota. Since usage was effectively reduced to only 1% of the peak 1980 usage, 2002 is considered the effective end year. This pattern is consistent with national and regional alachlor use estimates, by Scribner et al. (2005) and Frick and Dalton (2007), except their analyses show peak usage in 1975, which is 5 years earlier than the 1980 estimate of this study. The discrepancy is likely due to this analysis including 1980 usage information from Hanthorn et al. (1982a, b) that was not included in previous studies.

Fig. 5
figure 5

Estimated alachlor use and average rate of change in Minnesota from 1971 through 2015. Alachlor was first registered for use in 1969, with peak usage occurring in 1980 and its use effectively ended in 2002

Spring and well alachlor ESA concentrations and trends

Alachlor ESA was not sampled at three of the study springs, whereas of the remaining 13 study springs, alachlor ESA was detected in 79 to 100% of the samples. Average alachlor ESA concentrations at each spring ranged from 102 to 577 ng L−1. Trend analysis results are shown in Fig. 6a. Only Big Spring lacked a statistically significant trend (p-value =0.066), which is attributed to alachlor ESA concentrations at or near the peak. Peterson Spring was interpreted to be near peak concentration, while Hwy 76 spring is also likely near peak, but additional sampling years are needed. There is a clear relationship between alachlor ESA trends, hydrostratigraphy, and spring depth (Fig. 6b). All Galena and most PdC-Jordan springs have a decreasing alachlor ESA trend (p < 0.001). In contrast, all SL-LR springs, and one deep PdC-Jordan spring, Hwy 76, show significantly increasing alachlor ESA trends. Average concentrations at four study wells ranged from 151 to 926 ng L−1. Consistent with springs, alachlor ESA concentrations were decreasing in relatively shallow Galena and PdC wells. One well (15,747), with a likely groundwater contribution from the Jordan Aquifer, had an increasing trend.

Fig. 6
figure 6

Concentrations, depth, trends, and residence time estimates using the alachlor ESA tracer method for 13 study springs. a alachlor ESA in ng L−1, (n = 677). The height of the upper and lower box represents the interquartile range (middle 50%); the black, dotted line represents the median concentration. Sites with statistically significant trends at the 0.05 level are labeled percentages (orange = increasing, blue = decreasing). Percentages represent the average annual rate of change. b Spring depth: depth is the difference between median springshed elevation and spring elevation; c estimated residence time of study springs using a relationship between Minnesota alachlor use and ESA concentration trends in springs. Lanesboro SFH underpredicted (asterisks); adjusted using methods described further on; d nitrate–N in mg L−1 (n = 691) interquartile concentration range and trends. Charts a and d reflect sampling years 2003 to 2022

Residence time estimates

Residence time estimates of the modern component of groundwater for the 13 study springs using the alachlor ESA method ranged from 8 to 41 years (Fig. 6c; Table 1). Residence times using atmospheric methods, which include both modern and premodern components, are also listed in Table 1 in the ‘Age total’ column. Including the premodern component for the study springs increased ages by 48%, resulting in a range of 32 to 78 years. The footnotes in Table 1 indicate agreement or disagreement between independent tracer methods. Except for Lanesboro SFH Upper Spring, there was good agreement between tracer methods. At Hwy 76 Spring, paired comparisons of modern age estimates showed remarkable agreement with quantified ages of 38 and 41 years, respectively, (Table 1; Fig. 7a) for the atmospheric and alachlor ESA methods.

Fig. 7
figure 7

Relationship between mean apparent groundwater residence time and depth using independent age tracer methods within the Driftless Area of southeast Minnesota. a springs and wells by tracer method (n = 31); b separated by springs (n = 17) and wells (n = 14)

Anomalous ages were observed at four sites. Poor agreement at the Lanesboro SFH Upper Spring is based on acetochlor and desethylatrazine concentrations, suggesting an age substantially older than the one decade estimated using the alachlor ESA method. The one-decade estimate is also inconsistent with the predicted age based on hydrostratigraphic position and depth. While the cause of this anomaly is unclear, one possibility is that alachlor usage may have ended much earlier within the springshed or was not consistent with statewide usage. Frequent detections of desethylatrazine since 2003 suggest a residence time of less than 45 years (marker 1958 subtracted from year 2003), while an absence of acetochlor ESA detections through 2022 implies an age greater than 25 years. Therefore, a mean residence time minimum of 25 years for Lanesboro SFH is considered more plausible than the alachlor ESA estimate. Age estimates for Crystal Creek Spring and wells 219,159 and 227,011 using atmospheric tracers were also anomalous. Although the young age of 8 years for the modern fraction for Crystal Creek Spring is generally consistent with other Galena spring water dated using the alachlor ESA method, the total age of 40 years is much older. Only two atmospheric tracers were used at Crystal Creek Spring, which may have inflated the premodern component of the estimate. Additionally, there are challenges to using atmospheric tracers in karstic springs that may increase age estimates (Pérotin et al. 2021). The reason for the anomalously low modern age of 17 years for well 219,159 and 14 years for well 227,011 is unclear; incomplete grouting of these karstic Prairie du Chien wells or improper grouting or sealing of nearby wells may be factors.

Residence time relationships with hydrostratigraphic position and depth

Spring residence time is related to hydrostratigraphy and depth (Fig. 6b). Galena springs on the UC Plateau have the shortest residence time, of about one decade. Residence times for PdC and PdC-Jordan springs range from two to three decades, with longer times corresponding to springs emanating from the deepest stratigraphy (Big Spring and Hwy 76). These two springs also have high alachlor ESA concentrations (Fig. 6a), related to peak alachlor usage (Fig. 5). Residence times for the deepest hydrostratigraphic positions, the SL-LR springs, are three to four decades.

Residence times for wells show a similar relationship. Shallow wells of the PdC Plateau without an overlying aquitard (wells 19–21, Table 1) have residence times of two to three decades. Wells penetrating deeper into the Oneota aquitard and/or the underlying Jordan and St. Lawrence Formations (wells 22–30, excluding wells 24 and 28, Table 1) have residence times of three to six decades, with the longest residence times corresponding to wells drawing exclusively from the Jordan and St. Lawrence aquifers.

Statistical comparison of residence time to depth for study springs and wells (Fig. 7a, Table S10 of ESM 2) shows a significant correlation (p < 0.001, Ktau 0.45). The regression model indicates that for every 30.5 m (100 ft), residence time increases by 6 years (3–9 years, 95% confidence interval). Both the alachlor ESA and atmospheric tracer methods plotted consistently with this regression model. Paired comparisons between the modern mixture component at Hwy 76 Spring indicated there was only a 3-year difference (Fig. 7a). Figure 7b separates springs and wells. The correlation was significant for springs (p < 0.001, Ktau 0.67), but was not significant for wells (p = 0.162, Ktau 0.28). Wells 219,159 and 227,011 had anomalously low ages for their given depths. If these wells are excluded, the correlation becomes significant (p = 0.003, Ktau 0.65). At shallow depths between 30–60 m, spring residence times were typically 10–20 years shorter than wells at the same depth. These differences decreased with depth, with fewer differences below 100 m.

Premodern to millennia-aged well water

For several decades, Minnesota researchers have sampled wells for tritium to aid in residence time interpretations. Methods developed by USGS and later modified for application in Minnesota to account for a decline in atmospheric tritium deposition were applied to samples in this study (Michel et al. 2018; Lindsey et al. 2019; DNR and MDH 2020). Wells are classified as modern, mixed, or premodern (pre-1953 and mostly pre-1953) using this methodology.

Of the 233 well samples classified as premodern tritium in the study area, 97% had nondetectable to ≤1 mg L−1 nitrate. Well depth for the tritium-dated wells varied by geologic setting but averaged 131 m. A subset of 23 of the premodern wells with sufficient data for trend analysis showed about half were below detection and half had no trend. One well classified as mostly premodern had an increasing trend.

A subset of the premodern wells (n = 122) was dated using 14C methods (Alexander and Alexander 2018). The average age of 14C dated water is ~6,500 years, with 87% of wells completed below regional aquitards. Well depth for the 14C dated wells was variable by geologic setting, but averaged 132 m, with 95% of wells featuring nitrate concentrations ranging from below the limit of detection to <1 mg L−1 nitrate–N. Although these deep wells do not provide a direct comparison to the 14 relatively shallower wells with residence time determined in this study, the tritium and 14C datasets reveal the presence of much older and unimpacted groundwater beneath the regional aquitards of the study area.

Nitrate trend analysis

Spring nitrate and trends

Nitrate was detected in all samples from the 16 spring monitoring locations where the average nitrate concentrations ranged from 4.1 to 13.9 mg L−1. Galena springs and the Amherst PdC spring had the highest concentrations, above or just below the EPA 10 mg L−1 drinking water standard (Fig. 6d). Concentrations were notably lower for PdC-Jordan and SL-LR springs, with the latter averaging less than 5 mg L−1. There was a strong negative correlation between spring nitrate concentrations and depth (Fig. 6b,d), where concentrations decreased by 2.9 mg L−1 for every 30.5 m (100 ft) (n = 16, p < 0.001, Ktau −0.78).

Trend results are shown spatially in Fig. 8a and relative to hydrostratigraphy in Figs. 6d and 9b. Of the 16 sites, 25% (4) had no trend, while 63% (10) were increasing, and 13% (2) were decreasing. The four sites with no trends were all Group A Galena springs on the UC Plateau. The two springs with decreasing trends were Crystal Creek and Amherst. All remaining B and C springs had increasing trends and were in deeply incised valleys on the PdC Plateau.

Fig. 8
figure 8

Nitrate trends from time-series concentrations within the Driftless Area of southeast Minnesota: a springs (n = 16); b streams (n = 43); c wells: no trend (n = 541) or below detection (n = 374); d wells: increasing (n = 140) and decreasing (n = 42). Additional stream sites within the Upper Iowa Watershed in Iowa were also included (b). Trends were analyzed primarily between years 1997 to 2019 for wells and 2000 to 2021 for streams and springs

Fig. 9
figure 9

a box plots of spring, stream, and well nitrate concentrations evaluated for trends in the Driftless Area. Yellow boxes represent groundwater pathways that are primarily above regional aquitards, blue primarily below, and no color as mixed. The height of the upper and lower plots represents the interquartile range (middle 50%); the black horizontal line represents the median and the whiskers represent the bottom and top 25%. b nitrate trends as a percent of sites that are increasing (Inc.), decreasing (Dec.), no trend (NT), or below the lab detection limit (ND). Springs, streams, and wells are grouped by geologic setting and hydrostratigraphy. Data filtered to years 2000 to 2021 to allow equal comparison. Note that stream data represent a subset of watersheds (27 of 43) from Fig. 8b that were spatially contained within specific geologic settings. Similarly, well data represent a subset of 300 wells from Fig. 8c,d that were filtered for specific hydrostratigraphic and well construction properties

Stream nitrate trends

Nitrate was detected in all samples from the 43 stream sites in the Driftless Area. Average concentrations by watershed ranged from 1.9 to 12.2 mg L−1 (Fig. 8b displays trend results spatially). During the period 2000 to 2021, 47% (20) of the 43 stream sites showed no trend, 51% (22) were increasing, and 2% (1) were decreasing. More than 90% of the increasing stream sites were located on the PdC Plateau (Fig. 9, ID e), a result similar to springs.

For comparison, concentrations from 18 tile-drained watersheds in the glacial till landscape of the WCBP ranged from 3.2 to 14.7 mg L−1. In the WCBP, 56% (10) of the sites showed no trend, while 6% (1) were increasing, and 39% (7) were decreasing. Analysis of stream trends over different climatic periods showed there were more ‘no trend’ or ‘decreasing trends’ in the glacial till landscape during wet years compared to streams in the Driftless Area (Fig. S1 of ESM 2). With only one exception, all increasing sites were confined to the Driftless Area, regardless of different trend evaluation start and end periods. When accounting for variability caused by flow on a subset of stream sites between years 2012 and 2021, the only increasing trend site was on the PdC Plateau. The remaining seven sites with decreasing flow-weighted concentrations were in till-dominated watersheds.

Well nitrate and trends

Nitrate trend analysis of 1,097 individual wells revealed that 49% (541) had no trend, 34% (374) had concentrations below the lab detection level, 13% (140) were increasing, and 4% (42) were decreasing (Fig. 8 c and d). Among wells with increasing trends, 21% were in the UC Plateau while 79% were located within the PdC Plateau. The rate of nitrate change expressed as a median of the trend slope was +0.09 and –0.27 mg L−1 year−1, respectively, for increasing and decreasing trends across all spring, stream, and well sites. A subset of 300 wells filtered by hydrostratigraphic and well construction details (Fig. 9) revealed median nitrate concentrations are significantly lower for wells cased and grouted beneath regional aquitards. The greatest percentage of increasing trend wells, which ranged from 9 to 13% (Fig. 9, IDs g,h,j), were completed in PdC-Jordan aquifers in the PdC Plateau, while the greatest percentage of decreasing wells, 15% (Fig. 9, ID f), were wells located above regional aquitards in the Galena Group.

Residence time and nitrate trend relationship

The relationship between mean residence time and nitrate trends from 27 wells and springs listed in Table 1 is displayed graphically in Fig. 10. Three wells, denoted with NA in Table 1, were not included due to insufficient data for trend analysis. Older groundwater had a higher proportion of increasing nitrate trends. Seven percent of the decreasing sites had a mean age of 15 years, while 59% of the increasing sites had a mean modern age of 30 years. All increasing trends were from the PdC, PdC-Jordan, Jordan-SL, and SL-LR on the PdC Plateau. All decreasing and no trend sites were from the Galena and Spillville aquifers on the UC Plateau.

Fig. 10
figure 10

Relationship between modern mean residence time and nitrate trends for 27 springs and wells listed in Table 1. Older groundwater had a higher proportion of increasing trend sites with all increasing sites located in deeper hydrogeologic settings of the PdC Plateau

Supporting analysis

In addition to residence time, other variables, including current and historical cropping, nitrogen fertilizer use, and precipitation trends, were examined to assist with trend interpretations (see ESM 2 for additional details and figures).

Cropping trends and nitrogen fertilizer use

Understanding current and historical cropping trends is important because the fraction of land area in row crops correlates to fertilizer inputs and nitrate loss to surface and groundwater (Schilling and Spooner 2006; David et al. 2010). Historical county-scale cropping trends based on survey data (USDA 2017) across nine Driftless Area counties in southeastern Minnesota show that cropland planted for corn and soybeans increased beginning in the 1960s (Fig. S2a of ESM 2). A particularly rapid increase of 8.8% annually was observed from 1960–1980, and a 4.4% increase for soybeans from 1980 to 2000. Increases in row-crop land cover largely displaced land previously devoted to small grains and perennial crops such as oats, hay, and alfalfa. From 2000 to 2020, combined corn and soybean cropland remained steady, while perennials continued to decrease. Cropland totals using survey methods are inferred and represent a fraction of the total area. For this reason, geospatial raster data were also utilized.

Geospatial raster data available for years 2010 to 2022 (USDA 2022) indicate a relatively slow annual increase in row cropping of about 0.6% year–1 across the Driftless Area of southeast Minnesota and northeast Iowa and decreases of perennials by 2%. Analysis of the two plateaus indicates the UC Plateau and PdC plateaus had row crop increases of 0.8% and 2.6%, respectively, for the 2010 to 2022 timeframe. For comparison, the glacial till landscape to the west showed no significant land-use changes. Analysis of the specific land area recharging to each study spring (Table S2 of ESM 2) revealed an average increase of 2.8% year–1, with a range of 0.5–5.0% for sites with statistically significant increases. Irrespective of location, the magnitude of row-crop change over the past decade is relatively small compared to the 8.8% annual change between the 1960 and 1980s (Fig. S2c of ESM 2).

Increases in land used for corn production largely correspond to increased sales and use of nitrogen fertilizers within the state of Minnesota (Fig. S2b of ESM 2). Assuming nitrogen sales in the study area are similar to statewide sales, rapid adoption and widespread use of manufactured commercial nitrogen fertilizers began in the late 1940s and increased rapidly from the early 1960s to ~1980 (20.7% year–1). This increase was concurrent with the rapid increase in corn and soybean production (Fig. S2c of ESM 2). Between 1980 and 2020, the rate of increase ranged between 1.1–1.8% year–1, with fluctuating nitrogen fertilizer sales. Since 2012, nitrogen sales in Minnesota have generally plateaued. Estimated sources of nitrogen from livestock manure are not included, so nitrogen inputs are considered conservative; however, nutrient-input patterns and variations are consistent with Falcone (2021).

Precipitation trends

Variable precipitation and weather can influence nitrate trends, especially in streams and aquifers with short residence times. Consecutive wet years can decrease nitrate concentrations, while repeated dry conditions followed by a wet year can increase the transport of nitrate (Randall and Mulla 2001). Hydrologic conditions in southeast Minnesota are characterized using the monthly Palmer hydrologic drought index (PHDI, Palmer 1965; Fig. S3 of ESM 2). Southeast Minnesota has been trending towards wetter conditions since the 1970s, with the highest recorded index values occurring between 2010 and 2020. Since records began in 1895, the five highest annual index values occurred consecutively between 2016 to 2020, with all months containing positive values. Additionally, an examination of precipitation trends during different time periods indicated precipitation was increasing by an average rate of 3.0 mm per decade (0.12 in.) from 2000 to 2011, while years 2012 to 2021 were much wetter, with an average rate of 39.6 mm (1.56 in) per decade (DNR 2023b).

Conceptual model of groundwater residence time

A conceptual model of residence time in the Driftless Area of southeast Minnesota is presented in Fig. 11. Carbonate rocks that form the UC and PdC plateaus are karstified, allowing rapid recharge to the underlying water table. Rapid flow through karstic fractures and conduits, coupled with a very low matrix porosity, results in relatively short residence times of about one to two decades. The UC Plateau is underlain by a regionally extensive aquitard that includes approximately 12 m of shale (Decorah Shale), which greatly limits the downward migration of young anthropogenically-impacted waters. Groundwater age beneath this aquitard is predominantly pre-1953. In the PdC Plateau, the relatively leakier Oneota aquitard is the uppermost geological unit, which limits vertical recharge. The Jordan Sandstone, directly underlying the PdC Plateau, is variably impacted by nitrate leaking through the Oneota aquitard, with a range of residence times of two to five decades. Where streams dissect the plateaus, the integrity of underlying aquitards near valleys is diminished and can result in ‘stair-step’ subsurface leakage of contaminants to underlying aquifers (Runkel et al. 2018; Barry et al. 2018). Incision is especially pronounced across the PdC Plateau, which allows for the mixing of the anthropogenically impacted groundwater from unconfined water-table aquifers with confined and much older, pristine groundwaters devoid of anthropogenic contamination.

Fig. 11
figure 11

Conceptual model of residence time (RT) for springs and wells in southeast Minnesota’s Driftless Area. Mean residence time, nitrate concentration, and depth by hydrogeologic setting for the 30 study springs and wells are displayed in bars to the left of the figure. Pre-1953 to millennia-aged water below regional aquitards schematically shown using results of tritium and carbon-14 dating. Letters A–C group 16 springs into three hydrostratigraphic groups and numbers 17–30 correspond to 14 study wells (see Fig. 2). Groups B and C springs are commonly under artesian conditions and depicted with upward pointing arrows

Discussion

Nitrate concentrations, trends, and residence time

The analysis revealed consistent patterns between residence times, nitrate concentrations and trends, and hydrostratigraphy, as depicted in Figs. 6 and 9 and the conceptual model in Fig. 11. The position of springs and wells in their hydrogeologic setting dictates the relative proportions of younger and older water, which has a substantial impact on characteristic nitrate concentrations and trends resulting in nitrate concentrations that are inversely related to groundwater age. This improved understanding of nitrate in springs and wells also provides insight into variable concentrations and trends measured in streams across the area. The following discussion expounds on these relationships, including an evaluation of the roles of agricultural land use and precipitation in impacting nitrate concentrations and trends of water with mixed residence times.

Short residence time springs: unconfined karstic systems of the UC and PdC plateaus

All shallow Galena karst springs of the UC Plateau and PdC Amherst Spring of the PdC Plateau have residence times of about one to two decades. Median nitrate concentrations for individual springs within this group range from 9–15 mg L−1, and are 50% higher than springs located in deeper stratigraphic settings. In addition, Galena springs are generally consistent with expected soil–water nitrate concentrations (Kuehner et al. 2020) from row crop agriculture. This indicates nitrate-laden water has extensively penetrated the aquifer and concentrations may be in approximate stasis with agricultural activities at the land surface over the past 10–20 years.

Galena springs feature greater variability in nitrate concentrations relative to the deeper, longer residence time springs (Fig. 9a). This variability is consistent with large short-term changes in continuous temperature measurements in Galena springs (Luhmann et al. 2011), as well as episodic fluctuations in flow, temperature, nitrate, and stable isotopes following recharge found at the comprehensively monitored Galena Group Bear Spring (55A406) in Olmsted County (Barry et al. 2020) and Crystal (23A028) and Twin Springs (23A059) in Fillmore County (Aggarwal et al. 2024; Kuehner MDA, unpublished data). The variability reflects the behavior of an unconfined karstic aquifer that is highly responsive (hours to days) to local event-scale rain and snowmelt, as well as annual-scale precipitation patterns. Subdued variability in nitrate concentration at Amherst Spring may reflect a longer residence time, due to its deeper hydrostratigraphic position, which somewhat buffers episodic recharge.

Moderate changes in row cropping within some of these springsheds over the past 13 years have not resulted in expected corresponding nitrate trends (Table S2 and Fig. S4 of ESM 2). The springsheds of Fountain East, Amherst, and Crystal Creek springs with younger groundwater age mixtures had increasing corn and soybean production and corresponding decreases of perennial land covers; yet, Fountain East did not have a statistically significant nitrate trend, and concentrations decreased at Amherst and Crystal Creek springs. Rainy Springshed had marginally decreasing corn and soybean production (p = 0.067) and increasing perennial covers, but there was no corresponding decreasing nitrate trend for this spring. Differences in residence time, combined with increases in precipitation over the past several decades, could be masking the impact of land-use changes—for example, consecutive years of increased precipitation, especially pronounced since 2010 (Fig. S3 of ESM 2), may be diluting nitrate concentrations and obscuring water quality signals at spring sites with a greater proportion of younger groundwater. Another possible explanation is more efficient fertilizer use resulting in less residual nitrate and leaching.

Longer residence time springs: incised valleys of PdC Plateau

The PdC-Jordan and SL-LR springs are located within valleys of the PdC Plateau (group B and C springs) and emanate from lower hydrostratigraphic positions (greater depths) beneath one or more aquitards, with groundwater residence times of two to four decades with significant contribution of confined water, and much lower nitrate (Fig. 9). The contribution of confined groundwater discharge to the Lanesboro SFH Main Spring is as much as 50% (Barry et al. 2023) and is expected to be even greater in deeper-positioned springs with longer residence times. This renders the comparison of nitrate concentrations and trends to recent row crop activity (Table S2 of ESM 2) problematic because much of the discharge from these deeper springs infiltrated the land surface decades to centuries ago. The lower nitrate concentrations in these deep springs are more likely a result of dilution from water that infiltrated the land surface prior to significant increases in nitrogen fertilizer use (Fig. S2b of ESM 1) as opposed to the lower percentage of recent row cropping in their steep forested springsheds. This is evidenced by disproportionately low nitrate concentrations relative to the amount of row cropping in the recharge area—for example, the springsheds of Cold South and Spring Creek have slightly greater row cropping than those of Group A Galena springs Fountain West and Rainy, yet have ~50% lower nitrate concentrations, respectively (Fig. 12).

Fig. 12
figure 12

Relationship between mean spring nitrate concentrations and percent of springsheds planted to corn and soybeans (row-crop) for 16 study springs grouped by hydrostratigraphic position. Nitrate concentrations and row cropping represent an average between 2010 and 2022. Hydrostratigraphic groups (A–C) are significantly different

All middle to lower PdC-Jordan and SL-LR springs exhibited statistically lower nitrate variability and increasing trends relative to the shallower springs (Fig. 9). The low variability is consistent with stable temperatures recorded at these springs, which indicates buffering from episodic recharge events. Increasing nitrate trends also suggest that nitrate from the land surface has not fully penetrated these deeper groundwater systems, which are still adjusting to current and historical nitrate loading rates. The estimated mean residence times of two to four decades suggest that much of the water contributing to these springs recharged when commercial nitrogen inputs and the amount of cropland grown for corn and soybeans would have been rapidly increasing (Fig. S2 of ESM 2). As the fraction of older, nitrate-poor water is displaced with younger, nitrate-enriched water, it is expected that nitrate concentrations will continue to increase. These conditions account for the high percentage of springs with increasing trends on the PdC Plateau.

In relatively deep parts of these aquifer systems, anoxic conditions that lead to the conversion of nitrate to nitrogen gas through denitrification may also lower groundwater nitrate concentrations, especially for Group C, SL-LR springs. Samples collected through the DNR Groundwater Atlas Program suggest the Lone Rock Aquifer is anoxic where it is deeply buried and oxic near valley edges and in valleys where it is shallowly buried (DNR 2023a). Any such denitrification, and the degree to which it diminishes concentrations, requires additional investigation.

Wells: differences and similarities to springs

Unlike spring discharge, which represents an integrated mixture of flow paths that can include variable proportions of young water susceptible to contamination, most well water is made up of a more discrete part of the groundwater system and is intentionally designed to limit contamination. In Minnesota, such well designs have been mandated since the adoption of well construction codes in 1974 (MDH 2011b). As a result, wells generally have a less direct connection to the land surface than springs, resulting in longer residence times (Fig. 7b) and lower nitrate concentrations (Barry et al. 2018). This discussion of well residence times, nitrate concentrations, and trends clarifies some of the similarities and differences to springs. This section also provides insight into how to interpret nitrate monitoring data from wells and nitrate transport to deeper aquifers over time.

The shortest estimated residence times among the study wells (Table 1) is two decades, which includes samples from two relatively shallow wells drawing from unconfined karstic conditions on the UC Plateau. This residence time is twice as long as UC Plateau Galena springs, and average nitrate concentrations (4.8–6.8 mg L−1) are about half that of the springs; however, shallow wells in this setting are like Galena springs in having more variable concentrations. Both study wells and 73% of the Galena wells in the trend wells dataset (Fig. 9) had no statistically significant trends, due to either high variability of nitrate concentrations or stasis in nitrogen loading.

Prairie du Chien Aquifer wells vary significantly in residence time and apparent responsiveness to agricultural practices and climate, which is expected for a thick (up to ~100 m), karstic interval. The five PdC study wells range in residence times from two to four decades, similar to the age range of PdC-Jordan springs, with increasing depth generally corresponding to longer residence time (Table 1). Relatively high and variable nitrate concentrations across the PdC Plateau are exemplified by trend well data that include wells with no record of grouting (Fig. 9). Ungrouted wells are more directly connected to the land surface relative to grouted wells, which allows for a greater contribution of recently recharged water and therefore will be more responsive to climatic variation and land use. In contrast, PdC wells that are known to be grouted have markedly lower and less variable nitrate concentrations (Fig. 9), which is characteristic of decades-old water and buffering from conditions at the land surface.

The three to four-decade modern residence times for Jordan-SL wells across the PdC Plateau (Table 1) are similar to the estimated ages from lowermost PdC-Jordan and SL-LR springs. Increasing trends in five of the seven Jordan-SL wells (Table 1, sites 24–30) and 9–13% of the wells from the trend wells dataset (Fig. 9b, sites g, h, j) are analogous to PdC Plateau springs with similar residence times and increasing trends. Water recharged decades ago, when nitrate inputs were rapidly increasing (Fig. S2b of ESM 2), is gradually migrating downward across the leaky Oneota Aquitard to the Jordan Aquifer at these wells. This increasingly impacted lower PdC, and Jordan groundwater contributes to discharge at deep springs in the PdC Plateau and thus to their characteristic increasing trends.

Groundwater to surface-water connections and precipitation variability

Springs and distributed baseflow are significant sources of cool isothermal water to many streams in southeast Minnesota; therefore, any long-term contaminant trends observed in springs are likely to manifest in streams. Stream trend analysis showed differences across the landscape (Figs. 8b and 9b), which can be partly explained by short or long residence times and variable precipitation and weather. Watersheds in the till-dominated landscape and UC Plateau, where groundwater residence time is short because of tile drainage or karstic conditions, show more variable nitrate trends by watershed and have a high percentage of streams with no trend. This is suggestive of baseflow that is highly responsive to annual and interannual weather variability, which may challenge the ability to interpret stream nitrate trends (Danalatos et al. 2022). In these more responsive watersheds, a wet period from 2012 to 2021 showed more stream sites with no nitrate trend or decrease, when compared to a drier 2000 to 2011 period (Figs. S1b and S1c of ESM 2). Decreases in nitrate occurred at four glacial till and four UC Plateau watersheds. Row crop activity across these watersheds was either stable (glacial till landscape) or slightly increasing (UC Plateau) during the 2012 to 2021 wet period. The decreasing nitrate trends, therefore, are most likely in response to the consecutive wet periods, especially between years 2016 and 2019 (Fig. S3 of ESM 2), which flushed residual soil nitrate stored within row crop soils and resulted in dilution of subsurface tile nitrate concentrations transported to streams. It is also possible that the consistent pattern of decreasing or no trends in the till landscape, and some areas of the UC Plateau, are signaling improvements in nitrogen fertilizer management and use efficiency. Improvements in nitrogen use efficiency (NUE), or the amount of corn produced per unit of nitrogen fertilizer applied, have been documented (MDA 2019), but finer-scale farm practice data and analysis are needed to evaluate how this may reduce nitrate loss at various scales in the context of precipitation variability.

In contrast, all but one PdC Plateau stream had an increasing trend from 2012 to 2022 (Fig. S1 of ESM 2). Although row cropping was modestly increasing (2.6%) across this area at the same time, residence time to these streams and springs in deeply incised valleys obscures the signal of variable precipitation and land-use changes. A significant component of baseflow to these streams is derived from confined groundwater discharge, with longer residence times and increasing nitrate concentrations that reflect land-use practices from decades past. Upward hydraulic gradients towards streams from these deeper aquifers carrying such water have been previously documented across the PdC Plateau (Barry et al. 2023).

The relationship between row crop activity and stream nitrate concentrations is also impacted by variability in baseflow from older water. Watersheds with short groundwater residence times in glacial till and UC Plateau settings have relatively high correlation slopes of 0.17 and 0.15 mg L−1 change for every percentage point change in row cropping, respectively (Fig. 13). These values are close to the 0.195 mg L −1 estimate by Schilling and Spooner (2006) for a till dominated watershed in south-central Iowa. In stark contrast, watersheds with the longest residence times in the study area, those within the PdC Plateau, have significantly lower nitrate and correlation slopes that are only 0.05 mg L−1 or 67% less, reflecting the greater contribution of old, confined groundwater with low nitrate concentrations.

Fig. 13
figure 13

Relationship between mean stream nitrate concentrations and percent of watershed area planted to corn and soybeans (row-crop) by geologic setting. Nitrate concentrations and row cropping represent an average between years 2000 and 2021 across 55 watersheds in south-central and southeast Minnesota and northeast Iowa

Best management practice implications

Recognition of water quality vulnerability in the Driftless Area has prompted the acceleration of agricultural BMPs to reduce deleterious impacts on water quality. Quantifying residence time alone, however, may not entirely answer the question of how long it may take to detect water quality improvements of BMPs on the land surface (Meals and Dressing 2008). Groundwater discharging from a spring or pumped from a well is a composite of converging groundwater flow paths of differing groundwater ages (Focazio 1998), and it is the distribution of ages that influences water quality response (Lindsey 2003). Measuring water quality response at drinking water wells may be tenuous since they are intentionally constructed to avoid recent waters—for example, the estimated fraction of pre-1953 groundwater for well 777,654 is an estimated 61%, resulting in a total mean age of 78 years (Table 1).

Although drinking water wells are important in documenting current contamination levels that may impact human health and threaten ecosystems, many monitoring sites with residence times of several decades or more are unlikely to be useful for determining the efficacy of current agricultural BMPs. Just as the rising alachlor ESA levels characteristic of many deep springs (Figs. 4a and 6a) are reflective of pesticides applied decades ago, rising nitrate levels are reflective of fertilizer application from decades past (Fig. S2 of ESM 2). The effect of lag time applies to other contaminants, such as cyanazine, an herbicide that was used on corn crops from the early 1970s through the 1990s. Its usage was voluntarily stopped in 2002, but breakdown products in groundwater are still present in southeast Minnesota (MDA 2023, 2024). Decadal time lags have important implications for groundwater protection programs, policies, and managing expectations with water quality stakeholders. Moreover, decadal time lags invalidate many of the assumptions used by existing groundwater BMP prediction models (Ascott et al. 2021).

Measuring water quality response within short time frames is more probable in shallow springs, wells, and streams of the unconfined karstic aquifers on the UC and PdC plateaus, where residence times can be as little as 10 years—for example, a rapid water quality response was observed in a shallow (15.6 m, 52 ft) Galena monitoring well (695,883) located on the UC Plateau due to a sudden and dramatic land-use change (Fig. 14, Fig. S2 of ESM 2). From 2000 through 2005, the percentage of row cropping around the well decreased rapidly from ~80 to <2%, whereas between 2004 and 2022, nitrate concentrations declined significantly by 77%, from 17.6 to 4.1 mg L−1 over the 18-year period, with the most significant reductions occurring during the first 7 years. Trend interpretations for these shallow systems must also account for the influence of precipitation variability and the scale of land-use change or BMP adoption. Monitoring wells intentionally constructed to capture recent water can best provide the desired responsiveness if BMP adoption is consistent and substantial.

Fig. 14
figure 14

Rapidly decreasing nitrate–N concentrations after a sudden and substantial land-use change on the UC Plateau from predominantly agricultural to high-density housing in the catchment area of a shallow well (695,883)

Dominant factors affecting nitrate concentration trends in groundwater and surface water

Using the conceptual model and monitoring data presented in this report, Fig. 15 summarizes the dominant factors influencing nitrate concentration trends in groundwater and surface water. Dominant factors affecting nitrate trends above regional aquitards (above the black line in Fig. 15) that contain a greater proportion of local and younger groundwater mixtures (<20 years) include precipitation/weather variability, recent land-use changes, and recent nitrogen inputs and management. Alternatively, in areas below regional aquitards with groundwater mixtures that range from several decades to millennia, nitrate trends are largely controlled by dilution from regional and deeper groundwater mixtures, historical land-use changes, and historical nitrogen inputs.

Fig. 15
figure 15

Dominant factors affecting nitrate concentrations and trends in the Driftless Area of southeast Minnesota

Conclusions

The convergence of parallel results across multiple independent age tracer methods greatly improves the confidence of modern (post-1950s) residence time estimates in the karstic, multiaquifer system of southeast Minnesota. Results from this analysis indicate that groundwater from many springs and wells in southeast Minnesota have modern ages that typically range from 10 to 40 years. Residence time is correlated with depth and implies nitrate and alachlor ESA present within much of the shallow water-table aquifer, have not fully migrated into deeper aquifers where residence time mixtures can range from several decades to millennia.

Residence time estimates represent the integrated averages of the duration of time from when water infiltrates the land surface to its discharge at a spring or well. Residence time is greatly influenced by hydrostratigraphy, which leads to variable mixing of disparate water ages. Kraft et al. (2004) reasoned that nitrate will eventually reach a steady state in aquifers when modern groundwater with elevated nitrate penetrates the thickness of the entire aquifer. Due to groundwater mixtures with components that date back several decades, it is therefore predicted that nitrate measured in certain springs, streams, and wells, especially in the PdC Plateau, will continue to increase until an equilibrium has been reached with historical and modern land-use practices. Due to shorter residence times, nitrate concentrations from Galena and shallow PdC springs and wells located above aquitards are in approximate stasis and likely to have a faster response to land-use practices, especially when nitrogen input reductions are substantial. There are positive indicators that sites with high nitrate concentrations are relatively stable or in decline in these shallow aquifer settings, which may reflect improved practices. The influence of above-normal precipitation patterns resulting in a dilution effect is another key factor.

Residence time is an essential concept that resource managers and stakeholders should consider when calibrating groundwater models, interpreting groundwater monitoring trends, evaluating the effectiveness of BMPs, and defining realistic management goals and expectations. Separating the effect of current crop and fertilizer management practices on groundwater quality will require consideration of response time lags, precipitation/weather variability, and land use. This highlights the value of documenting the efficacy of nitrate management practices using specialized monitoring efforts, apart from those typically employed at springs, streams, and wells that can be confounded by varying residence times. Surrogate metrics, such as more detailed accounting of nitrogen use efficiency and row crop acres utilizing nitrogen BMPs, paired field-scale vadose zone monitoring networks, and improved groundwater models that integrate residence time estimates by hydrogeologic setting, may help improve water quality performance measures in the near term. In the long term, these results underscore the importance of efforts to minimize nitrate leaching at the field scale since current reductions will help reduce future loading to deeper aquifers. Further, this analysis underscores the value of maintaining long-term groundwater pesticide and nutrient monitoring programs and the ability to leverage these data for additional insights and applications.