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

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, δ18O, δ2H, δ13C, δ34S 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.3h to −16.8h δ18O, −131h to −127h δ2H, −12.9h to −10h δ13C, 18–23 ◦C) and a more altered signal likely derived from a shallower mixture of snowmelt, precipitation and surface water (−16.1h to −15.5h δ18O, −121h to −117h δ2H, −15.9h to −12.9h δ13C, 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.


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 (Figures 1 and 2) of the Basin [6,11,15].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 (Figures 1 and 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.

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; Table 2).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.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.

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; Table 2).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.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).

Groundwater Sample Analyses
Unfiltered and filtered water was collected from each well for analysis of physicochemical properties, major and minor element concentrations, δ 2 H and δ 18 O of water (water isotopes), δ 13 C of dissolved inorganic carbon and δ 34 S of dissolved sulfur (Table 3).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 2 SO 4 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 4 2− , Br − , Cl − , F − , PO 4 3− , NO 2 − and NO 3 2− ) 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 δ 18 O and ±0.1 for δ 2 H) located at Boise State University's Stable Isotope Laboratory.Values of δ 13 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 δ 13 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 4 , precipitated as BaSO 4 and δ 34 S was measured as SO 2 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 http://www.mdpi.com/2306-5338/6/1/15/s1.

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.1847NAD83; 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 δ 18 O and δ 2 H using traditional isotope-ratio mass spectrometry and(or) laser absorption spectroscopy (±0.1 for δ 18 O and ±0.8 for δ 2 H) [23][24][25].Precipitation δ 18 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 δ 18 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 δ 18 O values above the USNIP site.Modeled elevation values of δ 2 H were produced using the local meteoric water line (δ 2 H = 7.23*δ 18 O -5.56, R 2 = 0.95) based on the USNIP data.

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.8381NAD83, 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.

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.

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.

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 = δ 3 ) 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 1 + f 2 = 1) of the likely source waters given their average δ 18 O or δ 2 H values (δ 1 and δ 2 ).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., δ 18 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.
Following initial discrimination of likely source waters and inverse modeling using the water isotope values, groundwater samples were evaluated by their δ 13 C, δ 34 S and physicochemical properties to confirm potential source waters and evaluate subsystems and flow patterns within and between the aquifers.The carbon (δ 13 C) and sulfur (δ 34 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 δ 13 C and alkalinity values, but few model results were acceptable with δ 34 S and SO 4 values.

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 3 or mixed cation-HCO 3 except for deeper water (DGC1, DGC2) that indicate a Na-HCO 3 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 4 (13-83 mg/L) and relatively colder water temperatures (12-14 • C).Additionally, a few of these shallow wells (UGC1, LC1) indicate detectable NO 3 concentrations that was not present in other wells.The relatively higher concentrations of SO 4 and presence of NO 3 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 4 (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 4 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)-HCO3 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 SO4 (UGW1, GW4), a measurable quantity of NO3 (0.85 mg/L for GW4), and(or) relatively colder water temperatures (13.5-14.5˚C),which indicate 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 3 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 4 (UGW1, GW4), a measurable quantity of NO 3 (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) indicate that only SO 4 , 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 δ 18 O and δ 2 H and the non-conservative SO 4 were negatively correlated to alkalinity, temperature, δ 13 C and δ 34 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.

Precipitation and Surface Water Isotope Signals
Weekly measurements from 2006 to 2016 of δ 18 O and δ 2 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 δ 18 O (weighted average ± standard error of −13.2 ± 0.1 ) and −170 to −44 for δ 2 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 δ 18

Groundwater Source Discrimination
Groundwater from Central wells with an identified surface connection (UGC1 and LC1) also had δ 18 O values (−14.3‰ and −14.2‰, respectively) and δ 2 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-HCO3 water type and strongly negative δ 18 O (−17.3‰ to −16.8‰) and δ 2 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 δ 13 C (−12.9‰ to −10‰; Figure 5), δ 34 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.

Groundwater Source Discrimination
Groundwater from Central wells with an identified surface connection (UGC1 and LC1) also had δ 18 O values (−14.3 and −14.2 , respectively) and δ 2 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 3 water type and strongly negative δ 18 O (−17.3 to −16.8 ) and δ 2 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 δ 13 C (−12.9 to −10 ; Figure 5), δ 34 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.4).This "Shallow" source water also has similar δ 13 C values (−15.9 to −13.1 ; Figure 5), δ 34 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 δ 18 O and δ 2 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 of δ 13 C (Table 6) 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).

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 δ 18 O and −131 and −124 for δ 2 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 δ 13 C (−13.3 and −12.9 ) and δ 34 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 δ 13 C or δ 34 S.This lack of acceptable model fit likely is a product of increased alkalinity from water-rock interaction and decreased SO 4 because of reducing conditions.For example, a relatively low SO 4 concentration and an elevated δ 34 S isotope ratio in groundwater from GC1 indicate likely loss of SO 4 ; and this well contained noticeable H 2 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.

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) δ 18 O (−17.4 to −16.4 ), δ 2 H (−132 to −126 ), δ 13 C (−12.8 to −11.9 ) and δ 34 S (19 to 31.3 ) values and were a mixed cation (Ca and Na)-HCO 3 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 δ 34 S values indicative of SO 4 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 δ 13 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 (Tables 5 and 6: GW1, GW2, GW3, GW4 and UGW2).

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) (Figures 4-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 (Tables 5 and 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), δ 13 C and δ 34 S values (Figures 5 and 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).
groundwater flow from the intermediate-to-deep Central well locations moving westward.Testing of the water isotope mixtures (Table 5) by inverse modeling of δ 13 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 (Tables 5 and 6: GW1, GW2, GW3, GW4 and UGW2).

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) (Figures 4-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 (Tables 5 and 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), δ 13 C and δ 34 S values (Figures 5 and 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 (δ 18 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 (δ 18 O of −13.2‰) and Basin surface water (δ 18 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 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 (δ 18 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 (δ 18 O of −13.2 ) and Basin surface water (δ 18 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 δ 13 C and alkalinity pairs and most of the δ 34 S and SO 4 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 (Figures 4 and 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 (Tables 5 and 6).Groundwater from Western wells also showed elevated alkalinity and low SO 4 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.

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 δ 18 O, δ 2 H, δ 13 C and δ 34 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.

Figure 1 .
Figure 1.The South Fork Palouse River Basin of north central Idaho and eastern Washington (Wash.),USA.The South Fork Palouse River Basin is a subbasin in the Palouse River Basin.Included in the figure are the study's well locations shown as blue circles, a U.S. Network for Isotopes in Precipitation (USNIP) station (WA24) shown as a red square, a high elevation snow collection site shown as a red circle, the location of Pullman, Washington (P) and Moscow, Idaho (M) and the lineament of a geologic cross-section, A-A' (Figure 2).

Figure 1 .
Figure 1.The South Fork Palouse River Basin of north central Idaho and eastern Washington (Wash.),USA.The South Fork Palouse River Basin is a subbasin in the Palouse River Basin.Included in the figure are the study's well locations shown as blue circles, a U.S. Network for Isotopes in Precipitation (USNIP) station (WA24) shown as a red square, a high elevation snow collection site shown as a red circle, the location of Pullman, Washington (P) and Moscow, Idaho (M) and the lineament of a geologic cross-section, A-A' (Figure 2).

Figure 2 .
Figure 2. A west-to-east cross section of the South Fork Palouse River Basin from Pullman, Washington, to Moscow, Idaho (transect A-A' outlined in Figure 1) and select wells sampled during this study to indicate typical well depths and screened intervals in the three well groupings-Western, Central and Eastern.

Figure 2 .
Figure 2. A west-to-east cross section of the South Fork Palouse River Basin from Pullman, Washington, to Moscow, Idaho (transect A-A' outlined in Figure 1) and select wells sampled during this study to indicate typical well depths and screened intervals in the three well groupings-Western, Central and Eastern.

Hydrology 2018, 5 , 18 Figure 3 .
Figure 3. Piper diagram of major ion chemistry for groundwater samples collected from wells in the South Fork Palouse River Basin.

Figure 3 .
Figure 3. Piper diagram of major ion chemistry for groundwater samples collected from wells in the South Fork Palouse River Basin.
O and −15 to −5 for δ 2 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 δ 18 O and −146 ± 8 for δ 2 H compared to the modeled values of −15.3 for δ 18 O and −115 for δ 2 H for precipitation at the same elevation.Existing observations for δ 18 O values (δ 2 H not available) for mid-Basin creeks/rivers that convey water from the Palouse Range, such as Missouri Flat Creek (weighted average δ 18 O of −14.7 [31]) and the South Fork Palouse River (weighted average δ 18 O of −15.3 between October and March and −14.5 between April and September [34]), indicate δ 18 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 δ 18 O and δ 2 H values that may lessen the utility of the δ 2 H values for the inverse modeling.Hydrology 2018, 5, x FOR PEER REVIEW 10 of 18

Figure 4 .
Figure 4. Measured and modeled values of δ 18 O and δ 2 H for precipitation, surface water and groundwater in the South Fork Palouse River Basin.Measured values for summer (SP) and winter (WP) precipitation from USNIP station WA24, near Pullman, Washington, are plotted as weighted averages.Also included are published values for surface water (SW), analyzed values for snow (SM) atop the Palouse Range and modeled values for precipitation (MP) from atop the Palouse Range.

Figure 4 .
Figure 4. Measured and modeled values of δ 18 O and δ 2 H for precipitation, surface water and groundwater in the South Fork Palouse River Basin.Measured values for summer (SP) and winter (WP) precipitation from USNIP station WA24, near Pullman, Washington, are plotted as weighted averages.Also included are published values for surface water (SW), analyzed values for snow (SM) atop the Palouse Range and modeled values for precipitation (MP) from atop the Palouse Range.

Hydrology 2018, 5 , 18 Figure 5 .
Figure 5. Measured values of δ 13 C and alkalinity (as CaCO3) concentrations for groundwater samples collected from wells in the South Fork Palouse River Basin.

Figure 6 .
Figure 6.Measured values of δ 34 S and sulfate (SO4) concentrations for groundwater samples collected from wells in the South Fork Palouse River Basin.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 δ 18 O and δ 2 H values for these shallower wells with Ca-HCO3 type groundwater (WE1, WE2, LC2, WC1, WC2 and WC3) ranged from −16.1‰ to −15.3‰ for δ 18 O and −121‰ to −117‰ for δ 2 H and had well screens ranging from 12-228 m below land surface.This limited variation in relatively negative water isotope values for Ca-HCO3 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

Figure 5 . 18 Figure 5 .
Figure 5. Measured values of δ 13 C and alkalinity (as CaCO 3 ) concentrations for groundwater samples collected from wells in the South Fork Palouse River Basin.

Figure 6 .
Figure 6.Measured values of δ 34 S and sulfate (SO4) concentrations for groundwater samples collected from wells in the South Fork Palouse River Basin.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 δ 18 O and δ 2 H values for these shallower wells with Ca-HCO3 type groundwater (WE1, WE2, LC2, WC1, WC2 and WC3) ranged from −16.1‰ to −15.3‰ for δ 18 O and −121‰ to −117‰ for δ 2 H and had well screens ranging from 12-228 m below land surface.This limited variation in relatively negative water isotope values for Ca-HCO3 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

Figure 6 .
Figure 6.Measured values of δ 34 S and sulfate (SO 4 ) concentrations for groundwater samples collected from wells in the South Fork Palouse River Basin.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 δ 18 O and δ 2 H values for these shallower wells with Ca-HCO 3 type groundwater (WE1, WE2, LC2, WC1, WC2 and WC3) ranged from −16.1 to −15.3 for δ 18 O and −121 to −117 for δ 2 H and had well screens ranging

Figure 7 .
Figure 7. Representation of water sources, subsystems and flow paths through the aquifer system in the South Fork Palouse River Basin.

Figure 7 .
Figure 7. Representation of water sources, subsystems and flow paths through the aquifer system in the South Fork Palouse River Basin.

Table 1 .
Descriptions of wells east of Moscow (Eastern) and wells in and near Moscow (Central), Idaho, sampled as part of this study (Figure

Table 2 .
Descriptions of wells in and near Pullman (Western), Washington, sampled as part of this study (Figure1).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).

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

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.

Table 5 .
Inverse modeling (Equation (1)) results of δ 18 O and δ 2 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.

Table 6 .
Inverse modeling results for δ13C 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.