Daniel T. Feinstein

Groundwater quality changes in wells and streams lag behind changes to land use due to groundwater travel times. Two contaminant transport methods were compared to assess differences in their simulated travel time distributions (TTDs) to streams and wells in the Wisconsin Central Sands. MODPATH simulates advective groundwater flow with particle tracking, while MT3D simulates age‐mass using a finite difference solution without dispersion to allow for direct comparison of the two methods. MODPATH appropriately simulates groundwater TTDs from the water table to surface discharge but is subject to inaccuracies at weak‐sink well cells due to the flow‐model grid discretization and imprecise location of well discharge within well cells. MT3D better represents weak‐sink well cells since it removes mass in proportion to the prescribed pumping rate, although travel time within well cells is neglected. Conversely, MT3D's treatment of surface water boundary cells is not as accurate as MODPATH because mass should be removed from the water table rather than the full cell volume. MT3D simulations of TTDs can also be confounded by the instantaneous vertical distribution of mass introduced throughout recharge cells instead of at the water table, which initiates mass along deeper flow paths. We evaluated 9 MODPATH and 13 MT3D implementations, generating differences in median travel times of up to 18 years. Both methods have strengths and weaknesses, with MT3D better representing weak‐sink well cell behavior and MODPATH better representing surficial recharge and discharge. The effect of these characteristics on simulated TTDs, along with ideas for ameliorating method weaknesses, is discussed.

Groundwater is important as a drinking-water source and for maintaining base flow in rivers, streams, and lakes. Groundwater quality can be predicted, in part, by its residence time in the subsurface, but the residence-time distribution cannot be measured directly and must be inferred from models. This report compares residence-time distributions from four areas where groundwater flow and travel time were simulated with conventional simulation-inset models (IMs) and with a new automated model-construction method called general simulation models (GSMs). The comparison provides an opportunity to explore controls on travel time and improve the methods used in the creation of GSMs. These models can be useful for three main-use cases: (1) rapid testing of relationships that govern groundwater flow and age, (2) generation of consistent examples for training a machine-learning metamodel, and (3) serving as a starting point for more detailed models.Comparison of the GSMs to IMs indicated a qualified pattern of agreement for residence-time distributions as indicated by the Nash-Sutcliffe efficiency and Spearman’s correlation coefficient. The agreement was best for the median values of the simulated residence times in young fractions of groundwater (defined as the fractions of groundwater in samples less than 65 years old) at the scale of the eight-digit hydrologic-unit code. Generally, the median values of the young fractions in the IMs were correlated with the median values from the GSMs. The relative trends across the four areas also were similar for the other residence-time metrics. The medians of residence-time metrics at finer scales show a fair degree of scatter. The GSM results compared most poorly for median travel times in the older fraction of groundwater (older than 65 years).The GSM approach is intended as a flexible framework for developing models that can be useful individually as screening tools or collectively to support projects in statistical learning. Although one set of GSM algorithms was presented here, the approach can accommodate many types of data and also different categories of prior information. Comparison of GSMs and IMs suggests ways in which the GSMs, while remaining easy to construct and calibrate, can be improved for estimating groundwater travel times. IMs do not yield exact travel times, and matching GSMs to IMs does not guarantee an improvement; however, IMs provide a convenient benchmark against which to explore relations between physical characteristics of watersheds and the distribution of travel times within them.This effort was undertaken as part of the National Water Quality Program of the U.S. Geological Survey to assist in determining the susceptibility of groundwater in glacial aquifers to a variety of natural and anthropogenic contaminants.

To help support remedial efforts at the former Badger Army Ammunition Plant the U.S. Geological Survey built and calibrated a transient groundwater flow model using the Newton Raphson formulation (MODFLOW–NWT) of the U.S. Geological Survey’s modular three-dimensional finite-difference code. The model simulates the groundwater flow system at the site from 1984 to 2020. The former Badger Army Ammunition Plant is a 7,275-acre site in Sauk County, Wisconsin. The plant produced smokeless gunpower and solid rocket propellent as munitions components. Peak production periods were during World War II, the Korean War, and the Vietnam War. Subsequent groundwater contamination investigations have found four plumes at the site. A health risk assessment identified at least one contaminant of concern for human health risk present in three of the plumes: the propellant burning ground plume, the deterrent burning ground plume, and the central plume. A cooperative study began between the U.S. Army Environmental Command and U.S. Geological Survey to better understand the groundwater flow system at the former Badger Army Ammunition Plant. Field data, including aquifer tests, streamflow measurements, continuous groundwater elevations, and groundwater gradients with the Wisconsin River were collected and used to inform and calibrate the groundwater flow model. The model was used to assess the variability of the groundwater system over the study period, the components of the groundwater budget, and groundwater flow directions from identified source areas towards the Wisconsin River. Model performance assessment focused on using particle tracking to compare groundwater flowpaths that originate in the contaminant source areas to the observed plume footprints. This focus on plume behavior geometry should help constrain the advective component of a future groundwater transport model of the site.

This study compares the results of two regional steady-state U.S. Geological Survey Modular Three-Dimensional Finite-Difference Ground-Water Flow (MODFLOW) models constructed to quantify the hydrologic changes in the St. Louis River Basin from iron mining on the Mesabi Iron Range in northeastern Minnesota. The U.S. Geological Survey collaborated in this study with bands of the Minnesota Chippewa Tribe, and the Minnesota Pollution Control Agency to inform management decisions about aquatic resources in the St. Louis River Basin. A model constructed and calibrated to represent average 1995–2015 mining conditions produced regional groundwater heads and flows. A pre-mining scenario model was constructed from this mining model but had the land and bedrock surfaces restored to pre-mining topographies and had modeled mining features (mine pits, tailings basins, waste-rock piles, and mining-disturbed areas) eliminated to represent general pre-mining stratigraphy and hydrogeology. Many of the features important to the hydrology of this mining area (like individual mine pits) are difficult to represent in groundwater models and required the use of modeling tools to indirectly account for their effects. The difference between the results of these two models represents mining’s effects on the hydrology in the Mesabi Iron Range area of the St Louis River Basin. The mining and pre-mining regional models also can provide boundary conditions and initial properties for future local or site-specific groundwater-flow models in the area.Total groundwater flow through the mining model is 171 million cubic feet per day. Areal recharge is the largest source of groundwater (78 and 81 percent of total groundwater flow in the mining and pre-mining scenario models, respectively). Seepage from streams and lakes provides another 17 percent of the total groundwater flow through both models. Water leaves aquifers through seepage to streams (discharge as base flow, 43 percent in both models) and areal seepage to the land surface (surface seepage), for example to wetlands (45 and 49 percent, mining and pre-mining scenario models respectively).Comparison of the results from the mining and pre-mining scenario models shows that iron mining has produced measurable hydrologic changes in the St. Louis River Basin, but that most of those changes and the highest magnitude changes occur near the mining features. Flow changes to and from surface-water bodies like streams and wetlands were analyzed in detail because of their importance in sustaining surface waters and aquatic life. Overall, groundwater flow in the mining model was 3.62 million cubic feet per day (2.2 percent) greater than total pre-mining model groundwater flow. This was caused by an increase in recharge from tailings basins and a decrease in discharge from surface seepage. Groundwater discharge to mine pits was the largest change in groundwater flows between the models (a change representing 2.8 percent of total pre-mining model groundwater flow). Net recharge to groundwater from tailings basins (2.4 percent), net decrease in surface seepage from groundwater (2.7 percent), and net increase in seepage to streams (1.0 percent) were all in this same range of total pre-mining model groundwater flow. Groundwater lost through mine-pit withdrawals was nearly offset by groundwater gained through recharge from tailings basins. However, because losses and gains occurred in different areas, the effect of mining can have more substantial effects on local areas than the model-wide averages represent.

Future climate forecasts suggest atmospheric warming, with expected effects on aquatic systems (e.g., cold-water fisheries). Here we apply a recently published and computationally efficient approach for simulating unsaturated/saturated heat transport with coupled flow (MODFLOW) and transport (MT3D-USGS) models via a synthetic three-dimensional (3D) representation of a temperate watershed. Key aspects needed for realistic representation at the watershed-scale include climate drivers, a layering scheme, consideration of surface-water groundwater interactions, and evaluation of transport parameters influencing heat flux. The unsaturated zone (UZ), which is typically neglected in heat transport simulations, is a primary focus of the analysis. Results from three model versions are compared—one that neglects UZ heat-transport processes and two that simulate heat transport through a (1) moderately-thick UZ and (2) a UZ of approximately double thickness. The watershed heat transport is evaluated in terms of temperature patterns and trends in the UZ, at the water table, below the water table (in the groundwater system), and along a stream network. Major findings are: (1) Climate forcing is the product of infiltration temperatures and infiltration rates; they combine into a single heat inflow forcing function. (2) The UZ acts as a low-pass filter on heat pulses migrating downward, markedly dampening the warming recharge signal. (3) The effect of warming on the watershed is also buffered by the mixing of temperatures at discharge points where shallow and deep flow converge. (4) The lateral extent of the riparian zone, defined as where the water table is near land surface (<1 m), plays an important role in determining the short-term dynamics of the stream baseflow response to heat forcing. Runoff generated from riparian areas is particularly important in periods when rejected infiltration during warm and wet periods generates extra runoff from low-lying areas to surface water.

Changes in climate and land use will alter groundwater heat transport dynamics in the future. These changes will in turn affect watershed processes (e.g., nutrient cycling) as well as watershed characteristics (e.g., distribution and persistence of cold‐water habitat). Thus, groundwater flow and heat transport models at watershed scales that can characterize and quantify thermal impacts of surface temperature change on groundwater system temperatures may support forecasting changes to groundwater‐linked ecosystems in riparian zones, streams, and lakes. Including unsaturated zone processes has previously been shown to be important for properly determining the timing and magnitude of groundwater recharge (Hunt et al. 2008). Similarly, heat transport dynamics in the saturated‐zone, as well as connected surface‐water systems, can be appreciably influenced by unsaturated‐zone processes; in this way the unsaturated zone forms an inextricable link between the land surface where change occurs and the groundwater system that transmits that change. This paper presents new capabilities for the existing MT3D‐USGS transport simulator by adding functionality for simulating heat transport through the unsaturated zone. New simulation capabilities are verified through comparison of simulation results with those of the variably saturated heat transport simulator VS2DH under steady and transient conditions for both water and heat flow. The new capabilities are assessed using a number of conceptualizations and include evaluations of convective and conductive heat flow. These additional capabilities increase the utility for applied watershed‐scale simulations, which in turn may facilitate more realistic characterizations of temperature change on thermally sensitive ecosystems, such as stream habitat.

Although there is widespread agreement that future climates tend toward warming, the response of aquatic ecosystems to that warming is not well understood. This work, a continuation of companion research, explores the role of distinct watershed pathways in lagging and dampening climate-change signals. It subjects a synthetic flow and transport model to a 30-year warming signal based on climate projections, quantifying the heat breakthrough on a monthly time step along connected pathways. The system corresponds to a temperate watershed roughly 27 km on a side and consists of (a) land-surface processes of overland flow, (b) infiltration through an unsaturated zone (UZ) above an unconfined sandy aquifer overlying impermeable bedrock, and (c) groundwater flow along shallow and deep pathlines that converge as discharge to a surface-water network. Numerical simulations show that about 40% of the warming applied to watershed infiltration arrives at the water table and that the UZ stores a large fraction of the upward-trending heat signal. Additionally, once groundwater reaches the surface-water network after traveling through the saturated zone, only about 10% of the original warm-up signal is returned to streams by discharge. However, increases in the simulated streamflow temperatures are of similar magnitude to increases at the water table, due to the addition of heat by storm runoff, which bypasses UZ and groundwater storage and counteracts subsurface dampening. The synthetic modeling method and tentative findings reported here provide a potential workflow for real-world applications of climate-change modeling at the full watershed scale.

The Partridge River Basin (PRB) covers 156 square miles in northeastern Minnesota with headwaters in the Mesabi Iron Range. The basin is characterized by extensive wetlands, lakes, and streams in poorly drained and often thin glacial material overlying Proterozoic bedrock. To better understand the interaction between these extensive surface water features and the groundwater system, a three-dimensional, steady-state, groundwater-flow model of the PRB was developed by the U.S. Geological Survey in cooperation with the Great Lakes Indian Fish & Wildlife Commission using the finite-difference computer code MODFLOW-NWT. The model simulates steady-state base flow in streams and groundwater interactions using the streamflow routing (SFR2) package. Existing mining features including tailings basins, stockpiles, pumped mine pits, and flooded mine pits were simulated using either high hydraulic conductivity zones or the drain (DRN) package. The unsaturated zone flow (UZF) package was used to better represent the groundwater system in areas with a high water table and for wetlands often associated with such areas. UZF typically is used to represent unsaturated zone processes but also can simulate the rejection of recharge and groundwater discharge to the land surface when the water table is near land surface. The steady-state model used data from the 2011 to 2013 period when 2011 high-resolution land surface (light detecting and ranging [lidar]) data were available that reflected land-surface and water elevations from mining activity in the basin. The parameter-estimation software suite PEST_HP was used to obtain a best fit of the modeled to measured groundwater levels, streamflow, pit inflow rates, and mapped peat deposits. The PEST calibration used the target residuals from two models with the same model parameters and targets from two separate periods: (1) a 1995–2015 calibration model, which provided a larger number of calibration targets, and (2) a 2011–2013 mining conditions model, which included calibration targets that reflected conditions consistent with the modeled mine-workings topography.Calibration of the PRB model resulted in ranges of glacial horizontal hydraulic conductivity parameters that generally agreed with literature values and other models of the region. Horizontal hydraulic conductivity of the bedrock was higher in the upper bedrock layers where numerous and continuous fractures have been observed and lower in the deeper bedrock layers. Average basin-wide calibrated infiltration was 5.3 inches per year. An average of 4.6 inches per year of infiltration crosses the water table and becomes recharge and 0.7 inch per year is rejected by UZF due to saturated conditions at the land surface. Simulated groundwater runoff (the sum of rejected recharge and groundwater seepage to the land surface) can either be routed to streams or removed from the model as evapotranspiration. The calibrated model indicates relatively shallow groundwater-flow paths dominating and approximately 50 percent of the stream base flow coming from groundwater runoff.The 2011–2013 mining conditions model was then used to develop five model scenarios simulating the response of the groundwater and surface-water system to potential hydrologic stress. The purpose of these mine pit scenarios is to present a possible workflow to quantify a model’s uncertainty for a given model forecast and serve as a possible guide for initial data collection that may improve a future model’s ability to make such a forecast. The scenarios included one scenario with the currently existing Peter Mitchell pit at final buildout and flooded to an elevation of 1,500 feet, and four scenarios with a hypothetical, new mine pit plus the flooded Peter Mitchell at final buildout. The five model scenarios were used to forecast streamflow at six locations in the PRB, pit inflow rates for the new mine pits and the flooded Peter Mitchell pit, and the average depth to water in 12 wetlands. A linear uncertainty analysis was performed using information from the PEST calibration and tools in the PyEMU python package to assess model uncertainty propagation to the model forecasts. Streamflows generally were reduced with future mining and the greatest streamflow reductions occurred from the flooded Peter Mitchell Pit, probably due to its large size. Average depth to groundwater in wetlands was most affected the closer the wetland was to a new mine pit.Linear uncertainty methods were also used to evaluate data worth, which is the ability for potential new groundwater elevation observations to reduce the uncertainty in scenario forecasts. Data worth was performed for a grid of new hydraulic head observations. Overall, areas with nonnegligible data worth generally corresponded to wetland areas with no groundwater seepage to land surface from UZF. These model behaviors indicated that the land-surface boundary condition simulated by the UZF package was pinning the groundwater elevations to the land surface in areas with groundwater seepage (33 percent of the 2011–2013 base conditions model) such that the sensitivity to new observations in these areas was minimal. Therefore, representing wetlands as boundary conditions minimized the usefulness of data worth calculations because wetland areas were present over a large part of the model domain.Probabilistic capture zones were estimated for each of the mines in the model scenarios. A capture zone represents the area contributing recharge to a model feature, like a well or a mine pit, and can be calculated by forward tracking particles from the water table. By using Monte Carlo techniques, it is possible to generate estimated capture zones that include the probability of recharge capture given the uncertainty present in the model. Monte Carlo techniques use randomly generated model parameter sets sampled from a plausible parameter range to create many possible realizations....

Groundwater age is an important indicator of groundwater susceptibility to anthropogenic contamination and a key input to statistical models for forecasting water quality. Numerical models can provide estimates of groundwater age, enabling interpretation of measured age tracers. However, to extend to national‐scale groundwater systems where numerical models are not routinely available, a more efficient metamodeling approach can provide a less precise but widely applicable estimate of groundwater age, trained to make forecasts based on predictor variables that can be measured independent of numerical models. We trained gradient‐boosted regression tree statistical metamodels to MODFLOW/MODPATH‐derived groundwater age estimates in five inset models in the Lake Michigan Basin, USA. Using high‐throughput computing, we explored an exhaustive range of tuning parameters and tested metamodels through cross validation, a 20% holdout, and a round robin approach among the five inset models withholding each inset model from training and testing on the held‐out inset model. Forecast skill—measured by Nash Sutcliffe efficiency—was high for age‐related responses in the 20% hold‐out case (ranging from 0.73 to 0.84). The round robin analysis provided the opportunity to explore extending to unmodeled areas and a greater range of skill indicated the need to evaluate when it is appropriate to apply a metamodel from one region to another. We further explored the ramifications of metamodel simplification achieved through removing predictor variables based on their estimated importance. We found that similar metamodel performance was achievable with a fraction of the candidate set of predictor variables with well construction variables being most important.