The environmental and human toxicity of PFAS has led to numerous health studies that associate PFOA with testicular, kidney, prostate and ovarian cancer, and developmental diseases [Lopez-Espinosa et al., 2011; Vaughn et al., 2013; Vieira et al., 2013; Interstate Technology and Regulatory Council (ITRC), 2020]. PFOA enters the body primarily through oral or inhalation pathways and exists in serum albumin; it also may concentrate in the liver and kidney and be excreted in urine and feces (Kudo and Kawashima, 2003). PFOA was projected to be present in the blood, kidneys, and liver of every living person (Kannan et al., 2004; Emmett et al., 2006). Concentrations in factory workers, measured as PFO⁻ anion in blood, were found to be as high as 100 parts per million (ppm) (Ubel et al., 1980; Olsen et al., 1998; Kaiser et al., 2010). Residence time in the body can range from 3.5 to 8.5 years (Olsen et al., 2007). In 2006, the U.S. Environmental Protection Agency (USEPA)classified PFOA as a likely carcinogen.

PFOA and related PFAS are used as surfactants during emulsion polymerization in the production of fluoropolymers, notably Teflon™ (PTFE) [Environmental Protection Agency (EPA), 2016]. PFOA was first synthesized in 1947 by the 3M Corporation™ using electrochemical fluorination (ECF), a method used to create many perfluoroalkyl carboxylates, and added to PTFE production in 1951 [Environmental Protection Agency (EPA), 2016]. DuPont created Teflon™ in 1938, but only in 1951 did DuPont begin buying PFOA from 3M for use in fluoropolymer production at the Washington Works plant in Parkersburg, West Virginia (PFAS Project Lab, 2022). The chemical stability of PFOA and its ability to lower surface tension made it extremely useful for a variety of products in the post-war industrial boom. Commercial uses included cleaning agents, metal plating, fire-fighting foams, coating formulas, inks, varnish, polyurethane production, shampoo, food packaging, floor polish, vinyl polymerization, gasoline, lubricants, oil, and water repellants for paper, leather, and textiles [Holzapfel, 1966; OECD (Organisation for Economic Co-operation and Development), 2002; Boulanger et al., 2005; Lehmler, 2005; Prevedouros et al., 2006; Interstate Technology and Regulatory Council (ITRC), 2020]. The latest and most widespread use of PFOA has been in Gore-Tex™ and Teflon™ products (Steenland, 2010).

Worldwide annual production of PFOA was estimated to be 260 metric tons in 1999 (Prevedouros et al., 2006); since 2002, production and emission rates have fallen, replaced by newer products such as C6 (a 6-carbon chain PFAS). Mechanisms of PFOA contamination of the environment include smokestack emissions during production, leftover product disposal, and factory runoff (Barton, 2008; Xiao et al., 2015). During production, hot (>300 degrees Centigrade) (Schroeder et al., 2021) PFOA vapor is emitted from smokestacks (Step A-Figure 1), contacts cool air and is carried in the lower atmosphere (Step B-Figure 1), coalesces from vapor to form submicron particles (0.1–1 μm) that adsorb to aerosol surfaces, and is deposited onto the land surface (Step C, Figure 1), typically within 5–10 km of the emission source [Interstate Technology and Regulatory Council (ITRC), 2020; Schroeder et al., 2021]. Partitioning and transport of PFOA and other PFAS in soil and groundwater are largely controlled by factors such as topography, water and soil pH, aquifer geochemical properties, humidity, temperature, wind direction and soil type (Barton, 2008; Guo et al., 2020). Although PFOA is very soluble and has a low sorption potential (Ferry et al., 2012), leaching in the unsaturated zone above the water table (Step D, Figure 1) can be retarded by organic carbon, iron oxides and hydroxides, or retention at the water-air interface (Ferry et al., 2012; Guo et al., 2020). For example, adsorption to organic and inorganic particles in the unsaturated zone can lead to an 18-to-64-fold PFAS retardation relative to water transport rates, meaning that PFAS may also continue to occur in soils and sediments decades after release and air deposition (Xiao et al., 2015). Once transported downward to the saturated zone, further migration of PFAS occurs within the surficial aquifer in groundwater flowing through more permeable materials (i.e., sand and gravel) and to the FRA (Steps E and F, Figure 1), where angular or lateral movement along planar discontinuities such as beds, foliations, fractures, and faults occurs (Kim et al., 2014, 2016, 2018, 2019).

Bedrock mapping was completed using traditional field techniques that included use of Brunton Compass, Global Positioning System (GPS) satellite receiver, topographic base map, rock hammer for sampling, hand lens for mineral identification, and field notebook for recording all data. Bedding, foliation, and fracture data were plotted on rose diagrams and equal area nets using Daisy software (Salvini et al., 1999). The field equipment for surficial mapping included a trenching shovel, hand auger, trowel, Munsell Color Chart, and topographic maps generated from LiDAR elevation data. The collation of data sets and map cartography were completed using ArcMap software.

The State of Vermont maintains a database of groundwater wells drilled since ~1970 that has information on well depth, overburden thickness, driller's yield (gallons/minute), static water level (SWL), casing depth and length, and groundwater production zones (Vermont Agency of Natural Resources, 2023). The bedrock wells used to produce the Geographic Information Systems (GIS) maps for this study were located either by GPS or by correlation with E911 street addresses. The overburden data were used to produce bedrock surface contour and isopach maps, whereas the SWL data were used to generate potentiometric surface maps. Although the SWL data in this database do not necessarily represent equilibrium water levels, these data can be judiciously used at the scale of a town (303 km²) or greater, since groundwater flow directions at large scales generally follow topography. Ground surface elevations for wells were determined from LiDAR Digital Elevation Models (DEM).

We sought a relatively even distribution of wells throughout the field area, but were limited by our ability to obtain permission from well owners, as well as the constraint that wells were absent in zones served by public water supplies. Seventeen domestic bedrock wells were logged using the following Mt. Sopris Instruments tools: (1) temperature and resistivity (2PFA-10000), (2) natural gamma (2PGA poly gamma), (3) caliper (2PCA-1000 caliper), (4) acoustic borehole imager (ABI) (QL40), and (5) heat-pulse flowmeter (HFP2293). A continuous color video record of most wells was made using a GeoVISION borehole camera video system. The geophysical data and videos from these 17 wells are currently being collated by the Vermont Geological Survey with funding from the USGS National Geological and Geophysical Data Preservation Program (NGGDPP) and will be available on the Vermont Geological Survey website during late 2023.

The temperature/fluid resistivity probe was used in tandem with the gamma probe. Temperature and fluid resistivity measurements were recorded from the top of the water level downward at a 0.05 m interval as the probes were lowered. The logging speed of temperature and resistivity values was typically between 3.0 and 4.6 m/min. Well logging software converted fluid resistivity measurements to conductivity and specific conductivity (@25°C). Natural gamma measurements were acquired every 0.05 m, as the probes were raised. The logging speed was between 3.0 and 4.6 m/min. Caliper logs were completed by raising a three-prong caliper logger from the bottom of the well to the bottom of the well casing at a sample rate of 0.05 m with a logging speed between 3.0 and 4.6 m/min. Acoustical Borehole Imaging was completed by raising the imager from bottom of the well in order to maintain tension on the cable during imaging at a sample rate of 0.004 m. The logging speed was typically about 1.8 m/min. Logging was terminated at the bottom of the casing or the top of the water, whichever occurred first.

The velocity measurements using a heat-pulse flow meter were made above and below a major fracture. When the bedrock is highly fractured, measurements were made at various intervals. Care was taken to make certain the diverters were oriented properly, pointing in the direction of flow. Multiple measurements were made at each location until measurement values were duplicated. Upward velocity measurements with a peak time >45 secs were discarded because times beyond this reflect thermal density stratification from the heat pulse and not borehole groundwater flow. We assumed those peaks were a result of density stratification of the water due to heating.

Prior to logging, the water line to the home was disconnected and the pump and associated hoses and wiring were removed by a well driller and coiled onto a plastic tarp near the well. The wells were logged under ambient conditions. After logging was completed, the pump was replaced and run until the water ran clear. The pump-pulling and logging equipment was decontaminated with a 1:3 detergent/water mix when entering and after exiting the well. The logging data for each well were collated in a standard format using WellCAD™ software in the order discussed above. The attitudes of planar structures intersecting the well were calculated from Acoustic Televiewer and caliper logs using WellCAD™ software. Structural data for the wells were plotted on rose diagrams and contoured equal area nets using Daisy Software (Salvini et al., 1999) and compared with structural data from bedrock geologic maps.

Groundwater samples were collected over a 3-year period from rural groundwater wells that are the water supplies for private homes. These wells are 15 cm in diameter and steel-cased through surficial deposits to bedrock, and then open to the FRA below the base of casing. Water was sampled using in well production pumps; water treatment systems (e.g., water softeners and/or Point of Entry Treatment (POET) systems using activated carbon) were bypassed, if present. Most samples were obtained directly from the faucet on the pressure tank or from an outdoor spigot. A calibrated YSI or Hanna multiprobe was used to measure water temperature, specific conductivity, pH, and oxidation-reduction potential (ORP), —samples were collected once all parameters had attained steady values, typically after purging for 10 to 20 minutes.

Samples were collected in acid-washed, high-density polyethylene (HDPE) or glass bottles (for P analysis) and stored on ice. Concentrations of metals, metalloids and anions were analyzed at the VTDEC Laboratory in Burlington VT. Samples were acidified and digested before analysis, and methods by which metals and anions were determined are shown below. For representative precision, accuracy, and detection limit data, see this reference (Kim et al., 2014). The following methods apply: (1) EPA Method 6020 (ICPMS, Inductively Coupled Plasma Mass Spectrometry) for Al, As, Ba, Ca, Co, Cu, Fe, Pb, Mg, Mn, Mo, Ni, K, Na, Sr, Total Calculated Hardness (TCH), and U. (2) EPA Method 245.1 for Hg; (3) EPA 300.0 or SM 4500CL-G for Cl; (4) EPA Method 300.0 for SO₄ and NO₃; (5) SM 4500-PH for P; (6) SM 4500 SiO₂-F for dissolved SiO₂; and SM 2320B for alkalinity.

Groundwater chemical composition was analyzed using traditional approaches (e.g., Piper and Durov plots, geochemical spider plots, routine statistical analysis) and by a non-metric multidimensional scaling (NMDS) analysis. Spider diagrams plot the relative abundance of components useful in discrimination (e.g., TDS, Ca, Mg, Na, K, Sr, U, SiO₂, SO₄, Cl) associated with selected Bennington groundwater samples in comparison to their abundance in a reference water (in this study, the reference is the mean value of all groundwater samples analyzed in this study). The differences in the patterns of relative abundances (rather than the abundance values themselves) are what is useful in discrimination. NMDS is a common tool used for working with large datasets. The goal of NMDS is to collapse information from multiple dimensions (e.g., from dozens of chemical parameters) into just a few (e.g., two), so that similarity and difference can be interpreted. Unlike Principal Components Analysis (PCA), which primarily relies on Euclidean distances, NMDS uses rank orders, making it a more flexible technique that can accommodate a variety of datasets. To plot NMDS results on a map, a spreadsheet was created to spatially analyze sample ID and corresponding Universal Transverse Mercator (UTM) coordinates with a single-axis NMDS analysis value (the first dimensional axis). Initially NMDS values ranged from−0.046 to 0.135. These values were normalized with the lowest value equal to 0 and highest to 1, and subsequently multiplied by 100 for ease of use in ArcGIS. The resulting columns were added to ArcGIS as a feature class. The X and Y values of UTM were input as XY data and displayed spatially. Corresponding numerical NMDS data were added as unique values for each sample point, and visually shown with a color wheel (green to red). An empirical Bayesian kriging analysis was performed to create a contoured raster file. In all, three different Kriging models were run, with the Empirical Bayesian Kriging analysis best articulating the relation between points for this use. The Kriging model and sample points are displayed with a satellite base map in the background, along with a layer containing major rivers. In this analysis, slate-influenced waters plot as red and carbonate-dominated waters plot as green.

Samples for stable water isotopic analysis were collected in 15-ml HDPE bottles with tight fitting screw caps. In the laboratory, samples were filtered at 0.22 μm and placed into 2 ml glass auto-sampler vials with a PTFE septa. The isotopic composition (²H-H₂O, ¹⁸O-H₂O) of hydrogen and oxygen of the water molecule of surface water, precipitation and ground water were measured by wavelength scanned cavity ring-down spectrometry on un-acidified samples by Picarro L-1102i WS-CRDS analyzer (Picarro, Sunnyvale, CA). Samples were vaporized at 110°C. A 5 μl Hamilton glass syringe draws 1 μl of sample to inject into a heated vaporizer port (110°C). For each injection, the absorption spectra for each isotope were determined 20 times and averaged. Between injections, the needle was rinsed with 1-Methyl-2-pyrrolidinone and the sample chamber is flushed with dry nitrogen. Three standards that isotopically bracket the sample values were run alternately with the samples. Secondary lab reference waters (from Boulder, Colorado; Tallahassee, Florida; Amherst, Massachusetts) were calibrated with Greenland Ice Sheet Precipitation, Standard Light Antarctic Precipitation and Vienna Standard Mean Ocean Water (VSMOW) from the International Atomic Energy Agency (IAEA). Results were calculated based on a rolling calibration so that each sample is determined by the three standards run closest in time to that of the sample. Long-term averages of internal laboratory standard analytical results yield an instrumental precision of 0.51 ‰ for δ²H-H₂O and 0.08 ‰ for δ¹⁸O-H₂O.

Once field parameters such as temperature, conductivity, pH, and ORP had stabilized, most groundwater samples for chlorinated fluorocarbons (CFC11, CFC12, and CFC113), sulfur hexafluoride (SF₆) and tritium analysis were acquired using sampling methods and protocols of the Tritium Lab [The Tritium Lab at the Rosenstiel School of Marine and Atmospheric Sciences (RSMAS), n.d., n.d.], which calls for collection in 500 ml glass bottles with foil-lined screw caps (CFCs, SF₆) and 1 liter plastic bottles (tritium), respectively. Most samples were analyzed at this lab. Eight samples were collected by U.S. Geological Survey (USGS) personnel using USGS methods for tritium and chlorofluorocarbons (U.S. Geological Survey Groundwater Dating Lab methods), which are very similar to those of The Tritium Lab. Age tracers for USGS samples were analyzed at the USGS Groundwater Dating Laboratory in Reston, VA except for tritium, which was analyzed at Woods Hole Oceanographic Institute, Woods Hole, MA.

We had access to PFAS data from groundwater samples taken by consultants under contract to the VTDEC from private wells (and springs) throughout the contamination area (n = 628). This dataset was last accessed for this study in May 2019. If more than one PFOA result existed for a well, the maximum value was used. These samples were procured via the methods of Weston and Sampson Engineers Inc (2016). Additional samples referenced in this study were collected from springs and wells by faculty and students at Bennington College using the same methods as VTDEC above (Schroeder et al., 2021). All samples referenced in this study were analyzed at Eurofins Eaton Analytical Labs with USEPA method 537 or 537.1.

Many wells in the study area were tested multiple times for PFOA and the temporal variability of these results for an area near the former factory were assessed (Schroeder et al., 2021). Bennington College students and faculty studied the temporal variability of PFOA contamination by selecting three bedrock wells close to ChemFab with different initial PFOA concentrations (Schroeder et al., 2021). The three wells selected lie within 250 m of one another and are all within 330 m of the former ChemFab facility. They had initial PFOA concentrations of 7.1 ppt, 82 ppt, and 2,090 ppt on March 1, 2016. Each well was sampled between seven or eight additional times between March 2016 and July 2017 (Schroeder et al., 2021).

The concentration of PFOA varied over the 1.5-year sampling period, but the total variation of each well did not cross the ranges of the other two wells (Schroeder et al., 2021). These results indicate that spatial variation is much greater than temporal variation for PFOA concentrations in the FRA. The high spatial variation is likely due to long-term, structural variation in flow patterns rather than short-term turbulence in the system (Schroeder et al., 2021).

By correlating well locations with well driller logs (data on well depth, overburden thickness, well yield, static water levels, and casing depth) in the study area, bedrock surface contour, isopach, and potentiometric surface maps were produced using GIS software (Kim and Dowey, 2017). All these maps integrated LiDAR-derived ground-surface elevation data with the well logs.

The primary purpose of the bedrock-surface contour (BSC) map (Kim and Dowey, 2017) for the detailed area was to determine if linear features could be correlated with bedrock structures (Figure 6). To “ground-truth” whether the linear features on the BSC map could be correlated with brittle structures (fractures) and ductile structures (thrust faults), the statistically dominant steeply-dipping (> 60 degrees) fracture sets for geographic domains on the fracture map (Kim, 2017b) were superposed on this map along with the thrust faults (Figure 6).

The Walloomsac River valley, which represents the lowest elevations in the region, has distinct N-S and E-W segments (Figure 6), as well as intersecting N-S valleys, some that are filled with anomalously thick surficial deposits. Most orthogonal segments are parallel to either the N-S or E-W trending fracture sets, which form fracture intensification domains (FIDs). The NNE-SSW striking fracture set is parallel to the valley that is bounded by both segments of the Breese Hollow Thrust. Other fractures of this attitude are parallel to axial surfaces of major fold structures and both thrust faults (Figure 6). The NW-SE striking fracture set has a more subtle effect on the BSC surface by elongating some contours in that direction.

The isopach map for the detailed area shows how overburden thickness varies throughout the detailed field area (Figure 7A). Like the BSC map, the isopach map correlates strongly with the underlying bedrock structure. For example, thick (> 6 m) surficial deposits parallel the hanging wall of the Green Mountain Frontal Thrust and fill two fracture-controlled bedrock basins on the north side of the Walloomsac River. Thick surficial deposits are also present in the eastern one-third of the field area. Thicker surficial deposits may retard the downward movement of contaminants.

The potentiometric surface contour map for the FRA in the detailed field area (Figure 7B) shows that groundwater in the FRA generally follows topography and flows toward the Walloomsac River valley; this over simplifies groundwater flow paths and will be discussed further in later sections. The data points used to create Figures 6, 7 are found in Supplementary Table 1.

A total of 17 domestic bedrock wells were logged for this study under ambient conditions after the pump and associated piping and wiring had been removed. The wells ranged in depth from 38–183 m. Composite WellCad plots were assembled from the logs for all wells and included: (1) temperature, (2) specific conductivity, (3) gamma, (4) caliper, (5) acoustic televiewer (AT), and heat-pulse flowmeter. Continuous borehole camera video was recorded for most of these wells. When available, driller's logs for neighboring domestic wells that were not logged with geophysical equipment were consulted. The purpose of this section is to summarize the logs of the wells that were most informative to understanding the hydrogeology of the contamination area.

Figure 8 is a summary of wells logged from west to east that shows well depth, elevation relative to sea level, overburden thickness, depth of casing, bedrock formation, water-producing structures, flow directions, and PFOA concentrations. The wells are grouped according to the structural position they occupy on the bedrock geologic map and cross sections of Figure 5. The W-E groupings of wells are based on whether they are completed: (1) In gray phyllites of the New York Slice west of the Breese Hollow Thrust's surface trace, (2) Through limestones in the hanging wall of the Breese Hollow Thrust into the New York Slice (foot wall), (3) In limestones of the western North Bennington Slice, (4) In dolostones of the eastern North Bennington Slice, (5) Through phyllitic quartzites in the hanging wall of the Green Mountain Frontal Thrust into dolostones and limestones of the Bennington Slice (foot wall), and (6) In dolomitic quartzites of the Bennington Slice, where thick ocher zones are present.

Vertical flow in the wells was measured using the Heat Pulse Flowmeter under ambient conditions. Under ambient conditions, groundwater exchange between water in the borehole and individual fractures intersecting the borehole is a function of the transmissivity of the individual fracture and difference in hydraulic head within the fracture and the water elevation in the borehole (Paillet, 1998). Flow within the open borehole has little effect on the surrounding aquifer (Paillet, 1998) or mimics the regional groundwater flow.

Wells #3-4 and #14-15 were drilled through the Breese Hollow and Green Mt. Frontal thrusts, respectively, where upward wellbore flow directly from the fault or fractures below it produced upward wellbore flow, and static water levels tended to be near or flowing at the surface (flowing artesian). Faults can be recognized at depth through the abrupt lithologic changes observed on gamma, acoustic televiewer, and borehole camera logs. In Well #14 a steeply-dipping, E-W striking fracture ~10 m below the Green Mt. Frontal Thrust (GMFT), presumably related to the Walloomsac River valley fracture zone, produced upward wellbore flow. In Well #15, a 1.2 m cavern at a depth of 40 m produced upward wellbore flow. Wells #14-15 had anomalously thick overburden because they were completed in a filled N-S trending bedrock basin that intersects both the GMFT and the Walloomsac River valley fracture zone. Well #6 penetrated the eastern segment of the Breese Hollow Thrust (BHT-E) at 12.2 m and weak borehole groundwater flow was both upward and downward from this structure.

Wells completed in the internal parts of the New York, North Bennington, and Bennington slices produce groundwater from bedrock structures (primarily bedding) at various depths in the aquifer and the flow directions in these wells were upward, downward, or bidirectional. Well #13 is the deepest well in this group and produced upward-flowing water from a steeply-dipping, E-W striking fracture at 172 m, which is part of the Walloomsac River valley fracture zone.

Two of the wells logged internal to the Bennington Slice have 61 m+ thick ocher horizons, and both produce groundwater from the bottom of the well. The third well has a 5 m weathered bedrock horizon near the surface that may be ocher and produced groundwater from 141 m.

During this investigation, the VTDEC action level for PFOA in drinking water was 20 ppt (PFOA, n.d.). Of 628 wells (and springs) tested, 486 (77.4%) had detectable PFOA and 355 (56.5%) had PFOA concentrations > 20 ppt. The PFOA data were accessed in May 2019 from a VTDEC database (Vermont Dept. of Environmental Conservation, n.d). A reconnaissance-level examination of the complex fate and transport of PFOA in the contamination area is possible through the following sequence: (1) contour the PFOA results in groundwater and look for spatial patterns, (2) examine the spatial distribution of PFOA that occurred from air deposition from the ChemFab factory [air deposition envelope (ADE)], and (3) compare results between #1 and #2 (Figure 10) to see if the calculated PFOA air deposition patterns correlate with those in groundwater. Since > 80% of the wells in the Town of Bennington are bedrock wells, we are primarily concerned with the fractured rock aquifer (FRA) (Vermont Agency of Natural Resources, 2023).

The ADE is the approximated land area over which airborne wind-transported PFOA was deposited on the land surface and leached to groundwater in detectable concentrations. The envelope was delineated by extensive sampling of soil and shallow groundwater seeps/springs (Schroeder et al., 2021).

We compared the groundwater contours of PFOA with ADE and observed the following:

Although the “Air Deposition Envelope” is not a perfect representation of the PFOA distribution in the contamination area, observations 2–4, which are marked on the caption of Figure 10 suggest that heterogeneous hydrogeological processes in the surficial and fractured rock aquifers likely play a role in the fate and transport of PFOA, beyond the air transport and deposition geometry.

Waters analyzed in this study (Supplementary Table 2) are predominantly calcium-bicarbonate trending toward intermediate. Low total dissolved solids (TDS), low-NaCl groundwater from the anomalously low PFOA zone is distinguishable from surrounding groundwaters that contain higher TDS and NaCl (and PFOA). The graphical analyses that were most useful in discriminating water types are shown in Figure 11, including a bivariate plot of [alkalinity/(K+SiO₂)] vs. [Ca/SO₄], a spider diagram that shows relative differences of water types relative to the average of all wells in the study area (N = 59), and a ternary diagram with axes of [alkalinity/(K+SiO₂)], [Ca/SO₄], and [Sr/Mg].

Groundwater hydrochemical signatures that correspond to spatial location and bedrock type are distinguished in Figure 11 as follows:

All groundwater sample δ²H-H₂O and δ¹⁸O-H₂O values plot close to the regional meteoric water line (Cole and Boutt, 2021; Vermont Dept. of Environmental Conservation, n.d) with varying isotopic compositions from−12 to−8 δ¹⁸O-H₂O ‰ (Figure 12). The stable isotopic data can be found in Supplementary Table 3. See Figure 9C for all sample locations. The isotopic composition of the waters is distinguished by the varying amounts of heavy isotopes present in the waters. The Bennington College (BC) till well samples have the most enriched in heavy isotopic composition (more positive) with an internal sample range of ~1.5 δ¹⁸O-H₂O ‰. The GG1 and GG2 groups vary by a similar amount (1.5 δ¹⁸O-H₂O ‰) but are distinctly more depleted in the heavy isotopic composition. The GG4 samples have the smallest internal sample variability and cluster between the upper and lower end of the GG1/2 and BC till wells composition. The GG5 samples have the largest range (> 2 δ¹⁸O-H₂O ‰) and also cover the whole range of non-till groundwater samples. They do not cluster but rather are distributed across this range. Compared to this distribution, the GG3 wells tightly cluster between−11 and−10.5 δ¹⁸O-H₂O ‰ and are significantly more depleted than most other samples. The geochemical groups were defined based on elemental chemical components (not isotopic compositions), nevertheless there is some relation between the geochemical groups and the variability and isotopic composition of the waters.

All water samples (except for GG5) were collected within a relatively narrow surface elevation band above sea level (ASL) (171–287 m) and show no statistically significant trend in composition with elevation. The till wells are the most enriched, show seasonal variability, and because they are water table samples, they likely range in direct vertical recharge conditions across the lower portions of the river valley. The most depleted isotopic compositions of the till samples are interpreted to be reflective of average cold-season recharge conditions along the valley bottoms (lower elevations). The overlap in the GG4 group and the most depleted till well samples suggest that the GG4 wells receive local, low-elevation recharge waters. GG1 and GG2 samples range from this enriched end member to depleted values and are interpreted to be mixtures of low-elevation and perhaps high-elevation recharge waters. The GG5 group waters range in elevation from 317–792 m ASL with a median value of 658 m ASL, with most of these waters being spring samples. Thus, these waters likely reflect high-elevation recharge waters with a mixture of cold and warm season recharge conditions. The GG3 samples are from wells ranging from 171–226 m ASL and plot in a distinct cluster apart from other low elevation well samples (e.g., GG1, GG2, GG44). These waters are isotopically distinct from the rest, with the exception of high elevation spring water samples of GG5. Thus, these GG3 waters likely reflect a source of recharge in the high-elevation to the east or are waters that were recharged to the FRA locally, but during a period with more depleted (i.e., colder) low-elevation recharge water. The relation between precipitation water isotopic composition and elevation was not statistically valid in the region (Cole and Boutt, 2021). This is likely due to the absence of higher elevation isotope precipitation samples. A local elevation-dependent isotopic analysis does not exist for the region (Cole and Boutt, 2021). The stable H₂O isotope samples reflect source-water characteristics compared to the elemental geochemical information which is influenced by weathering and residence time characteristics.

CFC11, CFC12, CFC113, SF₆, and tritium concentrations were analyzed from 25 groundwater samples in the field area to determine recharge ages (Figure 9D). Recharge-ages were not analyzed for all major and trace element geochemistry samples. All recharge-age data are presented in Supplementary Table 4 whereas USGS data alone is in a data release (Shanley et al., 2023). Average/Apparent recharge ages were calculated from the CFC concentrations using appropriate recharge temperatures and elevations and converted to year of recharge (Plummer and Busenberg, 2000). Figure 13A shows the recharge years with standard deviations for the samples from each geochemical group. Tritium data in tritium units (TU) are also plotted for the samples in each geochemical group in Figure 13B. The CFC and tritium data in Figures 13A, B are raw/unprocessed.

Figure 13A generally shows “older” and “younger” CFC age groups, with the former being GG1 and GG3 and the latter GG 2 and GG4. Samples 12–14 in GG 3 and one samples in GG 1 plot below the approximate current tritium threshold for “modern” water (~0.2–0.4 TU, (Plummer and Busenberg, 2000), indicating recharge that predates the presence of tritium in rainfall from atmospheric atomic bomb testing in 1953 (Figure 13B). CFC data for these same samples indicate recharge around 1953. Pre-modern waters are older than 1953 and modern waters are younger than 1953.

Samples 3–4 also have CFC recharge ages near 1953, however, they also have much higher levels of Tritium (~4 TU), which suggests a modern source (s) of Tritium, possibly related, in part, to the 1963 Tritium peak in atmospheric precipitation (Figures 13A, B) (Solomon and Cook, 2000). Sample 1 has a much younger CFC age (1969) and a Tritium level of ~6 TU, indicating a significant component of modern water. Sample 15 has a CFC recharge-age of 1973, but a Tritium concentration of ~1 TU, indicating a significant proportion of modern water.

With one exception, all GG2 and GG4 samples have younger CFC recharge-ages that range from ~1980–1992, which fall within the timeframe of ChemFab plant operation (Barr Engineering, 2018) (Figure 13A). Tritium abundances of ~4–7 TU also indicate modern water. Sample 8 in GG2 has an older CFC recharge age of 1962 and Tritium level of ~1 TU that suggests a mixture of pre-modern and modern waters (Figure 13B).