Simulating water-quality trends in public-supply wells in transient flow systems.

Models need not be complex to be useful. An existing groundwater-flow model of Salt Lake Valley, Utah, was adapted for use with convolution-based advective particle tracking to explain broad spatial trends in dissolved solids. This model supports the hypothesis that water produced from wells is increasingly younger with higher proportions of surface sources as pumping changes in the basin over time. At individual wells, however, predicting specific water-quality changes remains challenging. The influence of pumping-induced transient groundwater flow on changes in mean age and source areas is significant. Mean age and source areas were mapped across the model domain to extend the results from observation wells to the entire aquifer to see where changes in concentrations of dissolved solids are expected to occur. The timing of these changes depends on accurate estimates of groundwater velocity. Calibration to tritium concentrations was used to estimate effective porosity and improve correlation between source area changes, age changes, and measured dissolved solids trends. Uncertainty in the model is due in part to spatial and temporal variations in tracer inputs, estimated tracer transport parameters, and in pumping stresses at sampling points. For tracers such as tritium, the presence of two-limbed input curves can be problematic because a single concentration can be associated with multiple disparate travel times. These shortcomings can be ameliorated by adding hydrologic and geologic detail to the model and by adding additional calibration data. However, the Salt Lake Valley model is useful even without such small-scale detail. Published 2014. This article is a U.S. Government work and is in the public domain in the USA. Groundwater published by Wiley Periodicals, Inc. on behalf of National Ground Water Association.

Introduction

Interpreting water-quality trends in hydrologically complex and (or) transient groundwater-flow systems is challenging. Groundwater modeling can help unravel complex interaction among variables by providing a framework for that interaction that is consistent with the physics of groundwater flow. In particular, a model can be used to compute spatial and time-series quantities that can be used to explain measured water-quality changes. The purpose of this study is to (1) evaluate the use of an existing transient groundwater-flow model to answer water-quality questions, (2) develop and test model-derived time-dependent variables in relation to measured water-quality trends, and (3) evaluate the usefulness of tracer calibration data for the purpose of explaining trends. The approach used in this study is to adapt an existing groundwater-flow model for advective transport simulation by adjusting effective porosity to match tritium concentrations in samples from wells. The model is then used to estimate changes over time in groundwater age and source area and to relate these changes to measured trends in dissolved solids concentrations in Salt Lake Valley, Utah. This study draws on previous work by the authors and will help guide future transport modeling of regional-scale groundwater-quality trends in transient flow systems.

Recent editorials/articles by Konikow (2011), Voss (2011)a), and Voss (2011)b) raise questions about appropriate complexity of transport models, including the role of macrodispersion (especially at the regional scale), meaningful values of effective porosity, and simultaneous calibration of flow and transport properties. Macrodispersion is a primary transport mechanism in groundwater; however, several studies of regional-scale transport have used advection-only to provide a first-order estimate of solute concentrations. For example, Sanford (2011) and Green et al. (2010) simulated groundwater age using only advection. Simulations of advective groundwater transport have been successfully used to interpret trends of solute concentrations (Kauffman et al. 2001; McMahon et al. 2008) and assess groundwater age (Wellman et al. 2012) in steady flow. One reason these efforts were successful is that groundwater at receptors such as streams or wells at the regional scale integrates a wide distribution of flow paths, source areas, and groundwater ages. These successes suggest that simulation of advection without consideration of macrodispersion would adequately represent large spatial scale transport in this regional system.

Effective porosity, which converts inter-cell flux to velocity (Konikow 2011), is required to adapt a groundwater-flow model for transport simulation. It accounts for the fact that transport does not occur evenly throughout a cross section of flow and does not include pore spaces that are not connected. For this reason prior estimates of effective porosity vary widely (Konikow 2011) and often are not available. In this study, effective porosity is adjusted as part of the calibration process so that simulated tritium concentrations match tritium concentrations measured in samples from pumping wells. Advective transport models have been calibrated to tritium/helium age dates (Reilly et al. 1994; Sheets et al. 1998; Murphy et al. 2011); however, fewer studies have relied solely on tritium concentrations.

Ideally all calibration data (hydraulic heads, streamflow, and tritium concentrations) could be simulated by a single flow and transport model and parameters of the flow model (e.g., hydraulic conductivity and storage coefficient) adjusted simultaneously with effective porosity. However, as Voss (2011)a, 2011b) points out, this may not be possible because of disparate transport processes—heads change instantaneously and uniformly throughout an aquifer, whereas mass changes discretely and nonuniformly through permeable pathways in an aquifer. In this study, the flow model is not modified; groundwater velocity is calibrated by adjustment of effective porosity.

Groundwater-flow systems are never really in steady state, and there are relatively few studies that consider transient flow in analyzing water-quality trends (Furlong et al. 2011). In transient flow, age distributions become time dependent (Manning et al. 2012) and the system's response to solute input also depends on time. Time-dependent source areas (Rock and Kupfersberger 2002; Masterson et al. 2004) possibly contribute different solutes to groundwater, resulting in a complex mixture of water at the receptor (Starn and Brown 2007). Simulation of transient flow can be used to unravel these processes by showing possible changes in source areas and groundwater ages. Short-term recharge variations cause transient flow, but the timescale of groundwater flow often is much longer, in which case the effect is minimal over time (Reilly and Pollock 1996). However, long-term recharge or discharge variations (longer than the timescale of groundwater flow) can lead to significant long-term changes in water quality. Transient flow caused by pumping variations is simulated in this study as a fundamentally important characteristic of the groundwater system in Salt Lake Valley.

Models are built for specific purposes, and using a model for a purpose other than for which it was intended can lead to misleading results. It is a significant question whether groundwater models built to answer questions about water levels and flow rates can be used to explain possible causes of water-quality trends. The model used in this study provides an example; it was created to answer questions about water levels and flow rates in Salt Lake Valley, Utah (Lambert 1995).

Salt Lake Valley Geology and Water Quality

In the Salt Lake Valley, Utah (Figure 1), the basin-fill aquifer is bounded by mountains to the east, south, and west and by Great Salt Lake to the north. Quaternary-age basin fill consists of unconsolidated to semi-consolidated sediments deposited in climate-driven cycles of lacustrine, deltaic, and alluvial environments (Stolp 2007). Sediments average 600 m thick in interbeds of clay, silt, sand, and gravel and lenses of sand and gravel. Sediments originated from the adjacent mountains and are generally coarser and less well sorted near the valley wall and finer near the center of the valley. Discontinuous fine-grained layers confine groundwater in the central part of the basin. These layers are not present near the basin edges, and recharge occurs there from losing stream reaches, infiltrating precipitation, and flow through fractured consolidated rock in the adjacent mountains. Recharge through the confining layers from an overlying shallow unconfined aquifer also occurs, but to a lesser degree. The simulated aquifer includes tertiary-age semi-consolidated sediments at its base where they are permeable enough to yield water to wells (Stolp 2007). Groundwater flow is from the recharge areas primarily to the Jordan River in the center of the valley (Figure 1) and to pumping wells. Groundwater from wells supplied about 29% of the water used for public supply in the basin in 2005 (Thiros and Spangler 2010). While most of the groundwater pumped is used for drinking, some is used for agricultural and industrial purposes (Burden 2009). Groundwater also discharges through springs and evapotranspiration.

Figure 1

Salt Lake Valley, Utah. Solid black line is extent of active model grid. Tan area is water ≤500 mg/L dissolved solids in 1998.

Thiros and Spangler (2010) documented increasing dissolved solids in wells over a wide area using long-term water-quality records from the 1930s to the 2000s. In the Salt Lake Valley aquifer, a plume of water containing less than 500 mg/L dissolved solids is surrounded by groundwater of higher dissolved solids concentrations (Figure 1). Note that here “plume” refers to water that is of better quality for drinking, not lower as the term is usually used. Sources of dissolved solids include mineral dissolution from aquifer sediments, infiltration of de-icing chemicals, and concentration by evaporation in surface-water sources to the Jordan River. Some parts of the aquifer contain older water that has been in contact with aquifer minerals for centuries or millennia, resulting in high-dissolved solids concentrations (Thiros and Spangler 2010). Carbon-14 and Helium-4 data support very old ages in some wells in Salt Lake Valley (Stephen Hinkle USGS 2012, unpublished data). Dissolved solids from mineral dissolution occur in several broad categories of water, including water containing: calcium, bicarbonate, and sulfate from dissolution of shale minerals, sulfate from oxidation of sulfide minerals; and sodium and chloride from desiccated paleolakes.

Dissolved solids in water from other wells, however, have fluctuated within a small range. The plume appears to be associated with water coming from terrain underlain by relatively erosion-resistant rock in the adjoining mountains. Sources of water to these wells probably include recharge from snowmelt and subsurface inflow from adjacent mountain block areas.

Population in the valley steadily increased from the 1930s to more than 1 million people in 2011. Although pumping at individual wells has increased or decreased over time, total withdrawals from the aquifer leveled off after about 1997. During the period of population increase, many wells for which there are long-term data had changes in dissolved solids concentration. Changes in dissolved solids concentrations over time can be related to changing human activities such as increased use of de-icing chemicals and changes in the location and magnitude of groundwater withdrawals from wells. Groundwater withdrawals can induce flow from new source areas by enhancing downward flow from areas affected by recent anthropogenic activity or horizontal flow from areas affected largely by natural mineral processes.

Groundwater withdrawals and changes in storage were relatively constant from 1964 to 1968, and a steady-state model for that time period was manually calibrated to measured water levels, groundwater flow to the Jordan River, and vertical hydraulic gradients in the aquifer (Lambert 1995). The steady-state simulation was used as the initial condition for a transient simulation in which recharge and pumping rates were varied at their estimated annual rates from 1969 to 1991. The transient model was manually calibrated to measured water-level changes and groundwater flow to the Jordan River (Lambert 1995). Lambert (1996) extended the original simulation for part of the modeled area to include the time period 1935 to 1964. Total water use and public-supply annual withdrawals were relatively constant from 1997 to 2009 (Burden 2009), although the withdrawal at individual wells may have changed, causing local water-level fluctuations. Stolp (2007) updated groundwater withdrawals and recharge in the model to reflect average 1997 to 2001 conditions.

The models constructed by Lambert (1995), 1996) and modified by Stolp (2007) cover 1,152 km2 of the Salt Lake Valley. The model grid cells are 563.27 m on each side in 94 rows and 62 columns. The aquifer is divided into seven layers; the top two layers represent a shallow unconfined aquifer and an underlying confining unit in the center of the valley. The thickness of the two layers is variable. Layers 3 through 7 represent the principal aquifer and are 46, 46, 46, 61 m, and greater than or equal to 61 m thick, respectively. Layer 7 is a maximum of 460 m thick in the deepest parts of the basin. Boundary conditions include recharge from precipitation, losing streams, and the mountain block and discharge through pumping wells, gaining stream reaches, and evapotranspiration. The base of the aquifer is a zero-flow boundary. More details on model properties, boundaries, and calibration are given by Lambert (1995), 1996) and Stolp (2007).

The original model by Stolp (2007) was modified from Lambert (1995), 1996) for MODFLOW-2000. For this project, the model was modified for use with MODFLOW-NWT (Niswonger et al. 2011) which necessitated changing from the original block-centered flow package to the layer property flow package (Harbaugh 2005). Also, for this work, the multi-node well package for MODFLOW (Konikow et al. 2009) was used to better simulate the public-supply wells. No other changes were made to the model. The differences between the flow budgets in all stress periods between the original model and the modified model were very small and were considered negligible.

Approach

Groundwater age distributions at observation wells were simulated using transient particle tracking (MODPATH; Pollock 2012) rather than grid-based methods (e.g., MT3DMS; Zheng and Wang 1999) for efficiency and to limit numerical dispersion. Age distributions were transformed to tritium concentration breakthrough curves (BTCs) using convolution-based particle tracking (CBPT; Robinson et al. 2010). Convolution transforms the impulse response of a system to a response of the system that reflects time-varying input (Małoszewski and Zuber 1982). Flow-model grid cells were relatively large, and an analytical equation was used to compute velocity near pumping wells as described by Zheng (1994) and adapted for use with MODPATH by Starn et al. (2012, 2013). Although CBPT is only applicable to steady-state flow systems, recent work extends CBPT to time-varying flow by assuming steady flow over small time intervals (Srinivasan et al. 2011; Starn et al. 2013).

Calibration of the Salt Lake Valley flow model is well documented (Lambert 1995, 1996; Stolp 2007), and no parameters of that model were changed in order to maintain the integrity of its calibration. Also, adjusting only porosity avoids difficulty with numerical instability in the flow model caused by the addition of correlated parameters (porosity and hydraulic conductivity).

To adapt the groundwater-flow model for water-quality simulation, the effective porosity is required, which converts the inter-cell flux calculated by the model to velocity. Effective porosity was adjusted using the inverse modeling program PEST++ (Welter et al. 2012) to minimize the sum of squared differences between CBPT simulated equivalent concentration and tritium concentrations measured in samples from pumping wells. In this study, tritium concentrations and tritium/helium apparent ages were available for model calibration, although the ages were not used in the final calibration.

Tracer Concentration Data

Tritium concentrations from 80 wells are available for model calibration, mainly in the area where dissolved solids changes have been measured. Some of the wells have been sampled multiple times for a total of 122 concentrations; 61 of these values are paired with a tritium/helium apparent age date. The timescale of tritium persistence in this aquifer system (decades) is comparable to the temporal scale of changes in dissolved concentrations; therefore, calibration to these data should serve to decrease uncertainty in the model. Most of the tritium data are in the USGS National Water Information System database (http://waterdata.usgs.gov/ut/nwis/qw/). Additional data are contained in the papers by Manning (2002), Thiros (2003), and from USGS sites reported by University of Utah and Lamont-Doherty Earth Observatory. Tritium/helium apparent ages are available from the papers by Thiros (2003) and Thiros and Manning (2004).

Tritium Input Functions

Two tritium input functions were used based on a precipitation-weighted tritium input curve for Salt Lake City (Manning et al. 2005). For 1953 to 1962, correlation with the Ottawa record was used, for 1963 to 1984 Salt Lake City data are used, and for 1985 to 2000, the Vienna record was used (Manning et al. 2005). For this work, 8 TU (tritium unit) was used as the value of pre-bomb and post-2000 tritium concentrations, as estimated by Thatcher (1962) for the central United States. This tritium input function is the estimated input to land surface, not to the water table. Analysis by Manning et al. (2005) provides evidence for generally negligible tritium residence time (<5 years) in the unsaturated zone in this aquifer system.

A second tritium input function was used for the input of water from mountain block recharge to the alluvium to account for tritium decay before the water reaches the alluvium. This study assumed exponential mixing (Małoszewski and Zuber 1982) with a mean residence time of 12.5 years in the mountain block (Manning and Solomon 2005). An assumption of exponential mixing was consistent with simplified flow system geometry in which mountain block recharge discharges into the alluvial basin at the bedrock/alluvium boundary. This input function was applied at specified flux nodes in layers 3 and 4 that correspond to the inflow of mountain block recharge into the primary public-supply depths of the aquifer.

Inverse Modeling

The inverse model PEST++ (Welter et al. 2012) was used to estimate the magnitude and distribution of effective porosity, given the velocity field from the flow model and the measured tritium concentrations. PEST++ uses nonlinear regression to minimize least-squared residuals with respect to the observations. Porosity distribution was calculated during model calibration using 37 pilot points distributed in layers 3 to 5, with the density of pilot points roughly corresponding to the density of tritium measurements. Few tritium data were available for the shallow (1 to 2) and deep (6 to 7) layers, and single porosity values were estimated to represent each of those layer groups for a total of 39 porosity parameters. The model calibration was not stable without applying regularization. Prior-known-value regularization is difficult because little prior data are available on effective porosity and because structural errors (with respect to transport processes) are incorporated into the porosity parameter. To avoid choosing a prior preferred value of porosity, 404 regularization equations were included in the calibration that penalized the solution for differences between nearest-neighbor pilot point values weighted inversely proportional to their separation distance. In other words, regularization imposed smoothness on the solution. A kriging algorithm then distributed pilot point values to grid cells prior to running the model. Tritium concentrations were weighted in the regression based on the inverse square of the product of the coefficient of variation and the simulated concentration. A coefficient of variation of 0.1 and a lower standard deviation threshold of 0.1 TU was used here in order not to lose the effect of low concentrations in the regression (Hill and Tiedeman 2007).

Effective Porosity Estimated from Tritium Concentration

Before calibration, forward simulations with low, medium, and high effective porosities were run (Figure 2), which demonstrated a possible dependence on initial effective porosity estimates and a generally low sensitivity of the model to the data. Thus, the estimated effective porosity is not unique. Based on a comparison of simulated BTCs to a representative subset of tritium and age data, an initial porosity of 0.10 was judged the best initial value. This choice maximized the sensitivity of the model to the data, allowed a reasonable estimated distribution of porosity, and produced the lowest sum-of-squared errors among the starting values tried. In inverse models using different initial values, the results had these similarities: (1) the effective porosity was always lower than the expected bulk porosity, (2) effective porosity in the center of the fresh water plume was always lower than surrounding values, and (3) small relative adjustments (a few percent) in effective porosity reduced the sum-of-squared errors by about 50%. Calibration also was attempted with data sets that included the tritium/helium age dates. Including ages in the calibration did not alleviate nonuniqueness, and ages were not used in the final calibration. Tracer age only corresponds to mean particle age if the flow system is approximated by piston flow, and in a transient flow system complicated by pumping, this was not deemed likely.

Figure 2

Sensitivity of simulated tritium breakthrough to different constant uniform porosity values. Solid line, effective porosity 0.002; long dashes, 0.02; short dashes, 0.20. Crosses are measured tritium concentrations.

Although most simulated BTCs fit the relative magnitude (Figure 3A) and trend of measured tritium (Figure 3C), the fit at some wells was not good (Figure 3B). Local variations in hydraulic properties and the coarseness of simulated pumping history (as in Figure 3B) were probably the largest contributors to poor fit. Existing porosity estimates for the east side of the valley (Freethey et al. 1994) show an increase from the east edge of the valley toward the Jordan River with zones of higher porosity beneath Big and Little Cottonwood Creeks. The estimated effective porosity distribution from using pilot points was similar, although the magnitudes were relatively low, ranging from 0.03 to 0.15. The smoothness criteria prevented a wider range of values. This range probably was lower than Freethey et al. (1994) estimates because the degree of connectedness of porosity was less than expected and because it accommodates lack of small-scale detail in the flow model.

Figure 3

Simulated tritium breakthrough at three wells. Crosses are sampled tritium concentrations. Depths of screen below land surface. (A) 186 to 290 m; (B) 43 to 116 m; (C) 52 to 77 m. Well B is affected by pumping.

Simulated Factors Related to Dissolved Solids Trends

The area of the Salt Lake Valley aquifer containing water less than 500 mg/L dissolved solids decreased from 1964 to 1998, and the concentration of dissolved solids in many long-term observation wells increased over the same period (Thiros and Spangler 2010). The largest increases are along the southern margin of an area surrounding Big and Little Cottonwood Creeks (shaded area in Figure 4). One hypothesis is that increased pumping has increased the proportion of water recharged from surface sources (as opposed to mountain block recharge) in these wells. Three model-derived explanatory factors that can be used to test the hypothesis were developed using the calibrated model. The first factor was the fraction of each water source contributing to a well over time created by tagging reverse-tracked particles starting in a well with their ending locations (corresponding to the source-water area). The second factor was the fraction of water supplied to the well from each model layer calculated by the Multi-Node Well package (Konikow et al. 2009) in MODFLOW. The third factor was the mean simulated age created by taking the mean age of particles reverse-tracked from the well.

Figure 4

Range of changes in dissolved solids in mg/L and extent of water with dissolved solids ≤500 mg/L. Ranges are depicted for wells for which time series were available. Additional wells were sampled to determine extent in 1964 (tan area) and 1998 (enclosed by solid black lines) of water with dissolved solids ≤500 mg/L.

The wells with the highest range of increased dissolved solids all had simulated sources area that changed over time with respect to the plume location (for example, Well D in Figure 5A). Pumping apparently increased in well D or in a nearby well(s) around 1980. At the same time, the mean age decreased (Figure 5B), the proportion of water coming from surface sources increased (Figure 5C), and water traveled to the well through shallower model layers (Figure 5D). Source areas in the center of the plume (Well E in Figure 6A) changed with time but remained within the plume. The initial age of water in this source area was much younger than in well D (Figure 6B) and also decreased as the proportion of surface sources increased (Figure 6C).

Figure 5

Results of simulation for Well D (water quality shown in Figure 4). Panel A, source area of particles released from well in year indicated by color; B, simulated mean age, measured dissolved solids dashed line with circles; C, source of water (MBR is mountain block recharge); D, layer contributions.

Figure 6

Results of simulation for Well E (water quality shown in Figure 4). (A) Source area of particles released from well in year indicated by color; (B) simulated mean age, measured dissolved solids dashed line with circles; (C) source of water (MBR is mountain block recharge); (D) layer contributions.

Broad spatial patterns of sources and mean ages correspond predictably to the measured spatial shifts in dissolved solids concentration. Changes in source area and mean age can be mapped by placing particles at the center of each model cell at specific times (1950 and 2020 shown in the Figures 7 and 8) and reverse tracking them to their source. The broad spatial trend is for less water from mountain block recharge and more water from leaking canals and river reaches (Figure 7). In the center of the model area, more water is coming from canal and river leakage; along the model edges, less water is coming from mountain sources over time. Similarly, the mean age (Figure 8) in model layer 3 decreases over time in the southern two-thirds of the model area, consistent with a downward shift of recently recharged water originating from surface sources.

Figure 7

Simulated change in sources of water to model layer 3 from 1950 to 2020. Layer 3 is at typical depth for public water supply.

Figure 8

Simulated change in mean age of water to model layer 3 from 1950 to 2020. Layer 3 is at typical depth for public water supply.

Conclusions

Models need not be complex to be useful. A complex model in this sense would be one that simulated macrodispersion, chemical reactions in the aquifer, and detailed aquifer heterogeneity, which this model does not do. Although prediction of solute concentration trends at individual wells in complex systems remains a challenging problem, this work shows that an existing flow model of Salt Lake Valley, Utah, adapted for transport simulation was useful in explaining measured dissolved solids trends. At individual wells, the source fraction, mean age, and depth (model layer) change in concert with dissolved solids but not always in an easily predictable way.

The model provided a physically consistent framework of model and data that support the hypothesis that water produced by wells is increasingly younger with higher proportions from surface sources as a result of pumping. The most useful predictive variables for changes in dissolved solids are the changes in source area location and mean age. These variables change over time, and transient changes in the flow system must be considered in order to explain trends. An accurate estimation of velocity (through inverse modeling of effective porosity in this case) is important for assessing the timing of those changes.

Changes in mean age and source areas can be mapped across the model domain to extend the results at wells to the entire aquifer area to see where changes are expected to occur. However, constituents such as dissolved solids are affected by many geochemical processes. To make full use of model-derived factors, one must incorporate how recharge solute chemistry changes at the land surface and how rock-water interactions alter the composition of groundwater.

Calibration to tracer data can improve correlation between model-derived factors and measured water-quality trends. However, effective porosity estimated by inverse modeling may not be optimal despite reasonable spatial patterns and substantial decreases in the sum-of-squared error. Better estimation of effective porosity by using concentrations of multiple tracers over time can increase confidence in porosity estimates. Despite uncertainty in estimated effective porosity, the groundwater-flow model of the Salt Lake Valley aquifer, designed for simulating large-scale flow, proves useful for interpreting groundwater-quality trends.

Atmospheric tracer concentrations (such as tritium) collected over large areas in heterogeneous aquifers and (or) in transient flow systems exhibit variations due to spatial and temporal changes in tracer inputs, in tracer transport in unsaturated and saturated zones, and in pumping stresses at sampling points. Additionally, for non-monotonic tracers such as tritium, the presence of two-limbed BTCs can be problematic, especially in light of these variations, because a single concentration can be associated with multiple disparate travel times. Time-series data that span the peak and (or) define at least one limb of the BTC are most useful. Use of multiple tracers alleviates this problem.

Calibrated effective porosity compensates for lack of small-scale detail in the flow model. This shortcoming can be ameliorated by adding hydrologic and geologic detail to the model, and by adding additional calibration data. However, the Salt Lake Valley model is useful even without such small-scale detail.

This study highlights the usefulness of atmospheric tracer data (and the need for long-term commitment to such data) and ancillary chemical data for illustrating and explaining groundwater-quality trends. However, more work is needed to understand the relation among numerical method, calibration method, and basin-scale predictions of groundwater quality. Areas for possible future work include providing examples of the effect of the types and amounts of tracer data on estimation of porosity and hydraulic conductivity in simultaneous calibration of flow and transport models, in particular in relation to basin-scale predictions.