Purpose and Scope

This study was undertaken as part of the National Water Quality Program of the U.S. Geological Survey to determine the susceptibility of groundwater in glacial aquifers to a variety of natural and anthropogenic contaminants (Eberts and others, 2013). The purpose of this report is to compare differences between two different methods (IMs and GSMs) that were used in travel-time calculations for glacial aquifers. This report compares the input and output from four areas where groundwater flow and travel time were simulated with IMs (Feinstein and others, 2018; Fienen and others, 2018) and GSMs (Starn and Belitz, 2018). The overlap of these modeled areas enables the comparison of the GSMs to the IMs and provides opportunities to explore controls on travel time and to improve the GSM method for glacial aquifers.

Significantly, neither of these sets of models was calibrated to tracer data, although the GSMs were shown to be consistent with tritium concentrations measured in samples from wells (Starn and Belitz, 2018). The comparison between GSMs and IMs can highlight how systematic differences in model structure affect the calculation of residence-time distributions (RTDs) but may not provide insight into which modeling approach better depicts reality.

The models in this study encompass four watershed areas near Lake Michigan in the States of Wisconsin, Illinois, and Michigan (fig. 1). The models and supporting data are available from Kauffman and others (2023). The watershed boundaries were derived from the Watershed Boundary Dataset (WBD) component of the National Hydrologic Dataset Plus (NHDPlus; U.S. Geological Survey, 2016). The WBD delineates watersheds and assigns them numeric hydrologic-unit codes (HUC). HUCs are hierarchical, with the leading digits indicating larger basins and trailing digits indicating subbasins. The four watershed areas represent four eight-digit HUC (HUC8) watersheds averaging 1,663 square miles (mi²) in area. Within these 4 HUC8 watersheds are 40 ten-digit HUC (HUC10) watersheds averaging 166 mi² in area and 230 twelve-digit HUC (HUC12) watersheds averaging 29 mi² in area (fig. 1; table 1). These models represent a large area that contains diverse glacial deposits in the Lake Michigan Basin (Feinstein and others, 2010), including

Background

Groundwater travel time (sometimes referred to as groundwater age) is defined as the elapsed time for a particle of water to travel from its point of entry into the saturated part of the aquifer, typically at the water table, to a point of interest within the aquifer. Groundwater that discharges to wells or surface-water bodies includes a mixture (or distribution) of different travel times because flow paths converge in these locations. This study is concerned primarily with the RTD of the groundwater in the aquifer, which represents a static snapshot of the ages of the groundwater in an aquifer. The related travel-time distribution (TTD) represents the ages of water exiting an aquifer, such as into a stream or well. TTDs reported for catchments often are much younger than their RTDs (Berghuijs and Kirchner, 2017). The difference is the result of the probability that catchments will drain younger water over older water in storage; however, the two are related through reservoir theory (Bolin and Rodhe, 1973; Eriksson, 1971).

Numerical simulation models of RTDs involve solving the master age-density equation of Ginn and others (2009) either directly on a grid or by tracking the locations of many particles moving through a simulated flow field. Recent research highlights the need for understanding how RTDs are related to landscape (Basu and others, 2012; Condon and Maxwell, 2015; Kirchner and others, 2001; Leray and others, 2016; Maxwell and others, 2016; McDonnell and others, 2010; McGuire and McDonnell, 2006; Remondi and others, 2019; Sprenger and others, 2019; Starn and Belitz, 2018). The relation between groundwater flow and hydrogeographic features has long been noted. For example, Carlston (1966) found that groundwater discharge is proportional to drainage density in streams in the northeastern United States, including nonglaciated areas. Groundwater gradients and flow-path lengths in heterogeneous aquifers are affected by drainage density (Broers, 2004; Goderniaux and others, 2013), and drainage density and groundwater discharge are partly a function of local and regional topography (Tóth, 1970; Winter, 1999, 2007). The position of the surface that separates local flow discharging to nearby streams from underflow to more distant streams is a function of drainage density and recharge (Goderniaux and others, 2013). The presence of local and regional multiscale flow systems and their associated stagnation points leads to fractal scaling of RTDs (Cardenas and Jiang, 2010; Jiang and others, 2012; Kollet and Maxwell, 2008). Aquifer geometry as well as spatial variation of recharge have been found to result in complex RTD shapes (Leray and others, 2016).

MODPATH Inputs for Simulating Groundwater Travel Time

MODPATH was used to calculate groundwater travel-time distributions for both GSMs and IMs by tracking particles backwards from their starting locations to their points of recharge. Two million particles were started in each model and distributed among model cells in proportion to the flux of water moving through the cell. In this way, the RTDs were weighted by the amount of flow associated with each particle (as opposed to weighting by the volume associated with each particle). This differentiation from a volume-based RTD, where each particle represents an equal volume of water in the aquifer, is important. The reason for a flow-based RTD in the aquifer was that it should represent the travel time of discharge to wells or streams at any point within the aquifer. Particles within each cell were started at random three-dimensional locations. TTDs for streams were calculated in a similar fashion, except that particles were started only in discharge cells and distributed randomly on the top cell faces. Particles in the IMs were started only within the GSM boundary corresponding to a HUC8 basin, although the IMs have active areas outside the GSM boundaries.

MODPATH requires additional input beyond what is required for GSMs or IMs: these two model-type inputs were made as consistent as possible. Sources and sinks of water (recharge, drains, streams, and lakes) were placed on the upper faces of model cells. A correction to travel time in surface-water cells was applied (Abrams, 2013): if both sources and sinks were in the same model cell, the net value of flux was used.

For the purpose of comparison, the effective porosity was set uniformly to 0.20. Although this is a reasonable value for effective porosity overall, it is probably high for the bedrock part of the system and likely low in parts of the glacial sediments.

The SFR (streamflow-routing package in MODFLOW) was used to represent streams in the IMs, which allows water to leak back into the groundwater system. Generally, the flow paths followed by the stream leakage are short, and any water lost to the aquifer generally discharges back into a downgradient stream reach within a short distance. Because most of the MODFLOW-simulated stream leakage into the groundwater system resulted from the way that streams were discretized and stream altitudes were assigned, these flow paths were not used in the computation of the RTD and corresponding travel-time metrics.

Model-Input Characteristics Used for Comparison

A set of 26 model-input characteristics (including the mean and standard deviation) was calculated for each model type and for each HUC8 watershed and their HUC12 subwatersheds (fig. 4; table 3). These factors can be divided into four groups related to (1) the volumetric amount of water available (in the form of recharge for GSMs and of recharge and stream losses for IMs) compared to the volume of the saturated aquifer, (2) the surface-water network, (3) the ease of lateral flow, and (4) the ease of vertical flow. The mean and standard deviation were calculated in log space for factors related to transmissivity, hydraulic conductivity, and vertical anisotropy (the ratio of horizontal to vertical hydraulic conductivity), so that the expected value would be the geometric mean, and the standard deviation would be the log-based standard deviation defined in National Institute of Standards and Technology (2003).

Some factors were calculated for the glacial part of the saturated-flow system, some for the bedrock, and some for the total-flow system. In the GSMs, bedrock refers to the single bedrock layer and any parts of overlying layers with the same bedrock hydraulic conductivity. For the IMs, bedrock refers to the multiple layers assigned to bedrock units above a depth where the ratio of thickness to the vertical hydraulic conductivity of a geologic unit is greater than 50,000 days. This depth, corresponding to shallow bedrock, varies on a cell-by-cell basis in the IMs.

After preliminary analysis, 17 of the 26 potential explanatory factors were shown to influence results sufficiently enough to warrant more extended correlation analysis. Because of the importance of the zonal structure in the GSMs, several factors were calculated for these zones in the corresponding areas in the IMs, including recharge, horizontal glacial hydraulic conductivity, and shallow-bedrock hydraulic conductivity.

The differences and similarities between the two model types were summarized by using the ratios of model-input parameters (GSM:IM; fig. 5). A value of 1 indicates that the parameter is the same in both model types, and a departure from 1 in either direction indicates less agreement between the characteristics of the two model types.

Characteristics related to saturated volume (fig. 5A) showed agreement for recharge and glacial tau (the ratio of saturated glacial volume to recharge rate, generally equivalent to the time needed to flush the glacial material and thus directly tied to residence time). Differences related to the surface-water network (fig. 5B) were attributable in large measure to the cell-size difference in the models; more GSM cells contained surface water. The small but consistent gradient difference in the computed GSM:IM ratio between the water table at locations within the HUC8 and the nearest downgradient stream also could have been a function of grid size and could cause slower velocities in the GSMs. Characteristics related to ease of flow (fig. 5C–D) differed between the two model types, with notable differences in horizontal hydraulic conductivity of bedrock (K_b) and three characteristics related to vertical flow.

Groundwater RTD Metrics and Comparison Metrics

Four RTD metrics were computed on the basis of particle-tracking results. A threshold of 65 years was used to split young and old water to correspond roughly with bomb testing, which released large amounts of tritium into the atmosphere and thus allowed tritium concentrations to become a travel-time indicator in the groundwater system. Also, this threshold (water with a travel time less than 65 years) was considered a general indicator of potential anthropogenic contamination. In addition to the RTD, the following travel-time metrics were computed:

The agreement between model types was evaluated by using the Nash-Sutcliffe efficiency coefficient (NSE). NSE is the proportion of variance in one model that is predicted by the other model and can have values from −∞ to 1.0. An NSE of 0 indicates that the mean of observations is the best predictor of the observations (as opposed to a more complicated model). Values less than zero indicate that predictions are worse than those predicted by the mean, and that values close to 1 indicate greater correlation among RTDs. The equation for the NSE is

Comparison of RTDs in Aquifers and TTDs in Streams by HUC8 Basin

The NSEs for simulated RTDs in the aquifers between IMs and GSMs ranged from 0.67 in the BOARD model to 0.91 in the UPFOX model (fig. 6). RTDs from the GSMs indicated generally older water than from the IMs for all four basins. The shapes of the curves generally are similar within a given model, the most obvious exceptions being in BOARD, where the IM had a broader and possibly bimodal peak (fig. 6B). RTDs from the GSMs were less symmetrical, or more skewed, than those from the IMs and had a greater density of travel times less than the peak travel time.

The NSEs for simulated TTDs for streams were less than for aquifer RTDs for models of the MANI and UPFOX basins and greater for models of the KALA and BOARD basins (fig. 7). GSMs generally simulated older water than the IMs, as with the aquifer RTDs. Three of the stream RTDs for both model types were generally unimodal with a well-defined mode (fig. 7A–C); the UPFOX GSM RTD suggests two modes (fig. 7D). Stream TTDs from the GSM, like the aquifer RTDs, are less symmetrical than those from the IMs and have a greater density of travel times less than the peak travel time, especially in the MANI and UPFOX models. The IM TTD peaks were earlier than the GSM peaks in all cases due to a faster rising limb, indicating that the IMs tended to produce large groundwater velocities and small flow paths. The 65-year travel-time cutoff (between young and old fractions) crosses the falling limb for three of four IMs, whereas it is near the peak for all GSMs (fig. 7).

NSE values between aquifer RTDs and stream TTDs ranged from 0.83 to 0.96 for IMs, indicating that the age distributions of water in the aquifer and in discharge were similar, and the effect of water flowing through bedrock to streams was minimal (fig. 8A, C, E, and G). NSE values for aquifer RTDs and stream RTDs were lower for GSMs, ranging from 0.57 to 0.85 and thus suggesting a more important role for bedrock (fig. 8B and H).

Comparison of RTD Metrics Aggregated Over HUC8 and HUC12 Watersheds

The median age for groundwater in both the HUC8 (fig. 9A) and HUC12 (table 4) basins typically was higher for the GSMs than the IMs, particularly in the case of BOARD; they are almost equal in the case of MANI. The mean young fraction (less than 65 years) was always greater for the IMs (fig. 9B), a finding also made evident by their RTDs (fig. 6). The median age was generally younger for the IMs (fig. 8A and C) than for the GSMs, except for the MANI basin (fig. 6C). The median ages of the “old” fraction of water (greater than 65 years) generally were the same for both types of models (fig. 9D).

Although the GSMs and IMs used in this study were constructed at the HUC8 watershed scale, age metrics were also calculated for HUC12 watersheds in this study to determine how they would compare to the age metrics of the IMs at smaller scales. The RTD metrics (fig. 10) showed substantial scatter and a different correlation between the GSM and IM for each modeled area. The largest differences were in the metrics for the BOARD basin, where the GSM predicted older median ages of young water in all and a smaller fraction of young water in all but one of the HUC12 basins. The correlations were positive in all cases, but the slope of the trendline was always less than one, possibly pointing to some bias in the GSM results compared to the results from IMs. The slope of the trendline in the predicted fraction of old water was close to 0 in the MANI basin (fig. 10D); this result could indicate that the IM produces a greater variability in the age of the old water than does the GSM—results that are consistent with differences observed in the RTDs.

Correlation of RTD Metrics With Model-Input Characteristics

Correlation between model-input characteristics and RTD metrics generally showed that characteristics strongly correlated with RTD metrics (either positively or negatively) for the IMs also were correlated strongly with RTD metrics for the GSMs (fig. 11). The median age of glacial sediments across HUC12 watersheds for both model types was negatively correlated with recharge rate and the magnitude of glacial hydraulic conductivity: the higher those parameters, the lower the median age. The median age was positively correlated with the vertical resistance of saturated glacial material and the ratio of the volume of saturated glacial material to the recharge rate (glacial tau, variable name GLAC_tau) (fig. 11A).

Seven characteristics had correlations that exceeded 0.25 for the fraction of young water. Total saturated thickness, glacial saturated thickness, vertical resistance of glacial material, and glacial tau were negatively correlated, with glacial tau having the highest correlation. Recharge and glacial hydraulic conductivity were less strong but positively correlated (fig. 11B).

The median age of young water was strongly and positively correlated to glacial tau for both model types (fig. 11C). Recharge and glacial hydraulic conductivity were negatively correlated to this metric, whereas glacial vertical resistance was positively correlated. The pattern is less consistent for average glacial vertical anisotropy and shallow bedrock resistance.

The median age of the old fraction of water in both model types was most strongly and positively correlated to glacial recharge and glacial and total saturated thickness. The direction of correlation was opposite for many other factors. The negative correlations in the cases of the IMs may indicate that high bedrock resistance tends to keep many of the long flow paths in the glacial part of the system yielding younger water for residence times greater than 65 years, whereas the positive correlation for the GSMs may be affected by the higher proportion of flow which enters the bedrock (fig. 11D), which, having entered, tends to be older.

Spearman’s ρ for all the potential factors showed that the two model types agreed better with respect to the relative explanatory power of each factor for the median age and young- fraction metrics (fig. 12A and B) than they did for the median age of young fraction and median age of old fraction (fig. 12C and D). The model types disagreed the most for the median age of old fraction, with almost no correlation between model inputs and that metric (fig. 12D). This result complements the finding that, between the two types of models, the simulated median age and young-fraction metrics tend to correlate better than do the median ages of young and old fractions for both RTDs (table 4).