Natural and Anthropogenic Geochemical Tracers to Investigate Residence Times and Groundwater–Surface-Water Interactions in an Urban Alluvial Aquifer

Abstract

A multi-component geochemical dataset was collected from groundwater and surface-water bodies associated with the urban Fountain Creek alluvial aquifer, Colorado, USA, to facilitate analysis of recharge sources, geochemical interactions, and groundwater-residence times. Results indicate that groundwater can be separated into three distinct geochemical zones based on location within the flow system and proximity to surface water, and these zones can be used to infer sources of recharge and groundwater movement through the aquifer. Rare-earth-element concentrations and detections of wastewater-indicator compounds indicate the presence of effluent from wastewater-treatment plants in both groundwater and surface water. Effluent presence in groundwater indicates that streams in the area lose to groundwater in some seasons and are a source of focused groundwater recharge. Distributions of pharmaceuticals and wastewater-indicator compounds also inform an understanding of groundwater–surface-water interactions. Noble-gas isotopes corroborate rare-earth-element data in indicating geochemical evolution within the aquifer from recharge area to discharge area and qualitatively indicate variable groundwater-residence times and mixing with pre-modern groundwater. Quantitative groundwater-residence times calculated from ³H/³He, SF₆, and lumped-parameter modeling generally are less than 20 years, but the presence of mixing with older groundwater of an unknown age is also indicated at selected locations. Future investigations would benefit by including groundwater-age tracers suited to quantification of mixing for both young (years to decades) and old (centuries and millennia) groundwater. This multi-faceted analysis facilitated development of a conceptual model for the investigated groundwater-flow system and illustrates the application of an encompassing suite of analytes in exploring hydrologic and geochemical interactions in complex systems.

1. Introduction

Groundwater is a major source of water used for domestic and industrial purposes worldwide, particularly in arid and semi-arid regions. Alluvial aquifers are important in many such regions because of their generally high transmissivity, correspondingly large production rates, and accessibility for public-water supplies. Because many alluvial aquifers are unconfined and located along present-day streams, they also can be affected by land use and anthropogenic chemical sources [1], possibly necessitating treatment of groundwater prior to use as a drinking-water supply. In addition to contaminants of anthropogenic origin that have been studied for decades, such as NO₃, contaminants of emerging concern (CECs) are becoming increasingly identified in alluvial aquifers and other settings [2,3]. These CECs can include fuel components, solvents, pharmaceutical and wastewater-indicator compounds [2,4], and per- and polyfluoroalkyl substances (PFAS) [5,6].

There are numerous techniques for evaluating the occurrence and distribution of solutes in the environment, including spatiotemporal trend analysis [7], stable isotopes of water and nitrogen [8], and statistical methods [9,10]. Although each of these methods provides insight into geochemical processes occurring in an aquifer, additional knowledge from environmental tracers can provide independent evidence to define recharge conditions, groundwater–surface-water interactions, groundwater mixing, and residence times [11,12,13,14], all of which help to control solute distributions in an aquifer.

This study used a large geochemical dataset (approximately 200 total parameters) to characterize recharge and solute sources, geochemical reactions, flow paths, and groundwater–surface-water interactions in an alluvial aquifer in Colorado, USA. This study was proposed, in part, because of the occurrence of PFAS compounds in the Fountain Creek alluvial aquifer [5]. Geochemical constituents used in the analysis include major ions, trace elements, rare earth elements (REEs), stable isotopes, noble-gas concentrations and isotopes, and pharmaceuticals and wastewater-indicator compounds. Recent studies applying subsets of analytes used in this study have highlighted the benefits of using various tracers in complex hydrologic systems [7,14]. This study builds on previous work, by combining analyses of multiple types of tracers that have not yet been utilized together in the literature (e.g., REEs, noble-gas isotopes, and pharmaceuticals). For example, Reference [15] analyzed noble gases, nitrate, and stable isotopes of water, to investigate sources of water to wells in a mountainous environment, whereas Reference [16] used pharmaceutical compounds and tritium content across the United States, to evaluate occurrence and distribution of anthropogenic compounds in groundwater, as did Reference [7], using REEs in a near-stream aquifer. However, previous studies, such as those noted, used less comprehensive datasets than those described here (e.g., Reference [7] does not include an evaluation of groundwater age). This study is novel in that a set of tracers is applied which allows a more complete understanding of complex aspects of this hydrologic system including the interaction and mixing of groundwater and surface water with a range of ages and compositions. Understanding the ranges in groundwater residence times is critical in this study area because of transient histories of potential contaminant loading from several distinct sources. One end goal of this analysis is to further inform ongoing groundwater modeling in the area, and the results provide an example of the benefit of an integrated approach to investigating sources of recharge for an alluvial aquifer.

The study area is the Fountain Creek alluvial aquifer, located near the city of Colorado Springs, Colorado, USA (Figure 1), on the eastern margin of the Southern Rocky Mountains. The alluvial aquifer is primarily composed of Quaternary-age sand and gravel fluvial sediments derived from weathering of Precambrian granitic rocks to the west and Mesozoic sedimentary rocks to the east [17,18]. The unconfined aquifer has been the subject of previous studies that focused on sources of solutes in the interconnected groundwater–surface-water system, as well as on groundwater-flow directions and water budgets [17,18,19,20,21,22]. Previous studies generally found that NO₃ was a constituent of concern and was derived from a variety of sources including fertilizers, effluent from wastewater treatment plants (WWTPs), and septic systems [8]. Work by Reference [22] indicated the influence of WWTPs on solute concentrations and flow volume of Fountain Creek and determined that resource extraction activities in the area contributed dissolved metals (Cu, Mn, Ni, Pb, and Zn) to the stream. In addition to anthropogenic compounds, the aquifer has naturally elevated concentrations of Se and U derived from the Cretaceous-age Pierre Shale, which underlies the alluvial aquifer [17]. Recently, concentrations of PFAS compounds greater than US Environmental Protection Agency (EPA) lifetime health advisory level of 70 nanograms per liter (ng/L) [23] have been found in the alluvial aquifer [5]. Potential sources of PFAS compounds include WWTPs, aqueous-film-forming foam (AFFF) used in firefighting operations, and metal-plating operations [6,10], each of which is documented in the area.

Hydrologically, Fountain Creek interacts with the shallow unconfined groundwater system by both gaining and losing flow depending on the relative elevation difference between the stream and the water table [17,20]. Effluent from WWTPs contributes a substantial portion of flow to the creek during base-flow periods [22]. In the study reach, there are diversions from the creek for agricultural purposes, and these diversions, as well as streamflow, are simulated by using a transit-loss model [29], which is used by water managers for water-supply planning [30]. Streamflow in the area generally peaks in late spring and early summer in response to snowmelt runoff from the mountainous areas to the west (based on streamflow measured at US Geological Survey site number 07105500, Fountain Creek at Colorado Springs, Colorado) [24]. Groundwater flow within the central portion of the alluvial aquifer is generally from north to south consistent with the topographic gradient and streamflow direction [18]. Several small paleochannels are saturated to the northeast of the central portion of the aquifer (and distant from Fountain Creek). These paleochannels have been termed the tributary alluvium and are hydraulically and topographically upgradient from the main valley of the Fountain Creek aquifer (Figure 1). Groundwater flow between the tributary alluvium and the main aquifer appears to be contiguous based on water-level observations (see Supplemental Materials).

The study area is semi-arid with a median annual precipitation rate of approximately 430 millimeters per year, determined based on gridded meteorological datasets based on the period from 1990 to 2020 [31]. Precipitation is seasonal with the greatest quantity generally falling in the summer as monsoon rainstorms. Because of the semi-arid climate, precipitation recharge is minimal, and much of the groundwater recharge is derived from streamflow losses and irrigation return flow [17,18]. Groundwater pumping in the study area has occurred for decades for domestic, industrial, and agricultural uses [20,32] with combined pumping withdrawals up to 19,000 cubic meters per day [20]. Groundwater pumping of this magnitude for prolonged periods of time has likely altered the flow field creating localized groundwater depressions around pumping wells. Patterns in groundwater-level elevations are complicated by pumping in the area, but wells with less pumping-induced disturbances display seasonal variations controlled by snowmelt runoff recharge, streamflow seasonality, and precipitation seasonality [24].

2. Materials and Methods

2.1. Sample Collection, Analysis, and Quality Control

Water-quality samples were collected from groundwater wells and surface-water features, illustrated in Figure 1, to characterize aquifer conditions, as well as the composition of end-member sources of recharge (i.e., recharge sources with distinctive geochemistry). Groundwater wells sampled for this study included monitoring, public-supply, agricultural-supply, and domestic-supply wells. Groundwater sites were selected based on location within the alluvial aquifer and accessibility. Where possible, dedicated monitoring wells from previous investigations [21] were selected for sampling. Monitoring wells were sampled, using either a portable Grundfos stainless steel pump or a Proactive Mini-Monsoon 12V low-flow pump, and pumps were thoroughly cleaned between wells following standard US Geological Survey (USGS) techniques [33]. Where dedicated monitoring wells were not available, samples were collected from public-supply, agricultural-supply, or domestic-supply wells. For wells with dedicated pumps (i.e., pumps that remain in place), samples were collected from available sample ports prior to any water treatment, thus providing complementary data to samples collected by using removable pumps which were also not subjected to any treatment prior to sample collection. Surface-water sampling sites were selected to characterize the surface-water chemistry of Fountain Creek, as well as Canal 4, which diverts water from Fountain Creek to eastern parts of the study area (Figure 1). Both Fountain Creek and Canal 4 are considered potential sources of groundwater recharge to the alluvial aquifer.

All samples were collected according to procedures described in the USGS National Field Manual [33]. Briefly, groundwater wells were pumped at a low rate until temperature, specific conductance, pH, dissolved oxygen (DO), and turbidity had stabilized to within respectively +/−0.2 ℃, +/−3%, +/−2 standard units (SU), +/−0.3 milligrams per liter (mg/L), and +/−0.5 turbidity units, or until three casing volumes were removed, whichever took longer. This ensured that water being sampled was representative of water within the aquifer. Surface-water samples were collected, using methods for isokinetic sampling [33], although necessary velocity to achieve true isokinetic conditions was not obtained in all instances. Alkalinity and pH were measured in the field for groundwater samples and at the USGS National Water Quality Laboratory (NWQL) in Lakewood, Colorado, for surface-water samples [34].

Samples for analysis of major and trace elements, pharmaceuticals, and wastewater-indicator compounds, and stable isotopes of NO₃ were collected in polyethylene bottles and filtered to 0.45 micrometers (µm). Nitrogen-isotope samples were then filtered a second time, using a 0.20 µm filter and frozen. Samples for metals analysis were preserved in the field with nitric acid (HNO₃) and refrigerated, and samples for anions were unpreserved but refrigerated. Samples for stable isotopes of water were collected, unfiltered, in glass bottles. Quality-control samples (field blanks and replicates) were collected for trace elements, pharmaceuticals and wastewater-indicator compounds, and stable isotopes of nitrate (replicate only, no blank) to evaluate sampling and/or analytical bias. Tritium (³H) samples were collected, unfiltered, in polyethylene bottles with no headspace by bottom filling after flushing with three volumes of sample water. Samples of the major dissolved gases O₂, N₂, Ar, CO₂, and CH₄ were collected, unfiltered, in replicate, in glass bottles sealed with no headspace by bottom filling. Samples of noble gases (He, Ne, Ar, Kr, and Xe) were collected, unfiltered, in replicate, in copper tubes [35] and flushed of all air bubbles before being sealed under positive pressure. Samples of sulfur hexafluoride (SF₆) were collected, unfiltered, in replicate, in glass bottles with no headspace by bottom filling after flushing with three volumes of sample water.

Anions were analyzed by ion-exchange chromatography, and cations were analyzed by inductively coupled plasma-atomic emission spectrometry at the NWQL [36,37]. Trace elements were analyzed by inductively coupled plasma-optical emission spectrometry and inductively coupled plasma-mass spectrometry at the USGS Analytical Trace Elements Chemistry Laboratory located in Boulder, Colorado [38]. Pharmaceutical compounds and wastewater-indicator compounds were analyzed by high-performance liquid chromatography tandem mass spectrometry [39] at the NWQL. Stable isotopes of nitrate (δ¹⁵N, δ¹⁸O_NO3) were analyzed by continuous-flow isotope-ratio mass spectrometry (IRMS) [40] at the USGS Reston Stable Isotope Laboratory (RSIL) located in Reston, Virginia, and are reported in standard stable-isotopic units of per mil (‰) [41] with uncertainties of +/−0.5‰. Stable isotopes of water (δ²H, δ¹⁸O_H2O) were analyzed by dual-inlet IRMS [42,43] at the RSIL with uncertainties of +/−2‰ and +/−0.2‰, respectively. Tritium analysis was conducted by using distillation and electrolytic enrichment, followed by liquid scintillation, at the USGS Menlo Park Tritium Laboratory, located in Menlo Park, California, and are reported in tritium units (TU) with analytical uncertainty ranging from 0.09 to 0.17 TU. Concentrations of major gases were analyzed, using gas chromatography at the USGS Reston Groundwater Dating Laboratory located in Reston, Virginia, with uncertainties of 0.002 mg/L (O₂), 0.001 mg/L (N₂), 0.003 mg/L (Ar), 0.04 (CO₂), and 0.0005 mg/L (CH₄). Concentrations and isotopes of noble gases were analyzed at the USGS Noble Gas Laboratory located in Denver, Colorado, using a magnetic-sector mass spectrometer and ultralow vacuum extraction line [44], and are reported in units of cubic centimeters at standard temperature and pressure per gram of water (ccSTP/gH₂O) and as ratios (e.g., ²⁰Ne/²²Ne), respectively, with analytical uncertainties of 1% (He), 2% (Ne), 2% (Ar), 3% (Kr), and 3% (Xe). Concentrations of SF₆ were analyzed at the USGS Groundwater Dating Laboratory, using purge and trap gas chromatography followed by an electron capture detector according to the method of Reference [45] and are reported in units of femtomoles per kilogram (fm/kg) with uncertainty ranging from 20% at the minimum reporting limit to 3% at the maximum reporting limit (https://water.usgs.gov/lab/sf6/lab/analytical_procedures/, accessed on 1 January 2021). All water-quality data are available from the USGS National Water Information System (NWIS) database [24], using the USGS site IDs provided in Supplementary Materials Table S1.

Data quality was evaluated through quality-control samples, calculation of charge balance for major-ion data, and comparison to historical data. This dataset included more than 10 percent field blanks and replicates. Source solution, equipment, and field blanks were collected to test for possible contamination [46]. For the blanks, there were few detections for trace elements where only 7% of these detections were greater than 1 microgram per liter (µg/L). These detections greater than 1 µg/L occurred for Al, Cu, and Zn. The potential bias in samples for these elements was assessed by comparing the maximum concentration reported in blanks to the average environmental concentration. Potential bias of Cu and Zn are generally low given that the maximum concentration detected in the blanks are 6% and 16%, respectively, of the average environmental concentration. For Al, the maximum concentration in blanks was 103% of the average environmental concentration, indicating greater potential for bias and resulting in Al being removed from subsequent statistical analyses. Variability between environmental and replicate samples, quantified by the median of relative percent difference (RPD), for elements with potential bias (Al, Cu, Zn) were respectively 17%, 18%, and 9%, indicating that there is relatively little variability in the dataset. Of the 11 pharmaceutical and wastewater-indicator compounds detected in blanks, none had values above the reporting level. Rare earth elements did not have detections in the blanks with exception of 3 samples, which were greater than 0.002 µg/L for Ce. The maximum concentration of detections in blanks for Ce are 9% of the average environmental concentration, meaning samples are likely not substantially biased for Ce. RPD calculated from field replicates was less than 10 percent for nutrient and major-ion data with exception of 1 replicate pair that had an alkalinity RPD of 14%. Charge balance for the dataset ranged from −10.1 to 12%, where all but three samples had charge balances within the +/−10% range suggested for many water-quality studies [47]. Given the infrequent detections in blanks and generally low RPDs, the dataset is considered representative of environmental conditions in the study area.

2.2. Multivariate Evaluation

Laboratory analytical results for many constituents (particularly pharmaceuticals, wastewater-indicator compounds, and trace elements) were commonly below laboratory reporting limits (LRLs), which inherently limits data analysis [16]. In subsequent interpretive analyses, all censored values below the LRL were imputed, using regression-on-order (ROS) statistical methods [48], implemented in the NADA package of R [49]. The ROS method allows for robust estimation of censored values by using the structure of the dataset above the LRL and projecting this data structure to represent the censored concentrations. Results of imputation by ROS were further filtered by removing all parameters for which greater than 80% of the results were below the LRL (i.e., when more than 80% of the dataset was censored for a given parameter). This was done to limit the multivariate evaluation to those parameters with more distributed spatial information and to limit the effects of the ROS imputation on the overall results of the study.

Principal component analysis (PCA) was applied on the modified dataset to reduce the dimensionality of the dataset and aid in interpretation of geochemical signatures of different water types [50]. Geochemical data types were grouped together as follows: major ions, trace elements, REEs, noble gases, noble-gas isotopes, and pharmaceutical and wastewater-indicator compounds. Reference [51] suggests that multivariate compositional data analysis (CoDA) may not be useful for isotopes, so isotope ratios of noble gases were included in the PCA to test that assumption. The group-specific dataset was scaled, using the centered-log ratio transform (clr), which allows closed compositions (i.e., compositions summing to a fixed value [52]) to be reliably evaluated, using statistical methods [50]. There are a variety of transformations that could be applied (isometric-log ratio, ilr; allometric-log ratio, alr; centered-log ratio, clr), and Reference [52] suggests that the ilr may be superior for multivariate evaluation. However, application of the results of ilr with PCA may be difficult because back-transformation is required prior to interpretation [9,50]. As an alternative the clr was used, as it has been established as a reasonable method for scaling prior to PCA [9,51]. PCA was then applied, using the prcomp function of R [49].

Iterative PCA using different combinations of parameters may be helpful to further define groupings and understand processes [9]. Therefore, in addition to PCA of the previously defined geochemical groups, a combined PCA was also completed, using a subset of parameters. These parameters were included based on previous group-specific PCA, prior geochemical knowledge of the system, and other research in the Fountain Creek alluvial aquifer [17,22]. All parameters were converted to the same units of µg/L prior to regrouping because this is required by PCA [50]. Following unit conversions, the dataset was again processed using the clr prior to completion of PCA.

2.3. Noble Gases and Groundwater Residence-Time Tracers

2.3.1. Dissolved Gas Modeling

Estimation of conditions during groundwater recharge (specifically water temperature, atmospheric pressure, excess air, and salinity) are required for attributing proportions of dissolved gases ³H, He, and SF₆, which are used in lumped parameter models [53] and apparent-age calculations [54,55]. Inversion of dissolved-gas concentrations are commonly used to estimate recharge conditions, and the inclusion of noble gases in addition to major gases allows for greater confidence because of the specific solubility characteristics of noble gases [35,56].

In this analysis, Ne, Ar, Kr, Xe, O₂, and CO₂ were used because they are primarily of atmospheric origin. Helium and N₂ may have non-atmospheric sources [8,12] and were only used when no viable models were obtained from the set of gases with primarily atmospheric origin. Helium and N₂ were used in six of the 19 models for sampled groundwater wells where noble-gas data were available. Although Ar can have terrigenic sources, use in gas modeling is reasonable if samples do not show terrigenic ⁴⁰Ar production [35]. Evaluation of the ratio of ³⁶Ar/⁴⁰Ar in samples versus that in the air (³⁶Ar_air/⁴⁰Ar_air = 0.003378) [35] indicates that subsurface radiogenic production of ⁴⁰Ar can be ignored (i.e., mean ³⁶Ar_sample/⁴⁰Ar_sample / ³⁶Ar_air/⁴⁰Ar_air = 0.998).

For the nine groundwater wells that were not sampled for noble gases, recharge temperature and excess air were computed by the USGS Groundwater Dating Laboratory from concentrations of O₂, N₂ (corrected for excess gas), Ar, and CO₂ according to methods described in References [57,58] and References [59,60].

Modeling with noble-gas data utilized an inverse approach similar to those described by References [61,62]. Input parameters included the recharge salinity from groundwater sampling and the elevation of the land surface in feet above the North American Vertical Datum of 1988 (NAVD88) at the well. Using these inputs, an inverse-error fitting procedure was applied wherein sets of simulated gas concentrations are produced and compared to the observed values, for each of the three most commonly used models of excess air: closed-system equilibration (CE), partial re-equilibration (PR), and unfractionated excess air (UA), as described by Reference [35]. The error in these sets of simulations is quantified by using the χ2 parameter, which accounts for model errors in each gas. The model returns a recharge temperature, excess air concentration, χ2 value, and probability for each of the CE, PR, and UA models. Recharge temperature and excess air are crucial for groundwater-age estimates discussed in the next section, and in of themselves, are useful for hydrologic characterization [63].

A sensitivity analysis was conducted for several wells to evaluate the influence of changing recharge elevation on calculated recharge characteristics. Three different situations were simulated: (1) recharge elevation set equal to wellhead elevation (minimum elevation possible; this is the base-case model described previously), (2) recharge elevation set equal to the median elevation within the watershed, and (3) recharge elevation set equal to the maximum possible elevation within the upgradient watershed. Results of the sensitivity analysis indicate that simulated recharge conditions are relatively insensitive to the input recharge elevation (see results in Supplementary Materials).

2.3.2. Apparent Ages and Lumped-Parameter Modeling

Concentrations of ³H, ³He, and SF₆ in groundwater were used as tracers to evaluate the presence of modern water [55,64], as well as to estimate residence-time distributions [65]. Tracer concentrations were normalized to the groundwater-recharge conditions (temperature and excess air) calculated from dissolved gases as described previously.

A full discussion of the complexities and nomenclature relating to estimation of the time since groundwater recharge has occurred (generally referred to as “groundwater age”) is not discussed herein but is provided by References [66,67]. Generally, groundwater age is conceptualized as having a distribution rather than a discrete value because of processes such as mixing and dispersion. Methods were applied in this work to estimate age distributions of modern groundwater (i.e., groundwater recharged since the 1950s). For consistency, the terminology of Reference [67] is used with respect to apparent tritium-helium age (³H/³He), idealized age (for SF₆), and mean-residence time (for lumped-parameter models). Apparent ³H/³He ages were calculated according to methods described in References [54,56]. The apparent ³H/³He age relies on the known decay rate of ³H (derived from atmospheric nuclear weapons testing) to ³He, combined with an understanding of He sources in the aquifer [54]. Piston-flow ages (i.e., idealized ages) were calculated from SF₆ concentrations in water samples compared to SF₆ concentrations in the atmosphere through time based on methods described in Reference [45]. Finally, residence-time distributions were evaluated, using the program TracerLPM [53] for the piston-flow model (PFM), exponential-mixing model (EMM), and dispersion model (DM). These models are consistent with simplified understanding of the groundwater-flow system based on aquifer thickness, groundwater-flow directions, and residence-time distribution functions [68].

In addition to quantitative estimates of groundwater age, the fractions of old and young groundwater can be calculated, using ³H data alone (herein referred to as the ³H thresholds method), which may be less ambiguous than application of lumped parameter models (LPMs) because the conceptual mixing relations need not be known [69]. The ³H thresholds method [64,69] was applied to calculate the fraction of water recharged since the advent of atmospheric nuclear testing in 1953 (F_post-1953), using the measured ³H concentrations and estimated ³H in precipitation [70,71,72]. The ³H threshold method allows three categories of groundwater to be differentiated; modern (recharged after the 1950s), pre-modern (recharged prior to the 1950s), or mixed. The method consists of comparing decayed concentrations of ³H in precipitation for the study area to the ³H concentrations in water samples. Based on different thresholds that are established, the sample is then categorized into one of the three groups. Differences between the precipitation input datasets of References [70,71] are assessed in Supplementary Material Section S4.

3. Results and Discussion

All data described in the following sections are available through the USGS NWIS database [24]. Results are compiled according to USGS site numbers and site names, provided in Supplementary Materials Table S1. Groundwater wells were categorized into three groups based on their spatial location (Figure 1). Tributary alluvium wells are in the paleochannels to the northeast and upgradient from the main aquifer, wells proximal to the creek are within approximately 200 meters of Fountain Creek, and the remainder of the wells were classified as central groundwater.

3.1. Major Ions and Trace Elements

The major-ion composition of groundwater and surface water covers a broad range but also shows some spatial associations. Based on the Piper diagram [73] for the study area surface water is typically Na+K-Cl type, whereas groundwater is typically Ca-SO₄ type or Ca-HCO₃ type (Supplementary Materials Figure S1). Fountain Creek and Canal 4 have appreciable Cl possibly indicative of road salting effects. The composition of central groundwater (Supplementary Materials Figure S1) is between tributary alluvium and surface water, indicating a mixed origin. Ratios of Cl/Br and concentrations of Cl in groundwater may be used to indicate mixing and inputs from anthropogenic sources [74,75,76]. An analysis of Cl/Br ratios and Cl concentrations (Supplementary Materials Figure S2) indicates that groundwater well 04-009, located in the tributary alluvium (Figure 1), is representative of native groundwater recharge and that the majority of groundwater at other wells can be explained by variable mixing between native groundwater recharge with WWTP effluent, septic tank discharge, and road salting. Similarly, surface waters have compositions that include WWTP effluent, road salting, and meteoric precipitation. This analysis indicates that mixing relations between surface water, groundwater, and anthropogenic sources are important for understanding solute distributions in the aquifer. These results are further corroborated and strengthened by REE data and PCA, as described below.

Oxidation-reduction (redox) classification of the waters, using the method of Reference [77] indicates that 26 of 31 groundwater sites are characterized by oxic conditions whereas the remaining five display mixed oxic-anoxic conditions (Supplementary Materials Table S2). Mixed conditions of this type are possibly consistent with NO₃ reduction or Mn reduction [77].

Trace-element (Al, Cu, Fe, Mn, Pb, and Zn) data show spatial patterns in groundwater and surface water (Supplementary Materials Figure S5) potentially linked to variable processes or sources controlling transport, though Al may be marginally affected by sample bias and is thus not interpreted further. For instance, Cu, Fe, Mn, and Zn tended to have greater concentrations in surface water than in groundwater (Supplementary Materials Figure S5), possibly attributed to runoff from of mine tailings in the area as documented by Reference [22]. In contrast, Pb tended to have greater concentrations in groundwater potentially linked to a subsurface source or redox conditions in the aquifer.

3.2. Stable Isotopes

The δ²H and δ¹⁸O_H2O composition of groundwater (Figure 2) is generally consistent with mixing of water sources from a snowpack-driven Rocky Mountain meteoric water line (RMWL-Snow) as described by Reference [78] and annual variation in precipitation for the site estimated by using the online isotopes in precipitation calculator (OIPC) [79]. The composition is shifted slightly to the right of the snow-dominated sites likely because the area receives both snowpack accumulation and recharge from summer rainstorms (as reflected in the relative positions of the data and the OIPC meteoric water line (MWL). Tributary alluvium groundwater wells have a large range in isotopic composition, covering nearly the entire spread in the dataset (Figure 2a). This range is most likely related to seasonality of recharge, though the datasets required to thoroughly analyze such trends are not available for the study site (e.g., see Reference [80]). There are no clear spatial associations of wells with more depleted versus enriched compositions. Surface-water samples have a slightly shallower slope (slope = 7.16) versus groundwater (slope = 8.19), consistent with some evaporation of surface water [81].

Stable isotopes of nitrate (δ¹⁵N, δ¹⁸O_NO3) provide evidence that denitrification is occurring in the aquifer, as several wells plot along a line sloping at 1:2 (δ¹⁵N:δ¹⁸O_NO3) in Figure 2c, a trend typical of denitrification [41,82]. These wells also show ammonia and nitrate + nitrite (NO₃ + NO₂) concentrations consistent with isotopic fractionation by microbial reduction (Figure 3). Several groundwater wells and streams plot in a region with generally elevated NO₃ + NO₂ and ammonia (circled sites on Figure 3b). These wells are in the general proximity to current or former WWTPs (Figure 1), and thus the coincidence of multiple elevated nitrogen species is likely linked to WWTP effluent [22]. It is interesting to note that the indication of some nitrate-reducing conditions by N isotopes and concentrations of NO₃ + NO₂ and ammonia is slightly inconsistent with primarily oxic conditions as indicated by the method of Reference [77]. Specifically, of the wells in Figure 3 noted to have stable-isotopic evidence of denitrification, only well U-7 was noted to have mixed oxic/anoxic redox conditions by the method of Reference [77]. Isotopes of N and O are ideally suited to identifying nitrate reduction [41,83] whereas the idealized method of redox-condition assignment is subject to a variety of limitations [77]. The prevalence of oxidizing conditions, as indicated in Supplementary Materials Table S2, is believed to be generally representative of the system, whereas results in Figure 3 indicate that nitrate-reducing conditions may exist in some locations of the aquifer. An alternative is that the observed isotopic fractionation occurred prior to the water entering the hydrologic system (such as within the WWTP). This seems somewhat unlikely given that denitrification is indicated in all three zones of the aquifer (Figure 2c) meaning that water from similarly denitrified sources would have had to migrate throughout the aquifer and enter via recharge zones that are not near one another.

3.3. Rare Earth Elements

Rare earth elements are useful indicators of redox processes [84], groundwater mixing [85], water-rock reaction [86], and for tracing anthropogenic water sources [87]. Groups of elements amongst the REEs are also useful for tracing geochemical processes, specifically gadolinium (Gd) may be indicative of anthropogenic sources because it is commonly used in magnetic resonance imaging [7,88,89], and cerium (Ce) and europium (Eu) are redox sensitive and more likely to be derived from rock weathering [86,90]. In this study, three REE anomalies are used: Gd anomaly calculated according to the method of Reference [89]; Ce anomaly and Eu anomaly, both calculated according to the method of Reference [90]. The calculation of each anomaly has a specific purpose, with the Gd anomaly being indicative of the influence of WWTP effluent [87,89] and anomalies of Ce and Eu being indicative of redox potential and geochemical reaction between water and the aquifer matrix, respectively [86,90].

The Gd anomaly and visual inspection of REE concentrations normalized to the North American Shale Composite (NASC) [91] (Figure 4) indicate that surface water is affected by WWTP effluent. The NASC is a reasonable proxy for use in this instance because the Pierre Shale forms the basal surface of the alluvial aquifer. Gadolinium anomalies (unitless) in surface water range from 0.99 to 124 and are evident in the sharp pattern of increased NASC-normalized Gd concentrations relative to nearby elements in Figure 4, specifically in groundwater wells proximal to Fountain Creek (Figure 4a) and in Fountain Creek (Figure 4c). Note however that the calculated Gd anomaly values are not displayed on Figure 4, only the visual pattern arising from Gd enrichment. For comparison, Gd anomalies of greater than 1.5 are typically indicative of WWTP inputs [87], and Reference [89] noted Gd anomalies up to nearly 140 in WWTP effluent. The range of Gd anomalies in groundwater and surface water in this study indicates that there are a range of WWTP effects, from non-affected (Gd anomaly less than 1.5 in a variety of groundwater and surface-water locations) to highly-affected (Gd anomaly up to 124 in Canals 4 and 46 in Fountain Creek). Groundwater wells in proximity to creeks also display a subdued Gd anomaly (Figure 4a), whereas groundwater in the tributary alluvium and central groundwater generally do not (Figure 4b,d). The Gd anomaly in wells proximal to the creeks indicates that streamflow losses likely contribute to recharge in these areas, whereas distal groundwater from surface-water bodies is not influenced by streamflow losses from the creek or WWTP effluent.

The Ce anomalies range from 0.02 to 9.9 and are substantially greater in surface water than groundwater (Supplementary Materials Figure S3). Lesser Ce anomalies in groundwater may be caused by precipitation of dissolved Ce upon transformation of Ce³⁺ to Ce⁴⁺ [86,90], which is in agreement with the generally oxidizing character of groundwater in the study area. Greater Ce anomalies in surface water are somewhat unexpected given the assumed oxic nature of creeks in the area (which would presumably scavenge Ce from the aqueous phase via adsorption and decrease Ce anomaly). Detailed research on REEs in surface waters has indicated that a complex set of in-stream processes may control Ce anomalies [84], and the dataset in this study does not allow for mechanistic understanding of these in-stream processes such as adsorption, mineral precipitation, or speciation. The Eu anomalies have a smaller range from 0.007 to 2.11 and display the opposite pattern from Ce anomalies in that groundwater wells have greater Eu anomalies than surface-water sites (Supplementary Materials Figure S3). Reference [90] illustrated that Eu anomalies were poorly correlated with redox potential, and Reference [86] found that Eu anomalies were useful as an indicator of geochemical reactions within the aquifer matrix. Given these previous research findings, the pattern of Eu anomalies in the Fountain Creek aquifer is interpreted as indicative of weathering reactions within the aquifer, whereas the Eu anomaly of surface water is likely imprinted from anthropogenic sources (as corroborated by Gd anomaly). The imprint of human activity on Eu anomaly is also evident given that NASC normalized Gd concentrations are used in calculating Eu anomaly, meaning these Gd and Eu anomalies are not entirely independent.

Although general REE patterns and anomalies of Gd, Ce, and Eu are helpful for evaluating the interconnectivity of groundwater and surface water and the presence of WWTP effluent [87], the visual relations in NASC normalized diagrams do not allow for quantitative mixing or other processes to be evaluated. The REE normalization plot [85,90,92], however, may be useful for such analyses (Figure 4e). This plot helps to further contextualize the results graphically displayed in NASC normalized diagrams. For this analysis, the following classifications of REEs were used; light rare earth elements (LREEs) = La, Ce, Pr, and Nd; middle rare earth elements (MREEs) = Sm, Eu, Gd, Tb, and Dy; and heavy rare earth elements (HREEs) = Ho, Er, Tm, Yb, and Lu. See References [85,90,92] for further information on the construction of such plots.

Results shown in the normalized REE plot indicate that surface water of Fountain Creek and Canal 4 with pronounced Gd anomalies tend to plot near the bottom center, whereas other surface water does not plot in a tight grouping. Groundwater plots in a wide range of compositions with tributary alluvium groundwater tending to plot closer to the MREE = LREE line and other groundwater plotting farther to the right. To illustrate the utility of normalized REE plots in understanding mixing of water sources and geochemical reactions, a synthetic mixing line was calculated between two locations with near end-member compositions of REEs, T02-MW006 (tributary alluvium groundwater) and Canal 4 at Headgate (surface water with apparent anthropogenic effects). These water-quality sampling locations were chosen based on their position in space on the normalized REE plot in Figure 4e, which is a distinct mixing analysis from that presented earlier for Cl/Br ratios (Supplementary Materials Figure S2). Points along the synthetic mixing line are labelled with the fraction of groundwater in the mixture (x). The mixing line is not linear but takes an arcuate path, and mixing fractions of equal interval are not plotted at equal intervals along the mixing line.

The purpose of the synthetic mixing line is to show that central groundwater falls outside of the zone that may be explained by simple conservative mixing between water types with distinctive REE chemistry (REE end members). This is somewhat unexpected considering central groundwater wells are located downgradient from the tributary alluvium sites and Canal 4 (Figure 1). There are two possible explanations for this: (1) A third unsampled end-member water composition exists that contributes to the chemistry of the central groundwater, or (2) water–rock interaction produces a shift in the REE composition. The second option appears likely given that the compositions of central groundwater are becoming more enriched in terms of REE concentrations compared to the NASC, which was noted by Reference [85] to indicate geochemical evolution. This analysis provides evidence that groundwater in the Fountain Creek alluvial aquifer is evolving geochemically along the flow path from recharge to discharge zones and that zones of focused groundwater recharge exist. Specifically, based on Gd anomaly in U-7 (Figure 4a), it appears that Fountain Creek loses flow to groundwater recharge in the northern part of the study area in the reach downstream from the gage at Fountain Creek below Janitell Road below Colorado Springs, Colorado (USGS site number 07105530), and that the tributary alluvium was recharged by flow lost from discharge of the WWTP on the East Fork Sand Creek (based on Gd anomaly of 1.75 in MW1-1). This REE mixing analysis is complementary to the mixing analysis of Cl/Br presented previously in that both analyses indicate that groundwater in the tributary alluvium and proximal to Fountain Creek have end-member compositions, whereas central groundwater has a composition that is a mixture between the two end members with the addition of external solutes (both from anthropogenic and natural sources).

3.4. Pharmaceutical Compounds and Wastewater-Indicator Compounds

Previous studies on CECs have commonly used detection frequency as opposed to absolute concentrations [4], and that approach is used here, along with spatial analysis of concentrations. The pharmaceutical and wastewater-indicator compounds most commonly detected from this study are listed in

Table 1. Summary of percent detections of contaminants of emerging concern in water-quality dataset from groundwater and surface water associated with the Fountain Creek alluvial aquifer and occurrence of compounds from similar studies within the United States.

Compound Name	Percent Detections in Groundwater Samples	Percent Detections in Surface-Water Samples	Occurrence in Groundwater Downgradient from Wastewater Effluent in Minnesota, USA 1	Occurrence in Surface Water near Denver, Colorado, USA 2
Methylbenzotriazole	70%	92%	X
Carbamazepine	57%	83%	X	X
Sulfamethoxazole	13%	83%	X	X
Fluconazole	40%	83%	X
Lidocaine	20%	83%	X	X
Hexamethylenetetramine 3	7%	92%
Bupropion	13%	75%	X	X 4
Gabapentin	7%	83%		X
Tramadol	10%	75%	X	X
Meprobamate	7%	83%	X	X
Metformin	3%	33%	X	X

¹ Reference [4]. ² Reference [2]. ³ One blank sample had a reported concentration above the laboratory reporting limit (LRL). The concentration in the blank sample was 5% of the average environmental concentration, indicating little possible bias from blank samples for this compound. ⁴ Specifically noted presence of metabolite of this compound, not the compound itself.

. Typical uses and sources of many of these compounds are summarized in References [2,93,94].

Several compounds show differing abundance in groundwater versus surface water whereas others are found in similar distributions of samples. Compounds that were more prominent in surface water than groundwater are meprobamate, gabapentin, tramadol, metformin, bupropion, sulfamethoxazole, hexamethylenetetramine, and lidocaine. Other compounds in Table 1 generally had similar incidences of detection between groundwater and surface water. The overall high incidence of detection of these compounds in both groundwater and surface water show that the Fountain Creek alluvial aquifer is affected by WWTP effluent, consistent with REE analysis.

Many of the most commonly detected compounds in the Fountain Creek alluvial aquifer have been noted to occur in other settings [2,4]. An exception is hexamethylenetetramine, which was not noted to occur in other studies. One blank sample had a value greater than the LRL for hexamethylenetetramine of 36 ng/L. Compared to the average environmental concentration of 720 ng/L, this small potential bias of 5% is considered unlikely to affect results and interpretation. Hexamethylenetetramine may be derived from agricultural operations and can be associated with some pesticides [95] or from a variety of pharmaceutical uses [96]. Hexamethylenetetramine is believed to be rapidly leached from soils to groundwater and to be both abiotically and biotically degraded in the environment [95]. The greater prevalence of hexamethylenetetramine in surface water compared to groundwater (Table 1) supports the assumption that rapid degradation of the compound is occurring in groundwater.

Maps of CEC compounds in groundwater and surface water (Supplementary Materials Figure S6) support the interconnected nature of the system. Discharge points of CECs into surface waters are the WWTPs shown as potential solute sources on Figure 1, both adjacent to Fountain Creek in the northern portion of the study reach and on the east fork of Sand Creek near groundwater monitoring well T02-MW006. Canal 4 is a conveyance for water partially sourced from WWTPs and is believed to lose flow to the groundwater system based on field observations. The spatial distribution of pharmaceutical compounds indicates that solutes may be added to the aquifer from seepage from Canal 4 (Supplementary Materials Figure S6). Based on the spatial distribution of CEC concentrations (Supplementary Materials Figure S6), different CECs may have different sources and controlling processes in the hydrologic system. For example, methylbenzotriazole and carbamazepine are relatively widely distributed in groundwater and are found nearly ubiquitously in surface waters, whereas other compounds, such as bupropion and gabapentin, are sparingly observed in groundwater but are abundant in surface water. The variable occurrences are likely linked to different sources, transport properties, and degradation of some compounds in the environment [2].

3.5. Principal Component Analysis

Several iterations of PCA were completed [9], first on groups of parameters (Supplementary Materials Figures S7–S9), then on a subset determined to be most useful for discriminating between water sources and geochemical reactions (Figure 5). The final subset of parameters included the following: Hexamethylenetetramine, methylbenzotriazole, carbamazepine, and fluconazole were included because of their high frequencies of detection and contrasting behavior between groundwater and surface water (see Section 3.4); U and Se were included because they are known geologically derived constituents in Fountain Creek [17]; Gd and Eu were included because Gd is associated with WWTP effluent and streams, and groundwater displays spatially variable Eu-anomaly (see Section 3.3 and Supplementary Materials Figure S3); Cd and Mn were included because of previous investigations showing point sources in the study area [22] and trace-element-specific PCA (Supplementary Materials Figure S8); Rb was included because of trace-element-specific PCA (Supplementary Materials Figure S8) and because it may be indicative of weathering reactions in the aquifer matrix [86]; and finally NO₃ + NO₂ was included because it is primarily of anthropogenic (WWTP) origin and because of the potential influence of denitrification in the aquifer.

Results of PCA are shown in Figure 5. Rays (black lines) correspond to each parameter, and the lengths of the rays are proportional to the parameter variability within the two plotted components PC1 and PC2; the angles between any two rays is a measure of the correlation between the two variables [50]. The position of the symbols relating to each water-quality sample is a function of the score of that sample against the two PCs. Rays for constituents with a clear anthropogenic source tend to plot in the lower left quadrant with rays of nearly equal length. Constituents of known geologic origin U and Se [17] plot nearly on top of one another and are located near the NO₃ + NO₂ ray on the right side of the plot. Manganese, Cd, and Rb plot in the upper quadrant whereas Eu plots close to the CEC rays in the lower left, but with a shorter-length ray.

The location of samples in component space supports previous analyses in the identification of distinct end members and mixing. Most surface waters (especially those affected by WWTP effluent) plot between the Cd ray and the CEC rays. Three of the five groundwater wells in proximity to the creek also plot in this region. Reference [22] noted the presence of some trace metals including Cd and Mn in stormwater stream samples, and this analysis indicates that these metals along with pharmaceutical and wastewater- indicator compounds are a tracer of surface waters and are indicative of surface-water infiltration where they are found in groundwater samples. Tributary alluvium groundwater all plot on the right half of the plot in the arc between the Rb and fluconazole rays, with the exception of T01-MW002. This well is located downgradient from the former WWTP on Sand Creek (Figure 1) and appears affected by denitrification and nitrate reduction (Figure 3) such that it plots opposite from the NO₃ + NO₂ ray (because the well is devoid of NO₃ + NO₂ and thus this parameter has little control on the composition of the well) and closer to the rays representative of WWTP effluent. Aside from T01-MW002, PCA for all other tributary alluvium groundwater indicates correlation with geologically derived constituents Rb, Se, and U. Central groundwater tend to plot in an intermediate zone between tributary alluvium groundwater and groundwater proximal to the creek and surface water. This is consistent both with mixing of water sources (as indicated previously by Cl/Br and REE analyses) and with the occurrence of anthropogenically and geologically derived constituents in the central zone of the aquifer.

3.6. Noble Gases and Groundwater Residence-Time Tracers

3.6.1. Noble-Gas Isotopes

Isotopes of noble gases can be informative as to residence times and mixing of distinct waters [12,13,35,97]. Helium isotopes (Figure 6a) indicate that groundwater proximal to the creek and in the tributary alluvium are generally close to equilibrium with air-saturated water (ASW) for the study area (calculated using the method of Reference [98] and assuming 14 °C and 1700 meters elevation NAVD88) based on ratios near 1 of ³He/⁴He in the sample to ³He/⁴He in air (R/R_a). Several samples in the tributary alluvium have R/R_a values greater than 1, indicating ³He production. Given the abundance of tritium in all groundwater samples (discussed below), the R/R_a values greater than 1 are likely attributable to tritiogenic ³He (³He_trit; [54]). Samples in the central groundwater zone display lesser R/R_a values (0.75 to 0.90; Figure 6a) and greater ratios of ⁴He/²⁰Ne. These patterns are consistent with production of ⁴He in the subsurface relative to atmospherically derived ²⁰Ne likely caused by increasing residence time and interaction with aquifer materials [13]. A notable exception to the pattern of lesser R/R_a values in central groundwater wells is NONUM-2, which although centrally located in the aquifer, is also downgradient from Johnson Reservoir (Figure 1). Noble-gas isotopes indicate that this water is not geochemically evolved and may be sourced from the reservoir.

Possible recharge from the mountain block to the west is indicated by the calculated noble-gas recharge temperature (NGRT) of 4.3 °C for well TH-5, substantially colder than other wells in the aquifer (Supplementary Materials Table S4). Well TH-5 is located near Fountain Creek and is aligned with a paleochannel sourced from higher elevations along the mountain front on the western part of the study area (Figure 1). The presence of abnormally cold NGRTs has been used in the past to recognize mountain-block recharge [62,63,97]. In this instance, the alluvial aquifer may have a recharge source to the west, which would have gone undetected without the application of noble-gas modeling.

Calculation of apparent ³H/³He ages and application of lumped-parameter models are subject to uncertainty regarding the ratio of ³He/⁴He in water that is no longer in equilibrium with the atmosphere [12]. To evaluate subsurface He production, the method of Reference [99] to derive site-specific ratios for the addition of terrigenic He components was applied (Figure 6b). These results corroborate those of ⁴He/²⁰Ne in indicating geochemical evolution in the aquifer. A full sensitivity analysis focused on subsurface He inputs is included in Supplementary Materials Section S5 of the (Supplementary Materials Table S3).

Interaction with aquifer materials is also illustrated by Xe isotopes (¹³²Xe/²⁰Ne; Figure 6c). Several wells near Fountain Creek (specifically U-7, A1, and TH-5) have ¹³²Xe concentrations above that explainable by equilibrium with ASW, consistent with either interaction of water with the underlying Pierre Shale, or inflow of groundwater sourced from the shale [13]. These wells are near a reach of the creek known to have variable groundwater–surface-water interactions [20], and based on the Xe enrichments in groundwater wells in this area, it is considered plausible that there is some groundwater circulation within the Pierre Shale at the base of the alluvial aquifer. Upward leakage from the shale (which is typically considered a confining layer) could be attributable to sustained groundwater pumping in the area, which is known to cause leakage from confining layers [100].

Figure 6. Noble-gas isotope compositions, namely (a) R/R_a versus ⁴He/²⁰Ne (R_a = 1.4 × 10⁻⁶), (b) helium evolution diagram (e.g., see Reference [99]), and (c) ¹³²Xe/²⁰Ne versus ⁴He/²⁰Ne. The helium evolution diagram shown in (b) illustrates the ratio of ³He/⁴He in samples corrected for air-bubble entrainment (R_c) versus ⁴He derived from solubility (⁴He_eq) over total corrected ⁴He (⁴He_c). Wells of interest are labeled with site alias, for discussion. The composition of air-saturated water (ASW) was calculated, using the method of Reference [98] and the atmospheric composition from Reference [101], solubility at 14 degrees Celsius (°C), and 1700 m above North American Vertical Datum of 1988 (NAVD88). Extents of analytical uncertainties are within the points on all plots.

Figure 6. Noble-gas isotope compositions, namely (a) R/R_a versus ⁴He/²⁰Ne (R_a = 1.4 × 10⁻⁶), (b) helium evolution diagram (e.g., see Reference [99]), and (c) ¹³²Xe/²⁰Ne versus ⁴He/²⁰Ne. The helium evolution diagram shown in (b) illustrates the ratio of ³He/⁴He in samples corrected for air-bubble entrainment (R_c) versus ⁴He derived from solubility (⁴He_eq) over total corrected ⁴He (⁴He_c). Wells of interest are labeled with site alias, for discussion. The composition of air-saturated water (ASW) was calculated, using the method of Reference [98] and the atmospheric composition from Reference [101], solubility at 14 degrees Celsius (°C), and 1700 m above North American Vertical Datum of 1988 (NAVD88). Extents of analytical uncertainties are within the points on all plots.

3.6.2. Apparent Ages and Presence of Atmospheric Tracer Gases

Apparent ³H/³He ages (e.g., see Reference [67]) were calculated for all wells in the dataset with applicable data and assuming He isotope-abundances derived from terrigenic interaction (³He/⁴He = 2 × 10⁻⁸). In groundwater with appreciable terrigenic He (as indicated by noble-gas isotopes), the ³H/³He apparent-age calculations can be biased [102]. Moreover, some of the wells are likely subject to denitrification (Figure 3), which could cause degassing and stripping of noble gases, ³H, ³He, and SF₆ in unequal quantities proportional to their solubility [103], which would also bias ages in a variety of ways. Therefore, these results should be viewed with discretion and compared with complementary datasets as discussed below. A summary of the results of all applicable methods of estimating groundwater-residence times is presented in

Table 2. Results of apparent-age calculations, concentrations of atmospheric tracer gases, and idealized piston-flow ages of groundwater wells in the Fountain Creek alluvial aquifer.

USGS Site ID	Site Alias	R/Ra	3H (TU)	3H-3He Age (Years) 1	Fpost-1953	SF6 Concentration (fmol/kg)	SF6 Concentration (pptv) 1	SF6 Piston-Flow Age (Years) 2,3
384824104405101	03-002	1.011	6.79	3.3	0.8	2.56	8.03	5.4
384710104431201	04-009	0.942	3.11	0	0.36	2.7	7.26	7.9
384233104425801	A1	0.948	6.07	0	0.72	35.4	110.4	MF
384718104463701	3DAA	--	8.95	--	1.06	28.08	90.79	MF
384509104435901	CO259-25	0.746	4.91	10.8	0.58	10.41	22.41	MF
384719104444701	CO259-26	0.998	5.8	1.2	0.68	2.74	9.13	1.9
384437104422601	DBA	--	5.8	--	0.68	--	--	--
384112104421301	FTN1	0.995	5.91	1.1	0.7	3.82	11.9	MF
384639104461401	BAC1	--	8.49	--	1.01	83.07	196.38	MF
384949104424501	MW1-1	1.001	6.18	4.6	0.73	2.71	8.13	4.9
384848104413901	MW2-4	1.021	4.21	4.9	0.49	2.28	6.84	9.9
384408104424701	NONUM-2	1.015	4.95	2.7	0.58	4.68	8.65	2.9
384617104455901	CAD1	0.844	19.9	3.2	2.38	369.71	714.58	MF
384929104431101	T01-MW002	--	8.33	--	0.99	3.54	7.88	5.9
384956104422801	T02-MW006	--	4.58	--	0.54	--	--	--
384917104422701	T04-MW004	--	4.94	--	0.58	2.31	5.58	14.4
384818104415701	T07-MW004	--	4.22	--	0.49	3.48	7.65	6.9
384758104422301	T07-MW006	--	4.31	--	0.5	2.97	6.81	9.9
384732104430901	T13-MW004	0.982	4.62	3.2	0.54	2.57	6.33	11.4
384648104454501	TH-22	0.895	5.32	1.1	0.63	3.69	12.25	MF
384540104453601	TH-46	--	4.79	--	0.56	219.2	616.75	MF
384503104451601	TH-5	0.998	5.27	0.6	0.62	25.28	62.35	MF
384636104465401	TH-52	0.97	6.48	0.7	0.77	13.57	51.94	MF
384534104450302	U-11	--	4.52	--	0.53	3.58	8.61	3.4
384513104445302	U-12	0.808	4.96	5.6	0.58	33.84	72.45	MF
384433104440701	U-14B	0.891	4.57	6.9	0.54	98.94	248.12	MF
384420104432601	U-15	0.968	4.95	0.6	0.58	8.93	31.76	MF
384652104465101	U-7	0.968	8.98	0	1.07	7.78	20.17	MF
384604104451502	U-9	--	6.66	--	0.79	3.54	8.46	3.9
384217104402901	JCC	0.802	5.16	21.5	0.61	2.66	6.05	12.4

-- Indicates no applicable data. USGS, US Geological Survey. TU, tritium units. ¹ Calculation used components that were corrected for excess air and recharge temperature. ² “MF” indicates a failure of the piston-flow model to explain the data because SF₆ concentration is greater than that explainable by the piston-flow model. ³ When available noble gases were used to compute excess air and recharge temperature, which was used in calculating piston-flow age.

. Additional details related to environmental-tracer modeling are summarized in Reference [104].

All wells contain measurable ³H, indicating that some recharge has occurred since the 1950s and the advent of atmospheric nuclear testing [54]. All but one well have ³H concentrations in the range of 4–10 tritium units (TU). The CAD1 well, however, has a substantially greater value of nearly 20 TU, indicating that this well may have recharged during the bomb pulse of atmospheric nuclear testing in the 1960s [72]. The apparent ³H/³He age for this well is only 3.2 years, thus indicating the complexity introduced into these calculations by uncertainty in terrigenic He-isotope systematics [99,102]. These isotope systematics may be influenced by upward leakage from the Pierre Shale, which could impart a unique He signature to waters in the vicinity of high-volume pumping wells. Apparent ³H/³He ages for wells with applicable data range from about 0 to 21 years, although a majority of wells have ³H/³He ages of less than 10 years. These results have less precision from wells with noble-gas isotopes indicative of longer residence times (Figure 6) because of uncertainties arising from input of terrigenic He with an unknown ³He/⁴He ratio. The JCC well has the oldest apparent ³H/³He age of 21.5 years and a noble-gas composition indicative of interaction with the aquifer matrix (Figure 6a). Although it can be stated with some certainty that this well is in contact with older water in the aquifer, the exact age is uncertain, and true advective age for a groundwater sample does not exist because groundwater samples are mixtures of water molecules arising from converging flow paths [67].

Although all wells clearly contain a component of modern water based on the presence of ³H, this component has the ability to mask older water fractions by producing a young bias on the apparent ages [65,67]. One way of evaluating this age bias is through the calculation of the fraction of water recharged since 1953 [69]. Mixing-fraction calculations of this type indicate a range of recharge after 1953 (F_post-1953) of 0.36 to greater than 1 (Table 2). Theoretically, the F_post-1953 values should have a maximum at 1, but values greater than 1 may be obtained depending on the times series of ³H in precipitation. This study used the function from Reference [70], as detailed in Supplementary Materials Section S4. The well with the lowest F_post-1953 (and therefore the highest fraction of water recharged prior to 1953) is 04-009, which is located at the head of a paleochannel in the central portion of the study area (Figure 1) and has a Cl/Br composition representative of meteoric precipitation (Supplementary Materials Figure S2). The high proportion of older water in this well combined with the juvenile composition may indicate that well 04-009 is representative of pre-anthropogenic hydrologic conditions for the study area.

Concentrations of SF₆ (both in terms of molality and in equilibrium with recharge conditions in the study area in parts per trillion by volume, pptv) vary widely. Idealized SF₆ piston-flow ages range from about 2 to 15 years (Table 2). Comparison of these idealized ages with ³H/³He apparent ages shows that the idealized SF₆ piston-flow age tends to be greater, except in the JCC well. Several wells contain concentrations that are several orders of magnitude greater than that explainable by simple atmospheric equilibrium with recharge parameters (temperature, excess air) derived from noble-gas modeling. This phenomenon has been noted to occur both in specific aquifer lithologies [105,106] and in spatial association with urban areas [107]. Given that sediments in the alluvial aquifer are partially derived from weathering of granitic rock, it is considered likely that local contamination occurred from terrigenic SF₆ due to U-series decay [108]. However, local variations are extreme, and additional analysis is warranted to separate terrigenic from potential other sources.

3.6.3. Lumped Parameter Models

Results of lumped parameter models (LPMs) can indicate the mean residence time of groundwater [65,67]. The piston-flow model (PFM), exponential-mixing model (EMM), and dispersion model (DM) were applied in this study. Selected results according to different LPMs are illustrated in Figure 7. Additional details related to environmental-tracer modeling are summarized in Reference [104].

Comparison of tracer-tracer plots indicate variable agreement with simplified LPMs. A number of wells in the central groundwater zone and several in the tributary alluvium are consistent with both the EMM or the DM and have a range of mean residence times from 5 to 10 years (CO259-25, CO259-26, NONUM-2, T13-MW004, TH-22, U-12, U-14B). These mean residence times are generally in agreement with ³H/³He apparent ages and idealized SF₆ piston-flow ages (Table 2). Data plot closer to the LPM lines when considering ³He_trit versus ³H (Figure 7a) as opposed to ³He_trit versus SF₆ (Figure 7b) or SF₆ versus ³H (Figure 7c). This is because of likely terrigenic and external SF₆ sources, which have caused excess SF₆ to be present in the system.

Several wells are consistent with a PFM, namely A1, MW2-4, and TH-5, with mean residence times of generally less than five years. Two of these wells are located proximal to Fountain Creek (A1, TH-5) and have Gd anomaly consistent with infiltration of surface water (Figure 4). The short mean residence times, agreement with the PFM, and presence of surface-water REE signatures indicate that groundwater in the immediate vicinity of Fountain Creek receive a majority of their recharge from creek losses and that the flow paths from the creek to the well are simple enough such that they approximate piston flow.

Several wells cannot be explained by any of the applied LPMs, specifically CAD1 and U-7. Further analysis illustrates several possible reasons for these wells not fitting LPMs. For instance, U-7 has ³H contents consistent with recent recharge, no ³He_trit, and excess SF₆. The water from this well also has one of the largest Gd anomalies of any groundwater well (Figure 4). A possible explanation for the lack of agreement with any LPM is that this well was very recently recharged by water lost from Fountain Creek (which is in agreement with the REE chemistry and the ³H/³He apparent age of zero years, Table 2). The excess SF₆ is more enigmatic, but one possible explanation is that a more complex model for the formation of excess air beneath the losing stream is required (i.e., aside from the CE, PR, and UA models considered) [35]. Reference [109] found that reasonable results of dissolved-gas modeling near areas of focused stream losses were not always properly linked to physical processes in the aquifer, and that more detailed analyses may be required under these conditions. Considering CAD1 (which has non-zero ³He_trit and bomb pulse ³H), the inconsistency may be caused by differential diffusion within the aquifer. Diffusion of ³He_trit is more rapid than ³H [103], and thus it is likely that some ³He_trit has been lost via diffusion during groundwater flow, offsetting the point in Figure 7a for CAD1 down and to the right, off the DM line. An alternative explanation for lack of fit for U-7 and CAD1 to LPMs is that binary mixing between old and young waters of substantially different ages is occurring [11]. Although this may be possible, it cannot be explored with the dataset given than no tracers of “old” groundwater (e.g., ¹⁴C) were collected as part of this study. Taken as a whole, the LPM results affirm the findings of previous research that non-uniqueness and uncertainty may be limited by applying multiple tracers [56,67,110].

A comparison of the three methods of estimating groundwater age (idealized SF₆ piston-flow age, ³H/³He apparent age, and mean residence time from LPMs) is shown in Figure 8, which allows age bias within the dataset to be explored and groundwater-residence times to be visualized. Based on Figure 8, ³H/³He apparent age may be biased low, as a majority of wells have apparent ages less than five years old. This may be caused by somewhat complicated He mass balance in the aquifer (see Figure 6), introducing uncertainty in the amount of ³He_trit used in the apparent-age calculations. Idealized SF₆ piston-flow ages tend to be greater than ³H/³He apparent age but span a smaller range. Mean residence time estimates from TracerLPM have a greater number of possible analyses (because TracerLPM may produce multiple age estimates for a single sample based on combinations of different tracers), and these residence-time estimates span a much greater range than either other method (yielding a maximum mean residence time of 62 years). It is important to consider that the range of mean residence-time estimates from TracerLPM is not a measure of the full age distribution of each sample (which could include waters much older than illustrated); rather, the estimates are the distribution of possible mean ages for the entire dataset based on the tracers used in this study.

3.7. Conceptual Hydrologic Model

Combining the insights gained by the array of tracers used in this study allows for a general conceptual model of groundwater flow and solute sourcing and transport (Figure 9) that can serve to guide future studies and modeling. The influence of anthropogenic solutes is clear throughout surface waters in the study area based on a Gd anomaly and the occurrence of pharmaceuticals and wastewater-indicator compounds. Active groundwater recharge occurs both in the tributary alluvium and near Fountain Creek, although recharge in these two locations is derived from different sources. Focused recharge from streamflow loss occurs near Fountain Creek and is recognized by a Gd anomaly in both surface water and groundwater, the distribution of pharmaceuticals and wastewater-indicator compounds, and the agreement of near-creek groundwater with the PFM. Contrastingly, distributed aerial recharge in the tributary alluvium is recognized by Cl/Br ratios, groundwater ages, and noble gases nearly in equilibrium with the atmosphere. Distributed aerial recharge also fits the PFM because the sampled wells are close to the recharge zone.

Between the areas of active recharge, there is a zone of groundwater mixing and flow path convergence. This central zone of the aquifer is recognized based on agreement with the DM, older groundwater ages, and noble gases that are out of equilibrium with the atmosphere. Mixing is also supported by PCA, showing intermediate compositions of central groundwater and Cl/Br ratios. Longer residence times in this zone of the aquifer promote enhanced water-rock reaction and biogeochemical reactions as evidenced by REE relations and stable isotopes of nitrogen.

Dissolved solutes display spatial patterns that are controlled by the presence of anthropogenic point sources, as well as groundwater-residence times. Concentrations of a subset of major ions and trace elements (Cl, Cr, and Pb), a subset of REEs (Ce and Gd), and wastewater-indicator compounds are related primarily to proximity to surface-water bodies that lose flow to groundwater. Other constituents are derived from water-rock reactions and weathering within the aquifer matrix, such as HREE, U, Se, and ⁴He. This latter group of constituents tends to be enriched in the central portion of the aquifer where groundwater-residence times are longest.

There are several implications of these relations for future studies focused on hydrologic characterization in complex settings. First, a comprehensive dataset may allow for the full range of hydrologic conditions and geochemical processes to be recognized (e.g., see Reference [14]), but investigators may wish to tailor analytes if specific processes are key to the investigation. For example, without the use of REEs in this study, it may have been possible to recognize recharge from streamflow loss (by pharmaceutical and wastewater-compound indicators), but recognition of geochemical evolution in the central part of the aquifer would have been less evident. Moreover, the use of REEs is much more straightforward for recognizing the influence of WWTP effluent throughout the hydrologic system as the Gd anomaly alone is an indicator of effluent migration whereas use of pharmaceuticals and wastewater-indicator compounds is costlier analytically and may be complicated by complex quality-assurance aspects [16]. A second implication is that tracers of groundwater-residence time (be they quantitative or qualitative) are highly useful. Quantitative groundwater-residence time estimates allow for direct calculation of recharge rate [111] and are indicators of transience in the hydrologic system [110,112]. More qualitative estimates of residence times can be used to discriminate different flow-path configurations (i.e., piston flow, dispersion, and convergence) and may indicate the presence of old groundwater via noble-gas isotopes. Noble gases are also unique in their applicability for investigating recharge conditions such as groundwater recharge temperature [63], which was useful in this study for indicating a possible groundwater recharge source to the west.

Ideally, this conceptual model would be used to evaluate the distribution of specific compounds to attribute point sources of these compounds. Potential compounds of interest in this aquifer include pharmaceutical and wastewater-indicator compounds and PFAS [5]. This study includes data on pharmaceutical and wastewater-indicator compounds but does not directly include data for PFAS. Instead potential sources of PFAS are evaluated based on groundwater ages and records on industrial activity.

The most easily discriminated potential solute source for pharmaceutical and wastewater-indicator compounds is the WWTPs throughout the area, given that these compounds occur ubiquitously in surface waters and are found in groundwater of varying ages (i.e., Table 2 and Supplementary Materials Figure S6). Sources of PFAS are more difficult to attribute based on groundwater-age data alone because of the range of calculated groundwater-residence times. All groundwater observed in this study has at least a proportion of modern recharge, but in order to attribute certain compounds to a single source area the active period of solute loading from the given source area must be known [113]. PFAS compounds in the area have at least three possible sources: WWTP effluent [22], former metal-plating operations [114], and AFFF use [5]. These solute sources generally overlap in time (from the 1970s to 1990s for metal-plating operations [114], from the 1970s through the present for AFFF use [27], and for WWTP discharge [28]). Because most groundwater is less than 10 years old, it is plausible that PFAS compounds could be derived from either WWTP or AFFF, though metal-plating operations are less likely. To more fully investigate PFAS sourcing, it would be necessary to evaluate PFAS distributions and compositions, using statistical techniques (e.g., see Reference [10]).

4. Conclusions

An extensive geochemical dataset was collected, to investigate the hydrology and geochemistry of the Fountain Creek alluvial aquifer. This dataset included major and trace elements, stable isotopes, REEs, pharmaceutical and wastewater-indicator compounds, noble gases, and groundwater-age tracers. Each of these constituents has shown use in tracing geochemical reactions or resolving interactions between groundwater and surface water [7,11,14,63].

This study supports previous studies in indicating that groundwater recharge in the Fountain Creek alluvial aquifer is influenced by losses from surface water [18]. Groundwater can be separated into three distinct geochemical zones, which are spatially associated with preponderance of aerial recharge, proximity to surface waters, or groundwater mixing. The tributary alluvium, located in upgradient recharge areas of the aquifer, has a juvenile geochemical composition reflective of recent equilibrium with the atmosphere, lack of interaction with the aquifer materials, and generally young groundwater. The central groundwater zone displays geochemical evolution from the recharge zone and contains groundwater with the longest residence time. Groundwater bodies near Fountain Creek, the primary surface-water feature in the area, were also recently recharged, but they have the geochemical character of creek water that is influenced by WWTP effluent, and thus are partially derived from streamflow losses. The use of REEs is essential to understanding groundwater–surface-water interactions affected by WWTP effluent.

Most groundwater in the hydrological system appears to have been recharged within the past 20 years. However, noble-gas isotopes, lumped-parameter modeling, and the fraction of water recharged since 1953 (F_post-1953) are consistent with a proportion of older groundwater in the system. The overall groundwater ages are likely biased young by the presence of recent recharge [65,67], and the presence of older groundwater may have gone undetected without the application of noble gases and their isotopes. More detailed evaluation of mixing of old and young waters would require an intermediate-age tracer, such as carbon-14, with the ability to define an old end member based on solute concentrations and detailed modeling (e.g., see Reference [11]).