Isotopic Discrimination of Aquifer Recharge Sources, Subsystem Connectivity and Flow Patterns in the South Fork Palouse River Basin, Idaho and Washington, USA

Abstract

Groundwater studies in the South Fork Palouse River Basin have been unable to determine recharge sources, subsystem connectivity and flow patterns due to the discontinuity of pathways in the heterogeneous and anisotropic aquifers located in Columbia River flood basalts and interbedded sediments. Major ion, δ¹⁸O, δ²H, δ¹³C, δ³⁴S and temperature for groundwater collected from 28 wells of varying depths indicate a primary recharge source dominated by snowmelt along the eastern basin margin. This recharge can be separated into two distinct sources—a deeper and relatively less altered snowmelt signal (−17.3‰ to −16.8‰ δ¹⁸O, −131‰ to −127‰ δ²H, −12.9‰ to −10‰ δ¹³C, 18–23 °C) and a more altered signal likely derived from a shallower mixture of snowmelt, precipitation and surface water (−16.1‰ to −15.5‰ δ¹⁸O, −121‰ to −117‰ δ²H, −15.9‰ to −12.9‰ δ¹³C, 12–19 °C). A mixing of the shallow and deep source waters is observed within the upper aquifer of the Grande Ronde Formation near Moscow, Idaho, which results in a homogenization of isotope ratios and geochemistry for groundwater at nearly any depth to the west of this mixing zone. This homogenized signal is prevalent in a likely primary productive zone of an intermediate depth in the overall aquifer system.

1. Introduction

Since 1935, groundwater levels have declined in aquifers of the Palouse River Basin, particularly in the South Fork Palouse River Basin (Basin) in north-central Idaho and eastern Washington (Figure 1) [1,2,3]. The Basin aquifers are contained in the fractured basalts of the Columbia River Basalt Group (CRBG) and interbedded sediments of the Latah Formation that compose much of the Basin [4]. Groundwater provides a primary source for industry and agriculture in the region [5] and is the sole source of municipal water in the Basin [2,6]. Extrapolation of the current trend in declining groundwater levels in the Basin indicates the possibility of insufficient groundwater resources to meet future community needs [6,7].

Multiple aquifers are contained in the sediments of the Latah Formation and the flood basalts of the Wanapum and Grande Ronde formations of the CRBG (Figure 2). The Latah Formation is composed of siliciclastic sediments varying from very fine clays and silts to cobbles and boulders [8]. This formation is thickest along the eastern margin near the mountain-front interface and becomes thin interbeds between successive basalt flows to the west [9]. The fractured basalts of the Wanapum Formation (Figure 2) and sediments of the Latah Formation contain the shallowest of the Basin aquifers [9,10,11]. The underlying Grande Ronde Formation is composed of dense, fractured basalt flows and interbedded sediments [4,9,12]. The Grande Ronde Formation lies unconformably atop granitic and metamorphic basement rocks, which also compose the surrounding elevated features such as the Palouse Range (Figure 1) [12]. The Grande Ronde Formation has been hypothesized to contain potentially two aquifers (upper and lower) that may be connected or independent depending on location in the Basin [13,14,15].

Previous studies have suggested that groundwater recharge is entering the Wanapum and Grande Ronde aquifers through pathways beginning in sediments along the eastern margin of the Basin [7,8,16]. Due to low soil permeability and limited connectivity between basalt flows [8,9,16], it is unlikely for diffuse recharge to occur with infiltration of precipitation. The combination of variable permeability, basalt fracture termination and discontinuity of basalt flows and interbedded sediments produces a heterogeneous and anisotropic aquifer(s) matrix. Past attempts to model the aquifers to predict the decline in groundwater levels have produced mixed results due to lack of resolution in recharge sources, subsystem connectivity and identifiable flow patterns [17]. Questions remain concerning surface-water recharge locations and connections between aquifers and within the aquifers, such as flow between the central and western portions (Figure 1 and Figure 2) of the Basin [6,11,15].

Municipal pumping in the Basin produces a widely varying potentiometric surface in the Wanapum and Grande Ronde aquifers as well as substantial well interference [14]. Adjacent wells screened at similar depths within CRBG aquifers may show orders of magnitude variation in hydraulic properties [18]. Past investigations have been able to characterize small portions of the Wanapum or Grande Ronde aquifers but these methods have failed to quantify characteristics that can be applied to individual aquifers or the entire aquifer system [19,20,21]. To overcome limitations of hydraulic tests, isotope and geochemical properties of groundwater from wells in and near Pullman and Moscow were analyzed to identify water sources, subsystem connectivity and flow patterns within this portion of the Basin. The goal of the analysis was a discrimination of source waters at each location from any upgradient surface water or groundwater that was contributing to the aquifer(s) at the sample location. Resolution of groundwater recharge will assist with future modeling attempts for preservation of groundwater resources in the Basin. Discrimination of source waters, subsystems and flow patterns will have utility for further investigation of aquifers within the entire CRBG province.

2. Materials and Methods

2.1. Wells, Sample Collection and Field Parameters

In 2017, groundwater was collected from 28 private and municipal wells following sufficient purging of each well casing (minimum of 3 well volumes). Wells were selected for sampling based on location in the Basin, presence of a dedicated pump and regular pumping of the well—seven wells had screens in the Latah Formation, six wells had screens in the Wanapum Formation and nineteen wells had screens in the Grande Ronde Formation (

Table 1. Descriptions of wells east of Moscow (Eastern) and wells in and near Moscow (Central), Idaho, sampled as part of this study (Figure 1). Subgroup identifiers relate to likely source waters and depth. Well identifiers are based on formation and depth of screen interval(s) as determined from the geologic layer model (W = Wanapum, G = Grande Ronde, UG = Upper Grande Ronde, DG = Deep Grande Ronde, L = Latah Formation) and location (E = Eastern, C = Central).

Study Group	Subgroup Identifier	Well Identifier	Location (Lat./Long., NAD83)	Well Elev. (m NAVD88)	Well Depth (m)	Screened Depth (m)
Eastern	Shallow	WE1	46.7241, −116.9428	794	82	20–44; 69–82
		WE2	46.7256, −116.9555	797	150	32–83
Central	Surface Connection	UGC1	46.7374, −117.0619	771	93	72–93
		LC1	46.7297, −117.0278	796	106	86–90; 95–101
	Shallow	LC2	46.7193, −117.0376	838	100	96–100
		WC1	46.7349, −117.0025	783	174	12–73
		WC2	46.7351, −117.0023	783	173	21–72
		UGC2	46.7351, −117.0249	778	224	209–228
		WC3	46.7533, −117.0646	796	78	66–78
	Mixing Zone	LC3	46.7435, −116.9724	797	105	73–105
		GC1	46.7346, −117.0324	779	382	193–205; 215–222; 239–244; 369–375
	Deep	DGC1	46.7404, −117.0132	798	444	320–444
		DGC2	46.7410, −116.9954	788	399	334–398
		DGC3	46.7370, −117.0210	782	407	201–236; 297–335; 366–407

;

Table 2. Descriptions of wells in and near Pullman (Western), Washington, sampled as part of this study (Figure 1). Subgroup identifiers relate to likely source waters and depth. Well identifiers are based on formation and depth of screen interval(s) as determined from the geologic layer model (G = Grande Ronde, UG = Upper Grande Ronde) and location (W = Western).

Study Group	Subgroup Identifier	Well Identifier	Location (Lat./Long., NAD83)	Well Elev. (m NAVD88)	Well Depth (m)	Screened Depth (m)
Western	Surface Connection	GW1	46.7357, −117.1763	715	70	12–70
		GW2	46.7134, −117.1824	746	216	205–216
		GW3	46.7475, −117.1730	739	158	72–158
		UGW1	46.7301, −117.1711	725	84	20–31; 63–84
		UGW2	46.6969, −117.1500	753	122	9–122
	Mixing Zone/Deep	GW4	46.7589, −117.1673	741	108	90–108
		GW5	46.7260, −117.1137	771	92	79–91
		GW6	46.7358, −117.1765	715	219	83–107
		GW7	46.7228, −117.1765	766	243	100–243
		GW8	46.6925, −117.2416	780	183	122–183
		GW9	46.6925, −117.2414	780	183	131–183
		GW10	46.7342, −117.1568	773	214	121–214
		GW11	46.7291, −117.1698	736	678	164–176; 203–215; 225–234; 274–286; 307–553
		GW12	46.7320, −117.1497	794	247	165–212; 229–244

). A Hanna HI9829 multi-parameter probe with in-line flow cell was used to measure pH (calibrated to pH of 4 and 7), oxidation-reduction potential (ORP, calibrated to a +200 mV vs. Ag/AgCl standard), temperature and conductivity (calibrated to a 1000 μS/cm at 25 °C standard) to ensure sufficient purging and stabilization of field measurements prior to sample collection.

2.2. Groundwater Sample Analyses

Unfiltered and filtered water was collected from each well for analysis of physicochemical properties, major and minor element concentrations, δ²H and δ¹⁸O of water (water isotopes), δ¹³C of dissolved inorganic carbon and δ³⁴S of dissolved sulfur (

Table 3. Geochemical and isotope analyses of groundwater samples collected from 28 wells in the South Fork Palouse River Basin.

Analyte	Collection	Filter	Preservation	Analysis Method
Cations (mg/L)	Acid-washed HDPE, 125 mL	0.45 µm	1 mL HNO3	ICP-MS
Anions (mg/L)	Nalgene, 125 mL	0.45 µm	Chilled	IC
Alkalinity	Nalgene, 125 mL	0.45 µm	Chilled	Titration: inflection pt.
δ13C (VPDB, ‰)	Glass/polyseal, 240 mL	0.20 µm	Chilled	Gas-ratio MS
δ34S (CDT, ‰)	3 L carboy	0.45 µm	Chilled	Continuous Flow IRMS
δ2H, δ18O (VSMOW, ‰)	Glass/polycone, 60 mL	Unfiltered	Chilled	Cavity ring-down spec.

). Alkalinity of filtered (0.45 µm) samples were analyzed by inflection point titration using a Hach digital titrator, Thermo Scientific Orion Star pH meter, stir plate and 0.16 N H₂SO₄ digital titrator cartridge. Filtered (0.45 µm) samples were analyzed for metal concentrations (Ba, Ca, Cd, Cr, Co, Cu, Fe, K, Mg, Mn, Mo, Na, Ni, V and Zn) and anion concentrations (SO₄²⁻, Br⁻, Cl⁻, F⁻, PO₄³⁻, NO₂⁻ and NO₃²⁻) at the University of Idaho’s Analytical Sciences Laboratory using inductively coupled plasma optical emission spectrometry (PerkinElmer Optima 8300) and ion chromatography (Thermo Scientific Dionex Aquion), respectively. Water isotope values of unfiltered samples were determined using a Picarro L2130-i Analyzer (±0.025‰ for δ¹⁸O and ±0.1‰ for δ²H) located at Boise State University’s Stable Isotope Laboratory. Values of δ¹³C were determined from dissolved (0.2-µm filtered) inorganic carbon at the University of Arizona’s Environmental Isotope Laboratory using an automated KIEL-III carbonate preparation device coupled to a Finnigan MAT 252 gas-ratio mass spectrometer (±0.08‰). The 0.2-µm filtering for δ¹³C determination was chosen to reduce the need for introduction of an anti-microbial agent for sample preservation [22]. Total sulfur in 0.45-µm filtered samples were converted to SO₄, precipitated as BaSO₄ and δ³⁴S was measured as SO₂ on continuous flow isotope ratio mass spectrometer at the University of Arizona’s Environmental Isotope Laboratory. Isotope values are presented as either individual values (groundwater samples) or averages (precipitation, snow, groups of groundwater samples) plus or minus (±) the standard error of the mean. Quality control and accuracy were checked with instrument blanks, replicate samples and calibration standards over the course of sample collection and analysis. Analyses of metal and anion concentrations were evaluated for acceptable (±5%) charge balance error. No samples were outside of the acceptance criteria. All groundwater data are included as supplementary data, and available online at https://www.mdpi.com/2306-5338/6/1/15/s1.

2.3. Precipitation Isotope Data

Basin precipitation was collected by personnel of the U.S. Network of Isotopes in Precipitation (USNIP) at the WA24 site, Palouse Conservation Farm, in Pullman, Washington (46.7606, −117.1847 NAD83; Figure 1) at an elevation of 766 m (NAVD88). Precipitation samples were collected using USNIP standard collection equipment, composited as monthly samples and analyzed for δ¹⁸O and δ²H using traditional isotope-ratio mass spectrometry and(or) laser absorption spectroscopy (±0.1‰ for δ¹⁸O and ±0.8‰ for δ²H) [23,24,25]. Precipitation δ¹⁸O values for the higher elevation Palouse Range bordering the eastern margin of the Basin (Figure 1) were modeled using elevation fractionation rates [26]. An elevation range of 1000–1500 m (NAVD88) was used to model likely δ¹⁸O values for the extent of the Palouse Range from its foothills to crest using the values from the USNIP site. The elevation of the precipitation monitoring station (766 m NAVD88) was subtracted from each 100 m interval of the selected Palouse Range elevation and the difference multiplied by the fractionation coefficient of 0.28 for a gradation of δ¹⁸O values above the USNIP site. Modeled elevation values of δ²H were produced using the local meteoric water line (δ²H = 7.23*δ¹⁸O – 5.56, R² = 0.95) based on the USNIP data.

2.4. Snow Isotope Data

Snow samples were collected in March 2018 near the peak of the Palouse Range (Figure 1) for analysis of water isotopes. Snow samples were collected in 1-L vacuum bags from a 1-m deep vertical profile excavated in the snowpack adjacent to a U.S. Natural Resources Conservation Service Snow Telemetry (SNOTEL) site (Station #989, 46.8062, −116.8381 NAD83, 1,350 m NAVD88). Samples were collected at 0.15-m intervals from the surface of the snowpack to the bottom of the profile. While frozen, each sample was vacuum sealed and allowed to melt while isolated from the surrounding environment. Each sample was transferred into 120-mL glass containers and submitted to the Stable Isotope Laboratory for analysis.

2.5. Model of Geologic Layers

A 3D geologic model of the subsurface of the eastern Basin was created with Golden Software’s Surfer (v. 15.5.382, Golden, CO, USA) and Voxler (v. 4.3.771, Golden, CO, USA) to evaluate well and screen locations in relation to geologic layers. Geospatial data were collected and compiled from existing well logs and each geologic layer was based on Bush et al.’s [9] interpretation of the Latah, Wanapum and Grande Ronde formations. Surfer was utilized for the interpretation of grid layers (kriging surface based on a linear variogram) for the top and bottom of each basalt flow. Voxler was used to create a 3D representation of each layer by rendering volumes between each surface using a trilinear method of interpolation between gridded data.

2.6. Principal Component Analysis

Principal components analysis, PCA (XLSTAT-Base software, v. 2018.5, New York, NY, USA) was used as a screening tool to identify linear combinations of data that may explain groundwater-composition variation across the study site. All concentration data were converted to mM/L for analysis. The dimensionality reduction was performed as a correlation matrix PCA, which are less influenced by differences in units of measurement [27]. Identification of correlated analytes assisted in identifying potential source waters and mixtures of source waters. Based on the percent of contribution by each analyte in describing the first component of the PCA, as well as the correlation between variables within the first component, the primary analytes of stable isotopes, water temperature and major ions (water type) were selected for further evaluation of source water discrimination, subsystem connection and flow patterns.

2.7. Inverse Modeling and Source Water Mixing

Isotope values of Basin precipitation, snowmelt and groundwater were the primary data for identifying source waters in the study area. Both surface water and groundwater were considered possible source waters contributing to a groundwater sample because of questions regarding connections and flow patterns within and between the eastern, central and western portions of the Basin. With the likely movement of groundwater from east to west, surface water sources were discriminated for the eastern Basin and resulting mixtures in the central Basin were used as source waters for groundwater in the western portion of the Basin. Mixing of possible water sources was considered for groundwater with sample values aligned between potential source waters, and source water contributions were determined through two-component inverse modeling (Equation (1)). Each groundwater sample was considered a mixture (mix = δ₃) of potential upgradient sources (e.g., snowmelt or groundwater) that were represented by the water isotope signals of the likely source waters discriminated by their water type, temperature and isotope values.

The inverse modeling of each groundwater sample was conducted with a perspective of potential groundwater evolution or step-wise progression of groundwater mixtures of different isotope signals from east to west. This inverse modeling was conducted in this step-wise progression because of a lack of confidence in upgradient-to-downgradient connections between portions of the fractured basalts and interbedded sediments. Using Equation (1), the inverse calculation allows the unmixing of the groundwater isotope signal by varying the possible fractions (f₁ + f₂ = 1) of the likely source waters given their average δ¹⁸O or δ²H values (δ₁ and δ₂). Microsoft Excel (Solver Tool) was used to perform the inverse calculation (precision of fraction contribution equal to 0.00001, convergence tolerance set at 0.0001). The inverse calculation is a best-fit scenario where the fractions of likely inputs (source waters) are varied concurrently to minimize the residual of a model solution compared to the actual value (e.g., δ¹⁸O). The convergence tolerance is the fit parameter that must be met for an output of fractional contributions to be estimated by each inverse calculation or the model was deemed unacceptable.

f₁δ₁ + f₂δ₂ = δ₃

(1)

Following initial discrimination of likely source waters and inverse modeling using the water isotope values, groundwater samples were evaluated by their δ¹³C, δ³⁴S and physicochemical properties to confirm potential source waters and evaluate subsystems and flow patterns within and between the aquifers. The carbon (δ¹³C) and sulfur (δ³⁴S) isotope values were evaluated similarly to the water isotope values but with the mass (m) of the solute included in the two-component inverse modeling equation (Equation (2)). Average values for mass and isotope ratios were created for each source water signal similar to averages for observed water isotope values of selected sources. Acceptable model results were produced for several of the groundwater samples using δ¹³C and alkalinity values, but few model results were acceptable with δ³⁴S and SO₄ values.

f₁m₁δ₁ + f₂m₂δ₂ = m₃δ₃

(2)

3. Results

3.1. Hydrochemistry Variation

Major ion concentrations of the collected groundwater indicate the largest variation in water types (hydrochemical facies; Figure 3) in and around Moscow (Eastern and Central wells). Present in this area are water types of Ca-HCO₃ or mixed cation-HCO₃ except for deeper water (DGC1, DGC2) that indicate a Na-HCO₃ water type (Figure 3). Groundwater from a majority of shallow Eastern/Central wells (UGC1, LC1, WE1, WE2, WC1, WC2) contained the largest concentrations of SO₄ (13–83 mg/L) and relatively colder water temperatures (12–14 °C). Additionally, a few of these shallow wells (UGC1, LC1) indicate detectable NO₃ concentrations that was not present in other wells. The relatively higher concentrations of SO₄ and presence of NO₃ in groundwater for certain wells in the Central and Eastern area likely indicates recharge from nearby surface water (surface connection) that may be influenced by agricultural areas [28,29]. Deeper water in the Central area (DGC1, DGC2, DGC3) contains greater Na (20–66 mg/L) and less SO₄ (2.9–5.6 mg/L). The deeper water likely has a longer subsurface travel path where Na from the basaltic host rock and siliciclastic sediments reduces Ca as the dominant cation [30], and SO₄ concentration likely is reduced because of changes in redox conditions (deep water ORP average of −240 mV compared to an average of −65 mV for groundwater from the shallower wells).

The water type variation with depth to groundwater for the Eastern and Central wells was not apparent for groundwater from the Western wells, which contained groundwater with a more homogeneous water type regardless of screen depth. Water in Western wells were a mixed cation (Ca and Na)-HCO₃ water type except for select wells that appear to have a nearby surface connection as observed with groundwater from select Eastern and Central wells. These select Western wells contained groundwater with greater concentrations of SO₄ (UGW1, GW4), a measurable quantity of NO₃ (0.85 mg/L for GW4), and(or) relatively colder water temperatures (13.5–14.5 °C), which indicate a likely nearby surface connection.

Results of the PCA (

Table 4. Principal components analysis results for groundwater data from wells in the South Fork Palouse River Basin. Listed parameters include only those that explained >5% of the variance of the first component.

Analyte	Correlation	First Component Contribution (%)
SO4 (mg/L)	−0.50	6%
Alkalinity (mg/L as CaCO3)	0.75	14%
Temperature (°C)	0.68	11%
δ2H (‰)	−0.85	18%
δ18O (‰)	−0.81	16%
δ13C (‰)	0.79	15%
δ34S (‰)	0.86	18%

) indicate that only SO₄, alkalinity, temperature and isotope ratios of H, O, C and S show significant (≥0.5) positive or negative correlation and contribution (>5%) to explaining the variance of the groundwater data in the study area. The conservative tracers of δ¹⁸O and δ²H and the non-conservative SO₄ were negatively correlated to alkalinity, temperature, δ¹³C and δ³⁴S. Each of these analytes accounted for a minor portion (<20%) of the variability of the fitted first component (initial orthogonal axis) suggestive of multiple source waters and multiple mixing scenarios as indicated by the variation in water type. From these PCA results, the water isotopes were selected for evaluation of primary water sources and possible unmixing of these sources at downgradient well locations. The remaining analytes identified by the PCA as contributing to groundwater data variation were evaluated in a similar manner to verify water isotope results and identify additional influences on water chemistry.

3.2. Precipitation and Surface Water Isotope Signals

Weekly measurements from 2006 to 2016 of δ¹⁸O and δ²H in precipitation collected by USNIP provide a view of the seasonal flux of precipitation isotopic signal in the center of the study area (valley floor at 766 m NAVD88). Fall/winter/spring precipitation/snowmelt is the dominant seasonal source of potential recharge because of limited rainfall and high evapotranspiration during summer [31]. Basin recharge is expected to occur primarily along the eastern margin where infiltration/percolation pathways can allow water to enter the Latah and CRBG formations [7,16]. The difference in elevation between the USNIP site and the Palouse Range would produce differences in the isotopic signal of potential recharge from precipitation and snowmelt [26]. After sufficient infiltration/percolation (minimization of evaporation), water isotope values can be considered conservative for tracing aquifer pathways [30].

The water isotope values for precipitation collected between October and March during 2006–2016 (winter precipitation (WP)) at the USNIP site ranged from −22‰ to −2‰ for δ¹⁸O (weighted average ± standard error of −13.2 ± 0.1‰) and −170‰ to −44‰ for δ²H (weighted average of −100 ± 1‰) (Figure 4). The crest of the Palouse Range is almost twice the elevation of the USNIP site and receives approximately twice as much annual precipitation (50 cm vs. 110 cm) [32,33]. Modeled values of the water isotopes for precipitation along the Palouse Range indicate a likely additional −2.1‰ to −0.7‰ for δ¹⁸O and −15‰ to −5‰ for δ²H with the increasing elevation from the foothills to the crest. Snow samples collected from the Palouse Range SNOTEL site indicate substantially more negative values of −19.9 ± 1.2‰ for δ¹⁸O and −146 ± 8‰ for δ²H compared to the modeled values of −15.3‰ for δ¹⁸O and −115‰ for δ²H for precipitation at the same elevation. Existing observations for δ¹⁸O values (δ²H not available) for mid-Basin creeks/rivers that convey water from the Palouse Range, such as Missouri Flat Creek (weighted average δ¹⁸O of −14.7‰ [31]) and the South Fork Palouse River (weighted average δ¹⁸O of −15.3‰ between October and March and −14.5‰ between April and September [34]), indicate δ¹⁸O values more negative than precipitation recorded at the USNIP site but more positive than the modeled precipitation and measured snow values. The isotope values of these potential source waters (Figure 4) indicate a variable evaporation effect on conveyed snowmelt/precipitation from the Palouse Range prior to recharge to the aquifers. This variation in the evaporation effect is visible in the 0.3 to 13.2‰ range of deuterium excess values for groundwater collected as part of this study. This range of deuterium excess indicates a possible divergence in the δ¹⁸O and δ²H values that may lessen the utility of the δ²H values for the inverse modeling.

3.3. Groundwater Source Discrimination

Groundwater from Central wells with an identified surface connection (UGC1 and LC1) also had δ¹⁸O values (−14.3‰ and −14.2‰, respectively) and δ²H values (−110‰ and −108‰, respectively) similar to Basin surface water (Figure 4). A nearby surface water connection agrees with past observations of surface water contribution to groundwater in wells near Paradise Creek [19,20,21] and the South Fork Palouse River near Pullman [34]. Groundwater from these wells is designated as “Surface Connection” for use as a source water to discriminate source contributions to the remaining wells.

The deepest Central wells (DGC1, DGC2, DGC3) have a mixed cation-HCO₃ water type and strongly negative δ¹⁸O (−17.3‰ to −16.8‰) and δ²H (−132‰ to −127‰) values compared to groundwater in Central wells screened at shallower depths. More negative water isotope values indicate a source water from a relatively colder or higher elevation location. This deeper signal is supported by similar δ¹³C (−12.9‰ to −10‰; Figure 5), δ³⁴S (32 ± 4‰; Figure 6), alkalinity (150 to 207 mg/L) and relatively warmer temperatures (18–23 °C) than groundwater from shallower Eastern and Central wells. These warmer temperatures can be achieved by a geothermal gradient of 2–3 °C per 100 m [35], which is within the standard gradient for the upper lithosphere. This temperature signal is attained without the water undergoing evaporation or mixing of additional surface water; thereby suggesting a primarily snowmelt isotope signal for this designated “Deep” source water, which is reflective of unaltered (or less altered) snowmelt infiltration.

Despite a few Central wells indicating a surface water source (Surface Connection) or primarily snowmelt source (Deep), groundwater from many of the shallower Eastern and Central wells appears more representative of evaporated/mixed recharge from the Palouse Range. The δ¹⁸O and δ²H values for these shallower wells with Ca-HCO₃ type groundwater (WE1, WE2, LC2, WC1, WC2 and WC3) ranged from −16.1‰ to −15.3‰ for δ¹⁸O and −121‰ to −117‰ for δ²H and had well screens ranging from 12–228 m below land surface. This limited variation in relatively negative water isotope values for Ca-HCO₃ groundwater in shallow-to-intermediate depth Central wells indicates a recharge signal more aligned with modified (evaporation effect) snowmelt or the modeled precipitation signal for the Palouse Range than the more positive Basin precipitation (USNIP), mid-Basin surface water (Surface Connection) or more negative Deep water (Figure 4). This “Shallow” source water also has similar δ¹³C values (−15.9‰ to −13.1‰; Figure 5), δ³⁴S values (9.7 ± 2.2‰; Figure 6), alkalinity values (110 to 150 mg/L) and groundwater temperatures (12–19 °C) between those of groundwater from Surface Connection and Deep wells. It is likely this group of wells represents a source water mixture influenced to varying degrees by the Deep source and the Surface Connection source. Inverse modeling of δ¹⁸O and δ²H values with end members of Surface Connection and Deep water was able to produce the observed values for groundwater in these Shallow wells (

Table 5. Inverse modeling (Equation (1)) results of δ¹⁸O and δ²H values for groundwater collected from wells in the South Fork Palouse River Basin. Two potential source waters (Surface, Shallow, Mixing Zone, Deep) were chosen according to perceived representation of groundwater from wells indicative of characteristics associated with those recharge pathways, such as detectable nitrate (Surface) or most negative isotope values (Deep). Source waters not included in the inverse model are shaded in gray.

Study Group	Subgroup Identifier	Well Identifier	δ18O				δ2H
			Surface	Shallow	Mixing Zone	Deep	Surface	Shallow	Mixing Zone	Deep
Eastern	Shallow	WE1	45%			55%	40%			60%
		WE2	40%			60%	30%			70%
Central	Shallow	LC2	55%			45%	25%			75%
		WC1	50%			50%	25%			75%
		WC2	55%			45%	35%			65%
		UGC2	55%			45%	25%			75%
		WC3	65%			35%	45%			55%
	Mixing Zone	LC3		45%		55%		25%		75%
		GC1		30%		70%		0%		100%
Western	Surface Connection	GW1	50%		50%		20%		80%
		GW2	60%		40%		10%		90%
		GW3	25%		75%		25%		75%
		GW4	10%		90%		0%		100%
		UGW1	40%	60%			0%	100%
		UGW2	20%			80%	0%			100%
	Mixing Zone/Deep	GW5			0%	100%			0%	100%
		GW6			100%	0%			0%	100%
		GW7			55%	45%			0%	100%
		GW8			50%	50%			0%	100%
		GW9			65%	35%			0%	100%
		GW10			0%	100%			20%	80%
		GW11			50%	50%			0%	100%
		GW12			0%	100%			0%	100%

). Inverse modeling of δ¹³C (

Table 6. Inverse modeling results for δ¹³C and alkalinity concentrations and evaluation of sulfate concentrations and ORP values for groundwater collected from wells in the South Fork Palouse River Basin. The same potential source waters (Surface, Shallow, Mixing Zone, Deep) used for the water isotope modeling were selected for the carbon modeling. Source waters not included in the inverse model are shaded in gray. Modeling attempts that could not produce an acceptable fit were denoted with an “―”. “…” represents no change in relative concentration (Conc.) or ORP compared to shallow recharge without agricultural influences.

Study Group	Subgroup Identifier	Well Identifier	δ13C and Alkalinity					Sulfate and Reduction
			Surface	Shallow	Mixing Zone	Deep	Conc.	Conc.	ORP
Eastern	Shallow	WE1	70%			30%	…	↑	↓
		WE2	85%			15%	…	↑	↓
Central	Shallow	LC2	40%			60%	…	…	…
		WC1	40%			60%	…	↑	…
		WC2	―			―	↑	↑	…
		UGC2	0%			100%	↑	…	↓
		WC3	30%			70%	…	…	↓
	Mixing Zone	LC3		―		―	↑	↑	…
		GC1		―		―	↑	…	↓
Western	Surface Connection	GW1	20%		80%		…	…	↓
		GW2	50%		50%		…	…	↓
		GW3	35%		65%		…	…	↓
		GW4	30%		70%		…	↑	…
		UGW1	―	―			↑	↑	↓
		UGW2	55%			45%	…	…	↓
	Mixing Zone/Deep	GW5			45%	55%	…	…	…
		GW6			0%	100%	↑	…	↓
		GW7			60%	40%	…	…	↓
		GW8			0%	100%	↑	…	↓
		GW9			―	―	…	…	↓
		GW10			0%	100%	…	…	↓
		GW11			―	―	…	…	↓
		GW12			0%	100%	↑	↑	↓

) generally supports the water isotope modeling results, with evidence of additional alkalinity in several shallow wells (WC2, UGC2) leading to a lack of acceptable model fit (Table 6).

3.4. Source Water Mixing in the Central Area

The water isotope values for groundwater from two of the Central (LC3, GC1) wells (−16.7‰ and −16.5‰ for δ¹⁸O and −131‰ and −124‰ for δ²H, respectively) and their groundwater temperatures (15–19 °C) and major ion concentrations (Figure 3) were between those of the Shallow and Deep source waters. Representative contributions produced through the water isotope inverse modeling indicate that groundwater from these wells can be recreated through mixing of Shallow and Deep groundwater (Table 5). The δ¹³C (−13.3‰ and −12.9‰) and δ³⁴S (3.6‰ and 34‰) values for these mixed groundwater samples (designated as “Mixing Zone”) also appear to be a mix of Shallow and Deep source waters, but modeling with the isotope-mass inverse calculation (Equation (2)) could not produce acceptable model results for either δ¹³C or δ³⁴S. This lack of acceptable model fit likely is a product of increased alkalinity from water-rock interaction and decreased SO₄ because of reducing conditions. For example, a relatively low SO₄ concentration and an elevated δ³⁴S isotope ratio in groundwater from GC1 indicate likely loss of SO₄; and this well contained noticeable H₂S during sampling. Additionally, there is an increase in alkalinity concentrations with deeper water as indicated by the largest alkalinity concentrations (150–207 mg/L) identified in Deep and Mixing Zone source waters.

3.5. Source Water Homogenization and Transport to the West

Groundwater from many of the Western wells (8 of 14: GW6, GW7, GW8, GW9, GW10, GW11, GW12, UGW2) contained similar (homogenized) δ¹⁸O (−17.4‰ to −16.4‰), δ²H (−132‰ to −126‰), δ¹³C (−12.8‰ to −11.9‰) and δ³⁴S (19‰ to 31.3‰) values and were a mixed cation (Ca and Na)-HCO₃ water type with similar alkalinity values (160 to 170 mg/L). Sulfate concentrations were very low for groundwater from these Western wells (<5 mg/L) and were accompanied by reducing conditions (ORP values ranging from −300 to −60 mV) and higher δ³⁴S values indicative of SO₄ reduction. Groundwater from several of these wells (UGW2, GW6, GW7, GW8, GW9) had similar water isotope values compared to the discriminated Mixing Zone source water, while other Western wells (GW5, GW10, GW11, GW12) contained groundwater more similar to the Deep source water. Inverse modeling of groundwater for most Western wells were completed using a mixture of Mixing Zone and Deep waters (Table 5). This geochemical similarity for groundwater from a majority of the Western wells and the likely Mixing Zone and Deep source waters indicate a strong connection of groundwater flow from the intermediate-to-deep Central well locations moving westward. Testing of the water isotope mixtures (Table 5) by inverse modeling of δ¹³C values also indicates a mix of Mixing Zone and Deep waters for many Western wells (Table 6: GW5, GW6, GW7, GW8, GW10, GW12). However, groundwater from several Western wells (GW1, GW2, GW3, GW4, UGW2) could not be inversely modeled because of lower alkalinity values (138 to 160 mg/L), suggesting Surface Connection water. Inverse modeling of groundwater from these wells required Surface Connection water to produce the respective isotope mixing scenarios (Table 5 and Table 6: GW1, GW2, GW3, GW4 and UGW2).

4. Discussion

Groundwater from select wells in the Eastern and Central parts of the study area indicate primary source water signals of relatively direct snowmelt (Deep source water), partially-evaporated snowmelt and precipitation likely from surface transport away from the mountain-front interface (Shallow source water), a Shallow and Deep water mixture (Mixing Zone source water) and nearby surface water (Surface Connection source water) (Figure 4, Figure 5, Figure 6 and Figure 7). These four source waters—Deep, Shallow, Mixing Zone and Surface Connection—were used to unmix all remaining groundwater samples from wells in the Eastern, Central and Western areas of the study area (Table 5 and Table 6). Deep water contains a depleted water isotope signal (Figure 4) indicative of a colder temperature or higher elevation source that undergoes relatively little evaporation or mixing with other surface waters prior to recharge. The Deep water designation is supported by similar major ion concentrations (Figure 3), δ¹³C and δ³⁴S values (Figure 5 and Figure 6) and water temperature. The Shallow water likely is a snowmelt runoff that undergoes more evaporation (time at the surface) and mixing (precipitation or surface water) than the Deep water prior to infiltration and recharge to the aquifers. The Shallow (longer surface path) and Deep (shorter surface path) waters appear to meet beneath the Central area and mix to form a new source water (Mixing Zone water) for downgradient wells. Throughout the study area, an additional surface water that has undergone even further evaporation (Surface Connection water = longest surface path) can enter the aquifer system and mix with Deep, Shallow and Mixing Zone waters (Figure 7).

Water isotope values within select shallow Eastern and Central wells indicate a dominant Surface Connection source water (Figure 5) with similar values to the modeled precipitation values for the Palouse Range (δ¹⁸O of −16.6‰ to −14.7‰). However, the observed water isotope values in groundwater from these wells are less negative than Deep water but more negative than mid-Basin precipitation (δ¹⁸O of −13.2‰) and Basin surface water (δ¹⁸O of −15.3‰ to −13.9‰). Thus, this Surface Connection source water is not explained solely by infiltrating precipitation but likely a snowmelt + precipitation + evaporation (surface water) signal that is altered with the longest transport across the Basin prior to infiltration and recharge to the aquifers. The location of Surface Connection water inputs correlates with the thinning of overlying sediments and the Wanapum Formation (Figure 2) to the west [8,9]. Surface connection input indicates infiltration of this source water to fractures or high permeability interbeds of the Wanapum and upper Grande Ronde aquifers at select locations in the Central and Western areas, which have been identified in previous studies [19,20,21,34].

The Shallow and Deep source waters in the Eastern and Central wells mix at varying ratios in the hypothesized upper Grande Ronde aquifer and form a homogenized isotope and geochemical signal (Mixing Zone) that becomes the primary groundwater extracted from wells in the Western area. Water type (Figure 4) and inverse modeling (Table 5) supports 1) recharge at the eastern margin of the Basin (Deep water) or away from the mountain-front interface (Shallow water) and 2) homogenization of isotope and geochemical characteristics with westward flow (Deep + Shallow = Mixing Zone). Discrimination of these source waters allows for inverse modeling (Table 5), or unmixing, of groundwater from Central wells, which indicate varying contributions of Deep and Shallow water to produce Mixing Zone groundwater and subsequent flow of this mixed water and additional Deep water to the Western area (Figure 7). The alteration of major ions (e.g., Ca to Na shift; Figure 4) and inability to inverse model some of the δ¹³C and alkalinity pairs and most of the δ³⁴S and SO₄ pairs (Table 6) indicate non-conservative behavior of solutes and groundwater chemistry evolution with water-rock interactions and changes in redox conditions with westward flow.

Groundwater chemistry in Western wells is similar to Mixing Zone water and(or) Deep water with some potential contributions of Surface Connection water (Figure 4 and Figure 5). Water isotope and δ¹³C values can be inversely modeled to produce groundwater from a mixture of Mixing Zone + Deep water or Mixing Zone + Deep + Surface Connection water (Table 5 and Table 6). Groundwater from Western wells also showed elevated alkalinity and low SO₄ concentrations similar to Mixing Zone and Deep water. This similar chemistry in the Western part of the study area is an indication of connectivity between the Central and Western portions of the Basin. This connection likely is the productive zone of the upper Grande Ronde aquifer that aligns with the depth of the Central Mixing Zone at about 100 to 200 m below land surface.

5. Conclusions

Evaluation of isotope and geochemical characteristics of groundwater collected from the study area indicates two major recharge sources to aquifers along the eastern margin of the Basin—a primarily snowmelt-dominant water (Deep source water) and a modified snowmelt/precipitation water (Shallow source water). These source waters are discriminated by variation in major ion chemistry, water temperature and values of δ¹⁸O, δ²H, δ¹³C and δ³⁴S. These two eastern Basin recharge sources mix near Moscow at intermediary depths and homogenize in a mixing zone (Mixing Zone source water) within the Central part of the study area. This Mixing Zone water moves westward to Pullman at a likely intermediary depth along with additional Deep water at a greater depth to produce a groundwater mixture in the Western area. Additionally, the same analytes that discriminated the Shallow, Mixing Zone and Deep source waters indicate an additional source water (Surface Connection) that is present in shallow groundwater in select areas throughout the study area. The two primary recharge sources—Shallow and Deep—can be used to unmix most groundwater found in the Central area. The Mixing Zone and Deep water can be used to unmix most groundwater found in the Western area. The utility of these isotope and geochemical relations allows for a better discrimination of recharge source waters, subsystem connections and flow patterns in the South Fork Palouse River Basin, which will allow for improved groundwater modeling to evaluate potential water resource management activities. Furthermore, the discrimination of these isotope and geochemical characteristics can be used for an evaluation of source waters in the larger Basin to incorporate potential changes in snowpack accumulations and additional municipal pumping. Given future changes in the regional climate [36], declining regional surface-water baseflows [37], uncertainty in regional recharge [38] and likely increases in groundwater withdrawals with increases in population, the spatial relations of the identified source waters likely will be altered and need to be evaluated under future conditions for protection of the resource.