Background and motivation

Model uncertainty in transport simulations is affected by how transport parameters are represented in the model and by the number and type of data available for model calibration. Typically, modelers attempt to reduce uncertainty by adding hydrologic and geologic detail to the model and/or by adding additional calibration data. In this paper, model results are compared by using simple and complex models calibrated with either tritium concentrations or tritium concentrations plus interpreted tritium/helium apparent ages. Tritium concentrations above background levels have persisted in the Salt Lake Valley aquifer, Utah (USA) for decades, a time period over which changes in total dissolved solids concentrations also have been observed (Starn et al. 2014); therefore, one would expect that calibration to tritium data should serve to decrease uncertainty in model predictions of changes in water quality.

Ideally, all calibration data are simulated by a single flow and transport model and parameters of the flow model (e.g., hydraulic conductivity and storage coefficient) should be adjusted simultaneously with effective transport parameters; however, this may not be possible because of disparate transport processes in groundwater flow (diffusion) and solute transport (advection; Voss 2011a, b). Perturbations in head propagate evenly throughout a model domain, whereas solute concentrations change discretely and non-uniformly through permeable pathways in an aquifer (Voss 2011a, b). This discrepancy can be more apparent in transient groundwater systems because a given observation may be affected by a head perturbation long before it is affected by a change in solute concentration. Only one parameter, effective transport porosity, is required to convert a groundwater-flow model for advective transport simulation; however, prior estimates of effective porosity vary widely and often are not available for most sites (Konikow 2011). Effective transport porosity values have been determined through calibration to tritium/helium apparent ages (Reilly et al. 1994; Sheets et al. 1998; Murphy et al. 2011), and, less often, directly to tritium concentrations (Herweijer et al. 1985; Engesgaard et al. 1996; Starn et al. 2014).

Large and long-term changes in recharge or discharge relative to groundwater residence times can lead to significant long-term changes in water quality (Reilly and Pollock 1996). Relatively few studies have considered short-term (annual) transient flow in analyzing water-quality trends (Furlong et al. 2011; Starn et al. 2014). In transient flow, age distributions are time dependent (Manning et al. 2012). Time-dependent source areas (Rock and Kupfersberger 2002; Masterson et al. 2004; Starn and Brown 2007; Starn et al. 2014) can contribute different solutes to groundwater, resulting in complex mixtures of water at pumped wells.

Compared to 20 years ago, automated parameter estimation methods are now in common use for groundwater flow and transport models, largely because of public-domain automatic model calibration software such as UCODE (Poeter et al. 2005) and PEST (Doherty 2010, 2014). The use of these methods is discussed at length by Hill and Tiedeman (2007). Each of these software approaches favors different model parameterization strategies. UCODE favors parsimonious models in which the modeler applies geologic knowledge to create piecewise-constant parameter zones. On the other hand, PEST favors highly parameterized models in which parameter fields can vary continuously, thus allowing nuances in the data to possibly be expressed in the solution. In this paper model predictions made by using each of the approaches are explored. Also explored is the effect of having different sets of observations for model calibration with each of these approaches.

In order to improve a model, a modeler can increase the complexity of the model in the hope of capturing additional physical processes that affect model results. Another way to improve a model is to add more or different types of data. In principle, a highly parameterized model should be better able to match the data because it has more parameters, and a model with more and different kinds of data should be able to capture more of the real physical variability in the system being modeled. Although this is not an exhaustive comparison of modeling approaches, the goal in this paper is to contribute to the constructive dialog on the relative merits of complex models and types of observations employed in model calibration. The question posed is the following: “Is added model complexity as useful for predictions as adding different types of observations?”

Approach

The approach used in this study is to convert an existing groundwater flow model to simulate advective transport in the Salt Lake Valley, Utah (Fig. 1), by adjusting effective porosity to match tritium concentrations in samples from wells and, alternatively, interpreted tritium/helium apparent ages (Starn et al. 2014). The modeling procedure includes simulation of groundwater flow (MODFLOW-NWT; Niswonger et al. 2011), simulation of advective transport (MODPATH; Pollock 2012), and a suite of pre- and post-processing steps to perform convolution-based particle tracking (CBPT; Robinson et al. 2010; Green et al. 2010; Srinivasan et al. 2011; Starn et al. 2012, 2013, 2014). To simplify interpretation, parameters of the transient groundwater flow model are held constant and only advective transport is considered. This is justified in part because parameters for groundwater flow may not correspond directly to parameters for transport. Two scenarios of porosity distribution are calibrated using (1) a complex model with highly parameterized porosity distributions and (2) a simple model with parsimoniously zoned porosity distributions. These are described in detail in the next two sections.

Fig. 1
figure 1

Salt Lake Valley, Utah. Yellow area is extent of active model grid. Hydrograph locations (A, B, C) are shown with blue triangles

Full size image

Each porosity scenario is calibrated to two data sets for a total of four simulations that are compared. The complex scenario is calibrated to the “full data set,” which consists of 185 samples from 80 wells (Fig. 2) mainly in the area where increased dissolved solids has been observed. Some of the wells were sampled multiple times over a period of decades, yielding 122 observations of tritium concentration. Interpreted groundwater residence times from tritium/helium ratios are available from 52 of the 80 wells (63 observations because some wells were sampled multiple times). The two data sets are (1) tritium concentrations plus apparent ages and (2) only tritium concentrations without apparent ages.

Fig. 2
figure 2

Highly parameterized model, Salt Lake Valley, Utah. Yellow area is extent of active model grid. Circles represent pilot points and triangles represent tritium samples

Full size image

The simple scenario is calibrated to a “reduced data set” that was extracted from the “full data set” by including only wells for which a time series of tritium concentrations (more than one sample over time) and tritium/helium apparent ages were available. A further criterion was that the observation be obtained from wells where the well screens corresponded to at most two model layers (Fig. 3). Application of these criteria winnowed the data to 48 samples from 12 wells consisting of 31 observations of tritium concentration and 17 observations of interpreted groundwater residence times from tritium/helium ratios. The reason that the two scenarios were calibrated to different data sets is that a modeler would typically choose model complexity in part based on the amount of data available. One would be less likely to create the complex model if all one had available was the reduced data set. Also, the reduced data set did not contain sufficient information to calibrate the complex model.

Fig. 3
figure 3

Parsimonious model, Salt Lake Valley, Utah. Yellow area is extent of active model grid. Triangles represent tritium sample locations

Full size image

Most of the tritium data are in the USGS National Water Information System (NWIS) database (USGS 2015). Additional data were obtained from Manning (2002) and Thiros (2003), and tritium data from USGS sites reported by University of Utah and Lamont-Doherty Earth Observatory.

Hydrogeologic setting and simulation of groundwater flow

The basin-fill aquifer in the Salt Lake Valley, Utah (Fig. 1), is bounded by mountains and by the Great Salt Lake. Quaternary-age basin fill consists of unconsolidated to semi-consolidated sediments, which average 600 m thick with interbeds of clay, silt, sand, and gravel and lenses of sand and gravel (Stolp 2007). 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, where recharge occurs from losing stream reaches, infiltrating precipitation, and flow through fractured consolidated rock in the adjacent mountains.

The basin fill includes Tertiary-age semi-consolidated sediments at its base permeable enough to yield water to wells (Stolp 2007). Groundwater flows from the recharge areas to the Jordan River in the center of the valley (Fig. 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). 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, unpublished data, 2012).

Population in the valley has steadily increased from the 1930s to more than 1 million people in 2011. Although pumping at individual wells has been variable over time, total withdrawals from the aquifer leveled off after about 1997. During the period of population increase, many wells exhibited long-term increases in dissolved solids concentration. The increasing dissolved solids are related to the location and magnitude of groundwater withdrawals from wells. Increased groundwater withdrawals have induced flow from different source areas, by enhancing both downward flow from areas affected by anthropogenic activity and horizontal flow from areas containing mineralized water from natural processes.

A steady-state model of the Salt Lake Valley was manually calibrated to measured water levels, groundwater flow to the Jordan River, and vertical hydraulic gradients from 1964 to 1968, when groundwater withdrawals and changes in storage were relatively constant (Lambert 1995). Heads computed by the steady-state simulation were 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, and later extended for part of the modeled area to include the time period 1935 to 1964 (Lambert 1995, 1996). Stolp (2007) then updated this model to reflect average recharge and withdrawals during 1997 to 2001.

The model covers 1,152 km2 of the Salt Lake Valley. The model grid cells are 563 m on each side in 94 rows and 62 columns. The aquifer is divided into 7 layers; the top two layers represent the shallow unconfined aquifer and the underlying confining unit in the center of the valley. The thickness of each layer is variable. Layers 3–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, infiltration from losing streams and adjacent mountain blocks, and discharge through pumping wells, gaining stream reaches, and evapotranspiration. The base of the aquifer is a zero-flow boundary.

Simulation and calibration of advective transport

The 2007 groundwater flow model was modified for this study to use MODFLOW-NWT (Niswonger et al. 2011), which necessitated substituting the upstream weighting (UPW) flow package for the block-centered flow (BCF) package (Harbaugh 2005). The multi-node well (MNW2) package for MODFLOW (Konikow et al. 2009) also was used to better simulate pumping from long-screened public-supply wells. Differences in heads and flows between the 2007 model and modified models were very small and were considered negligible.

Although macrodispersion can be an important transport mechanism in groundwater, studies of regional-scale transport have used advective transport simulations to provide first-order estimates of solute concentrations with success (Sanford 2011). Advective transport simulations also have been used successfully to interpret trends of solute concentrations (Kauffman et al. 2001; McMahon et al. 2008) and to estimate groundwater residence times (Wellman et al. 2012) in steady flow.

Advective transport was simulated by using convolution-based particle tracking (CBPT; Starn et al. 2012, 2013, 2014). This method computes either resident solute concentration (for samples from unpumped wells and monitoring wells) or flux-weighted concentration (for pumped wells). The simulated equivalent of apparent age is taken to be the mean particle age. Multi-layer flow rates for flux-weighted concentrations are simulated by the MNW2 package. Flow within a cell containing a pumped well is simulated by using an analytical solution (Zheng 1994) and embedded into a pre-processor by Starn et al. (2012). Particle advection is simulated by using MODPATH and converted into concentration break-through curves by using the CBPT method cited earlier.

The computer code PEST++ (Welter et al. 2012) is used as an inverse simulator. PEST++ is used for both the complex and simple models. In the latter case, it is set up to estimate parameters defined in piece-wise constant zones. The highly parameterized inverse model uses pilot points, regularization, and truncated singular value decomposition as described in guidelines by Doherty and Hunt (2010a) and Welter et al. (2012). Porosity parameters were defined for each of 3 layer groups (layers 1–2, 3–5, and 6–7). A single porosity was estimated for each of the shallowest and deepest layer groups. Porosity in the middle layer group, where most of the pumping occurs, was estimated by using 37 pilot points (Fig. 2). Continuously varying porosity fields for each layer in which there are pilot points are created by interpolating from the pilot points to the center of each finite-difference cell. Homogeneity in the estimated porosity field is favored by regularization that minimizes differences in porosity values between each pilot point and its nearest neighbors Doherty and Hunt (2010a). In this way parameter estimates at pilot points are constrained by a total of 202 regularization equations. Weights on regularization equations are proportional to their separation distance; the constant of proportionality is calculated by PEST++ using constrained optimization (Doherty and Hunt 2010a). Singular value decomposition is a technique whereby parameters are combined into unique linear combinations called super-parameters. Simulation runs do not need to be done for insensitive super-parameters (Doherty and Hunt 2010a), so this technique can reduce model run time.

The simple inverse model was parameterized into three zones, and one effective porosity value was estimated for each zone. The first zone consists of layers 1 and 2; the second zone consists of layers 3–5, and the third zone consists of layers 6 and 7. Pilot points were not used in the simple model; regularization was used, but differently than in the complex model. In the simple model, a preferred regularization value of porosity of 0.10 was assigned to each zone. As in the complex model, the constant of proportionality for these preferred values was determined by using constrained optimization (Doherty and Hunt 2010a).

Weights specified for all observations were calculated with a coefficient of variation of 0.1 and a threshold of either 10 tritium units (TU) for concentration or 10 years for residence time as

$$ W= \min \left\{\begin{array}{c}\hfill {\left(0.1\times O\right)}^{-1}\hfill \\ {}\hfill 10.\hfill \end{array}\right. $$
(1)

where W is the inverse of the standard deviation and O is the observed value. Using the coefficient of variation (standard deviation divided by the mean) allows larger values to have smaller weights, thereby ensuring that small observations have an effect on the inverse simulation (Hill and Tiedeman 2007). The threshold value is used so that very small concentrations, which might have relatively small standard deviations, do not dominate the solution (Hill and Tiedeman 2007). Although weights on tritium concentration and the apparent age that is derived from tritium concentration are correlated, these correlations are not expected to strongly affect the interpretation of these results (Tiedeman and Green 2013).

Results and discussion

Although many improvements such as revisiting the weighting scheme based on an analysis of influential observations, could be made to the calibration, the intent here is to perform the calibrations as consistently as possible and evaluate the relative merits of the four model/scenario combinations. The goodness-of-fit of the calibrations cannot be judged by the objective function value alone because of the different numbers of parameters, different calibration data sets, and because PEST++ adjusts weights among observation groups automatically between iterations. With that caveat in mind, calibration improved the objective functions by small amounts—by factors of 0.81 and 0.41 for the complex model with and without apparent age data, respectively, and by factors of 0.73 and 0.71 for the simple model with and without apparent age data. The pattern of the estimated porosity distribution for the complex model with age data (Fig. 4) is similar to that of previous estimates for the east side of the valley (Freethey et al. 1994), in that porosity is lowest near the margin of the valley, increases toward the center of the basin, and is higher beneath reaches of tributary streams such as Big Cottonwood Creek and Little Cottonwood Creek (Fig. 1). Estimated porosities for the complex model without age data remained close to their initial values throughout calibration, ranged from 0.09 to 0.13, and exhibited a pattern similar to the complex model but with a smaller range. Porosities for the simple model with age data were 0.11, 0.09, and 0.09 for the shallow, medium, and deep layers respectively, and was uniformly 0.11 for all layers in the simple model without age data.

Fig. 4
figure 4

Porosity values interpolated from pilot points, Salt Lake Valley, Utah. Circles represent pilot points

Full size image

When evaluating the model fit, it is instructive to compare an observation to the entire break-through curve (BTC) rather than to simulated estimates at the one time of the observation. Small changes in porosity shift the simulated curves left or right and, because of the time-dependence of tritium decay, also changes the height of the peak (see Starn et al. 2014 for a graphic example). Wells A, B, and C were selected for display because they are typical fits, are at increasing distance from the basin edge, and have multiple tritium samples and an associated tritium/helium apparent age (Fig. 5). All the models tend to over-predict age, which means that tritium has had more time to decay, and therefore the models also tend to under-predict tritium concentrations. In order to improve the fit, porosities would have to be more than an order of magnitude smaller than those estimated in this study. In general, all the models predicted a temporal trend similar to the trend in the observations.

Well B (Fig. 5) deserves a more detailed discussion. Well B is a pumped well, and during the time tritium samples were being collected (approximately 2000–2010), pumping was halted such that the source of water to the well switched several times between layers 3 and 4. Possible reasons for the low simulated values are (1) the coarse temporal resolution (yearly stress periods) of the groundwater flow model were insufficient to capture flow conditions at the time of sampling; (2) water containing tritium was sequestered in low-permeability lenses within the aquifer, which would essentially be an unsimulated source of tritium when the well was not pumped (multiple transport domains); and (or) (3) the mixture of water from layers 3 and 4 was inaccurately simulated because of unsimulated heterogeneity in those layers. A separate simulation (not discussed in detail here) with a finite-difference transport simulator, which tended to smear temporal fluctuations, did indeed produce a much better fit at well B.

Fig. 5
figure 5

Tritium and apparent age break-through curves. Blue lines are for the complex model. Orange lines are for the simple model. Dashed lines indicate no apparent age data used in calibration. X indicates tritium sample

Full size image

A qualitative way to examine the information content of the four models/scenarios is through the identifiability (IDENT) and relative parameter error reduction statistics (RPER; Doherty and Hunt 2009; comment by Hill 2010; response by Doherty and Hunt 2010b). IDENT and RPER are similar to the composite-scaled sensitivity (CSS) statistic discussed by Hill and Tiedeman (2007) and are appropriate for use in highly parameterized models. Calculation and evaluation of IDENT and RPER is complex and beyond the scope of this article. Simply put, they vary between 0 and 1: a parameter for which IDENT = 0 cannot be estimated with the current model and data set, and a parameter for which IDENT = 1 means that any error associated with the parameter estimate is associated only with measurement noise and not with model inadequacy. Both statistics are sensitive to a value chosen by the modeler (in reference to Doherty and Hunt 2009, it is the number of singular values assigned to model calibration space); however, a poor choice for this number results in values of the statistics less than 0, which did not occur in this case. The number of singular values assigned to the solution space is a function of the number of unique linear orthogonal combinations of parameter values that reduce model error. RPER is similar to IDENT except that it reduces the effect of measurement noise by comparing post-calibration error to pre-calibration error.

For the complex model, which has 39 parameters, 16 parameters have IDENT and RPER greater than 0.80 when using all data (24 singular values assigned to solution space), but when using tritium concentrations, only five parameters have IDENT greater than 0.80 and two with RPER greater than 0.80 (18 singular values assigned to solution space). On this basis, the complex model has significantly more estimable parameters when using all data than with tritium alone. Although not as meaningful because of the small number of parameters, each simple model had one parameter with IDENT and RPER greater than 0.80.

Multiple factors can hinder the ability of a model to reproduce concentrations (Konikow 2011). In this study, these include possible unsimulated processes such as transfer between multiple transport domains (multi-rate mass transfer), unsimulated variations of aquifer properties within flow-model grid cells, and coarse temporal resolution of pumping. With regard to the last factor, in transient flow, groundwater residence time is itself a function of time and varies in relation to stresses on the aquifer. In addition to these factors, there is the difficulty in model calibration of which limb on the tritium BTC to match. Equal fits can be achieved on rising and falling limbs. Apparent ages could help guide the model toward the appropriate limb. On the other hand, tritium/helium apparent ages are based on a piston-flow model and do not account for mixtures of water of different ages, as can be induced in pumping wells and in transient flow regardless of geologic heterogeneity. Tritium/helium apparent age can be significantly less than mean simulated particle residence time; that the amount is less depends on the ratio and age difference of the mixture, which is unknown. Thus, instead of helping choose a limb of the tritium BTC, apparent age could add disinformation. (It is well known that multiple tracers can sort out the ratio of the mixture, but these were not available in this study.)

Conclusions

Two model parameterizations and two data sets were evaluated to inform future modeling effort. Of the models tested, the complex model with tritium concentrations and tritium/helium apparent ages set performed best. Complex models require more time to construct and run, and can be complex to interpret; thus, a simple model might be preferred if it performs adequately. In this study, tritium break-through curves (BTC) simulated by complex and simple models were very generally similar; therefore, there is value in the simple model. The complex model reveals structure in the estimated porosity field that agrees with that found independently by Freethey et al. (1994), and this fact helps support the credibility. The degree to which this structure is expressed increases when apparent age data are included in model calibration. Apparently the potential for disinformation by including apparent age data was less than the information added by guiding the calibration toward the appropriate limb of the BTC. This conclusion is supported by the complex model having more estimable parameters when apparent age data are included than when they are not. Less improvement was noted when apparent age data were added to the simple model.

Two data sets were evaluated. The reduced data set (only used with the simple model) was constructed from the full data set because of the potential problems of matching a point on two-limbed BTCs and of sensitivity to mixing of water from different layers in response to pumping variations. The contention that model calibration might be improved if only measurements that limited those two complications were used was not realized. Even with carefully chosen data, unsimulated features may render that data less than useful, as seen in the BTC to well B (Fig. 5).

Despite many factors that contribute to shortcomings of both the models and the data, there was useful information obtained from all the models evaluated. Although any particular prediction of tritium concentration may have large errors, overall, the models mimic observed trends in tritium breakthrough. Care must be taken to interpret model calibrations, and multiple data types should be used in model calibration.