Badger Army Ammunition Plant History

The former BAAP is a 7,275-acre site in Sauk County, Wisconsin (fig. 1). The BAAP was constructed in 1942 and 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; during these times, the plant employed as many as 7,500 workers. Between peak production periods, the site was maintained with a standby status until a 1997 decision to close the plant permanently (SpecPro Professional Services, LLC, 2019).

During plant operations, waste disposal included burning and burial that, together with wastewater conveyance through open ditches and leaky sewer pipes, have affected the soil and groundwater quality at the former BAAP. Site investigations and restoration began in 1977, and most buildings were removed between 2002 and 2012. Groundwater contamination investigations began in 1980; various remediation actions have occurred since then, and monitoring is ongoing (SpecPro Professional Services, LLC, 2019). Groundwater contamination was discovered in the following plumes (fig. 1):

The PBG plume (fig. 1) was first detected in 1982 and has migrated off site from the southwestern corner of the former BAAP site. The contamination source for this plume includes a landfill, several waste pits, and open burning areas (fig. 1; landfill 1, propellant burning ground waste pits, 1949 pit, and the racetrack area). Past remedial activities at the source areas include capping of the waste pits with an impermeable membrane, excavating and disposing contaminated soil, operating a soil vapor extraction system, and operating a bioremediation treatment system. Despite the remedial activities at the source areas, some soil contamination likely remains beneath the caps. In addition to source area remediation, an interim remedial measure (IRM) pump and treat system began operation in 1990 and was expanded with the modified interim remedial measure (MIRM) system in 1996; and additional extraction wells were added in 2005. Together, these systems included varying combinations of source control wells and extraction and boundary control wells. By 2015, the contaminant mass being removed by the pump and treat system had diminished and the system was shut down (SpecPro Professional Services, LLC, 2019).

The DBG plume (fig. 1) has migrated off site from the northeastern corner of the former BAAP site. DBG contamination source areas are a landfill and burn pits where surplus propellant and other site waste were burned during plant operation (fig. 1; deterrent burning ground and landfill 3). Past remedial activities at the DBG source areas include capping of the waste pits with an impermeable membrane, excavating and disposing of some contaminated soil, and operating a bioremediation system. Soil contamination likely remains at the DBG source area beneath the caps and is thought to be more than 20 feet (ft) above the water table (SpecPro Professional Services, LLC, 2019).

The central plume (fig. 1) was discovered in 2004. A defined source area for the central plume is unknown, but the plume is thought to have originated from the north-central part of the former BAAP near the rocket paste production area (fig. 1), which had open ditches that could have leached production wash water through the unsaturated zone to the water table. Past remedial activities include general soil excavation activities near production buildings and along ditches and drainages from the rocket paste and nitrocellulose production areas. Since remediation, it is assumed but not confirmed that the central plume source was eliminated (SpecPro Professional Services, LLC, 2019).

An additional small plume was discovered in 2007 in the northwest corner of the former BAAP. This NC plume (fig. 1) is thought to have originated broadly across the 800-acre smokeless powder and nitrocellulose production areas. Soil testing in this NC production area found DNT contamination likely associated with leaky sewer pipes carrying contaminated wash water and with slab foundation cracks that may have allowed floor wash water to migrate down into the unsaturated zone and then to the water table. Remedial activities have included the removal of the sewer pipes, building slabs, and associated soils in areas with confirmed DNT (SpecPro Professional Services, LLC, 2019).

Contaminants associated with each plume are discussed extensively in SpecPro Professional Services, LLC (2019), a report prepared for the USAEC. The plume contaminants of concern generally include some combination of the six isomers of DNT, total DNT, and various volatile organic carbon compounds. As part of SpecPro Professional Services, LLC (2019), a groundwater human health risk assessment was completed for each contaminant of concern for onsite and offsite groundwater. Based on this risk assessment, the study identified the following human health risk contaminants of concern for each plume:

SpecPro Professional Services, LLC (2019) identified potential remedial alternatives for each plume. Based on the remedial alternatives considered, the volatile organic carbon compounds are expected to continue to reduce to acceptable levels through monitored natural attenuation. Active remedial options are being considered by the USAEC for 2,6-DNT in the PBG and central plumes and for total DNT in the DBG plume.

Previous Investigations

Previous groundwater modeling investigations in this area include a GFLOW analytic element groundwater flow model of Sauk County (Gotkowitz and others, 2005), a MODFLOW–NWT model of Columbia County (Gotkowitz and others, 2021), and a MODFLOW model and MT3DMS transport model of the former BAAP site (Zeiler, 2002). The former BAAP site (fig. 1) is at the edges of the Sauk County GFLOW model and the Columbia County MODFLOW–NWT model. Although these county-scale models provide important information about the regional groundwater system, they lack the detail and temporal variability to simulate groundwater flow at the former BAAP at the scale necessary for transport modeling. The Zeiler (2002) MODFLOW and MT3DMS models focused on the transport of carbon tetrachloride during the plant production period through 1998 but did not simulate more-recent groundwater flow conditions or have the large dataset of groundwater elevations collected in the site monitoring wells in the 2000s.

Previous investigations have characterized parts of the groundwater flow system and contamination at the former BAAP. The past remedial activities to reduce the extent and concentration of the groundwater contamination plumes were supported by various studies for the design and implementation of these actions. Notable site studies are summarized in detail in SpecPro Professional Services, LLC (2019). Well construction and groundwater elevation data from site wells are summarized in a Microsoft Access database discussed in appendix 1 (Joel Janssen, Project Manager, SpecPro Professional Services, LLC, written commun., December 19, 2019 [These data are not publicly available. Contact SpecPro Professional Services, LLC, for further information.]). Other previous work, as related to developing the site conceptual model and model layering is discussed in the “Site Description and Hydrologic Setting” section.

Study Overview and Objectives

The USAEC started this study with the USGS to better understand the groundwater flow system at the former BAAP. Specifically, the goals of this study were to (1) create a tool to communicate information about groundwater at the site to local stakeholders, (2) inform the advective component of a future groundwater transport model, and (3) ultimately support future remedial design and implementation at the site. The overall study has two parts. First, create and calibrate a transient groundwater flow model, as summarized in this report. Second, use the flow model to simulate contaminant transport at the site under the USAEC-chosen remedial strategy for each of the plumes. The groundwater flow model was informed by and calibrated to collected field data, including aquifer tests, streamflow measurements, continuous groundwater elevation data, and groundwater gradients with the Wisconsin River. The model simulates the groundwater flow system from 1984 to 2020, a period that begins about when groundwater contamination was discovered on the site and calibration data collection started. The model focuses on the PBG, DBG, and central plumes where active remediation is being considered for DNT contamination. The model was used to assess the variability of the groundwater system over time and the basic components of the groundwater budget.

Purpose and Scope

This report presents the results of an investigation of the groundwater flow system near the former BAAP, including the collection and interpretation of field data and development of a finite difference groundwater flow model. The report includes (1) a description of the hydrologic setting and geology for the groundwater flow model domain (study area); (2) a summary of the field data collection methods; (3) presentation and discussion of field data results; (4) a description of the groundwater modeling approach, including model construction and parameter estimation; (5) a summary of the model results; and (6) a discussion of model limitations and assumptions made during this study. The main body of the text is intended to summarize the study and highlight the key findings. Additional details are provided in the appendixes. All USGS data and models for this study are available in a USGS database. The slug test data (Corson-Dosch and Haserodt, 2023), the Soil-Water-Balance model (Nielsen, 2023), and the MODFLOW–NWT model (Reeves and Corson-Dosch, 2023) are in data releases. Streamflow data and well water level data collected for this study are in the USGS National Water Information System (NWIS; USGS, 2021).

Bedrock and Surficial Geology

Previous reports on the bedrock and surficial geology of the region informed the geologic model of the former BAAP site and surrounding area (Dalziel and Dott, 1970; Clayton and Attig, 1990; Clayton and Attig, 1997; Brown and others, 2013a,b; Ostrom, 1965; Parsen and others, 2016; Carson and others, 2019; Lundqvist and others, 1993; Krohelski and others, 2000; Gotkowitz and others, 2005). Additional geologic data came from USAEC monitoring well logs (Joel Janssen, Project Manager, SpecPro Professional Services, LLC, written commun., February 1, 2021), Wisconsin Department of Natural Resources county well reports for Dane, Columbia, and Sauk Counties (Wisconsin Department of Natural Resources, 2021), and well logs from Wisconsin Geological and Natural History Survey (WGNHS) online datasets in the area (available through https://wgnhs.wisc.edu/catalog/publication), including the Columbia County groundwater flow model (Gotkowitz and others, 2021), which has the former BAAP site in its northwest corner. The general geology of the site is described in this section of the report, and the conceptualized hydrostratigraphic units for the groundwater flow model are described in the next section. A simplified map of the surficial geology near the site is shown on figure 2. The units in figure 2 are primarily derived from three surficial, glacial, and bedrock geologic maps that cover the study area (Attig and Clayton, 2003; Hooyer and others, 2015; Clayton and Attig, 1997) because no single map covered the surficial units in the study area.

The former BAAP site sits at the edge of the Driftless Area of Wisconsin (Carson and others, 2019), an area with no evidence of glaciers since the deposition of the Cambrian and Ordovician sediments. The bedrock surface under the former BAAP site and the Wisconsin River valley in this area was cut by an earlier river system that was thought to flow away from the Mississippi River (not shown) rather than towards it as the Wisconsin River does today (Carson and others, 2019). The Johnstown Moraine (Lundqvist and others, 1993; Clayton and Attig, 1990), the farthest western edge of the last glaciation, is in the middle of the former BAAP site (fig. 2). Pro-glacial sediments derived from glacial runoff fill the pre-glacial valley to the west of the terminal moraine (Clayton and Attig, 1990). These sediments and the overlying moraine sediments (Lundqvist and others, 1993; Clayton and Attig, 1990) consist of very coarse- to fine-grained grained outwash sediments, and a sandy till cap that composes the Johnstown Moraine. To the east and north, some of the surficial glacial sediments were deposited in meltwater lakes to the west of the moraine and are more clay-rich (Clayton and Attig, 1990) as shown on the geologic map (Clayton and Attig, 1990, plate 1) and in well logs.

The bedrock geology of the region is composed of the Precambrian Baraboo Quartzite in the Baraboo Hills (figs. 2, 3; Dalziel and Dott, 1970; Clayton and Attig, 1990), an erosional remnant of the Canadian Shield rocks that underlie all of Wisconsin. These older rocks are flanked by Cambrian- and Ordovician-aged sedimentary rocks deposited in a shallow sea above the Precambrian surface (Dalziel and Dott, 1970; Clayton and Attig, 1990). These younger strata, as described in Dalziel and Dott (1970) and Clayton and Attig (1990) unless otherwise noted are, from oldest to youngest; the Parfreys Glen Formation (Clayton and Attig, 1990), a quartzite conglomerate and sandstone that laps against the Baraboo Quartzite; the Mount Simon Formation (fig. 3), a sandstone that does not crop out in the area (based on drill cuttings, the unit underlies the shallower bedrock units and is known to be greater than 200 ft thick); the Eau Claire Formation of the Elk Mound Group (hereafter referred to as the Eau Claire Formation), a heterogeneous unit composed of siltstone, shale, dolostone, and fine sandstone that also does not crop out in the study area but is interpreted from drilling; the Wonewoc Formation of the Elk Mound Group (hereafter referred to as the Wonewoc Formation) (figs. 2, 3; Ostrom, 1965), the oldest of the sedimentary units that crops out at land surface in the study area, that is a well sorted sandstone exposed under the bottom of the pre-glacial river valley; the Tunnel City Group, a dolomitic fine sandstone with thin dolomite layers about 100 to 150 ft thick; and the rocks of the Trempealeau Group, which most commonly are exposed on the hills on the east and west side of the Wisconsin River valley. The Trempealeau Group includes dolomite and siltstone of the St. Lawrence Formation and the Jordan Formation sandstone above. The topmost stratigraphic unit of the study area is the Oneota Formation, a dolomite that forms the caprock of many of the hills. Bedrock units that crop out in the study area are shown on the surficial geology map (fig. 2). The glacial and bedrock sequence is shown in a stratigraphic chart in figure 3.

Many of the bedrock units also have a weathered saprolite layer where they have been exposed to the surface and chemical erosion, especially to the west of the terminal moraine. The saprolite zones are not listed specifically in figure 3 but were added to the group of hydrostratigraphic units represented in the model.

Hydrostratigraphic Units

For the former BAAP groundwater model, the unconsolidated and geologic units at the site and the surrounding area were divided up into hydrostratigraphic units based on the lithology, depositional history, and water-bearing characteristics of the units. The bedrock units were generally assigned to hydrostratigraphic units based on the specific formations, but the surficial units were broken down into texture-based hydrofacies representing different sediment sizes.

Texture-based hydrofacies were used to classify the surficial hydrostratigraphic units for two reasons: (1) because well logs, the primary source of information used in the classification, describe the sediments using textural descriptions, and (2) this type of classification provided convenient categories to facilitate the model calibration of hydraulic conductivity for the upper unconsolidated model layers. Seven texture-based hydrofacies were used in the surficial hydrostratigraphic unit classification: gravel, sand and gravel, sand, mixed-coarse, mixed-fine, silt and clay, and clay. The mixed-coarse and mixed-fine facies were used when an interval in the log contained a mix of coarse-grained and fine-grained sediments but either the coarse fraction or fine fraction was dominant. The other facies were used when the log interval description was relatively uniform. Additional detail about the hydrofacies model development is presented in appendix 3.

The hydrofacies across the model domain are dominated by coarse-grained sediments, such as gravel, mixed sand and gravel, and sand, which is a good representation of the outwash sand and gravel deposits that dominate the study area (fig. 2). Because the water table across much of the former BAAP site is 25–50 ft below the bottom of the surficial sandy till of the Johnstown Moraine, the uppermost hydraulically active hydrofacies are the coarser-grained outwash deposits that are found under the surficial till. Fine-grained silt and clay deposits are most common to the north of Lake Wisconsin as well as northwest and east of the former BAAP site. A series of three-dimensional views of the hydrofacies model (fig. 4) illustrates the distribution of the units across the study area. The method used to classify the surficial geology for the groundwater flow model is described in appendix 3.

The hydrostratigraphic units (fig. 5) in the bedrock were simpler to assign because they correspond to the geologic formations. Because of their position in the landscape (forming the hills on either side of the Wisconsin River valley) and the depth to groundwater within them, the Tunnel City Group, Trempealeau Group, and Oneota Formation are lumped together as one unit for the purposes of the groundwater flow model. These form the uppermost bedrock hydrostratigraphic unit called the “upper bedrock” unit. Below the Tunnel City Group, the Wonewoc Formation is the next hydrostratigraphic unit and is present under much of the surficial sediments in the bedrock valley, although it is a relatively thin unit; and the bedrock surface is scoured below the bottom of the Wonewoc in the deepest parts of the bedrock valley below the Wisconsin River. Below the Wonewoc, the Eau Claire Formation is a thick unit above the lowest unit, the sandstone Mount Simon Formation. The Eau Claire Formation contains layers of shale, siltstone, and dolomite that form an effective confining layer between the upper sandstone of the Wonewoc and the lower Mount Simon Formation.

Mini-Piezometers

Six mini-piezometers were installed along the Wisconsin River (fig. 6) to assess the vertical hydraulic gradient between the shallow groundwater system and the Wisconsin River within the study area. The mini-piezometer locations were selected to assess spatial variability in the vertical hydraulic gradients upstream and downstream from the Prairie du Sac hydroelectric dam.

Continuous Groundwater Elevations

Otter Creek Synoptic Streamflow Measurements

Results

Data from the synoptic streamflow survey are available in NWIS (USGS, 2021) using the site numbers (figs. 6, 7). Some sites had fewer measurements because of unfavorable measuring conditions. The highest streamflow values were from the spring 2021 and 2022 baseflow measurements. Streamflow generally increased in the headwaters of Otter Creek until the stream section between sites 05406240 and 05406247, where streamflow was relatively flat or slightly losing. Flows then decreased between sites 05406247 and 05406252, which indicates that Otter Creek is losing water to the groundwater system somewhere along this stretch. Overall, Otter Creek seems to gain water in the upper part of the watershed and then lose flow in at least parts of the middle to lower watershed.

Slug Tests

Slug Test Analysis

Slug test water level data were analyzed with the AQTESOLV program (version 4.5; Duffield, 2007) using the solution method judged to be most appropriate based on well construction (all wells were partially penetrating) and the water level response to the slug (oscillatory or nonoscillatory). The solution methods used included the Kansas Geological Survey method for unconfined or confined settings (Hyder and others, 1994), and the Springer-Gelhar method (Springer and Gelhar, 1991) for unconfined. The Kansas Geological Survey method was used to analyze nonoscillatory data, and the Springer-Gelhar method was used to analyze oscillatory data. Nonoscillatory responses are common during slug tests in wells screened in aquifers with relatively low or medium hydraulic conductivity, whereas oscillatory responses may be observed during slug tests in wells screened in aquifers with high hydraulic conductivity (Duffield, 2007). Well construction information from well construction logs (Joel Janssen, Project Manager, SpecPro Professional Services, LLC, written commun., December 19, 2019) was used by AQTESOLV for both methods of analysis. Example AQTESOLV results from well PBN–1304C are shown in figure 8. The plot shows normalized head (observed water-level displacement divided by the initial displacement) with respect to elapsed time since the start of the slug test and the Kansas Geological Survey model fit to those data.

Discretization

The Badger model has a regular grid with a uniform horizontal cell size of 150 ft by 150 ft. The model area is subdivided into 392 rows and 362 columns with 14 numerical layers. Model cells covering the Baraboo Hills were set to inactive because of limited groundwater flow potential through the Precambrian bedrock. Of the total 1,986,656 model cells, 214,788 are inactive. MODFLOW convertible layers are used so that the model calculates the altitude of the water table while accounting for unconfined conditions.

Model layers 1–10 represent the unconsolidated materials above the bedrock, and layers 11–14 represent the bedrock units. Table 3 provides a summary of the geologic units represented by each model layer. Figure 9 has a map of bedrock surface and cross sections of the model layers across the model domain. Layer elevations were assigned using the following logic:

Stress Periods

The transient Badger model spans the period from October 1, 1984, through September 30, 2020; this period was selected based on the availability of site well data. In a MODFLOW simulation, a stress period is a part of the model during which boundary conditions such as recharge, pumping, or stream and river stages are held steady. The model has a total of 50 stress periods, starting with a steady-state period representing average conditions from 1980 to 1984, followed by 49 transient stress periods (fig. 11). Each of the 49 transient stress periods are subdivided into 5 time steps with a short step at the beginning of the stress period followed by progressively larger time steps for the rest of the period. This subdivision is done to improve the numerical accuracy of the model. Most stress periods represent a single water year (for example: October 1, 1984, through September 30, 1985). However, some of the stress periods are broken down more finely. For example, water year 2005 is broken into two periods to properly represent the response of the system to pumping associated with remedial activities (MIRM system) at the site in the second half of that water year. Additionally, water years 2017 through 2020 are broken into quarters to better capture water table dynamics, including a recent rise in the water table, and because there are enough data in recent years to support calibration to seasonal variations. Better understanding of these recent dynamics may help support remedial system design and interpretation of recent spikes in concentrations at some monitoring locations.

Boundary Conditions

As noted above, the Badger model simulates groundwater flow in unconsolidated sediments and sedimentary bedrock. Hydraulic boundary conditions were applied to represent flows into and out of the model. The bottom boundary of the model is a no-flow boundary representing the contact between Paleozoic sedimentary bedrock and Precambrian crystalline bedrock units. The top boundary is a specified flux representing recharge to the groundwater system. The lateral boundary conditions are a combination of specified heads and specified fluxes extracted from a regional model for Columbia County (Gotkowitz and others, 2021) that are forced using the MODFLOW Well (WEL) package. These lateral boundaries are set far away from the former BAAP site and the identified plumes to minimize their effect on computed flow directions and rates.

Potential recharge to the groundwater system was applied to the surface of the model active area using estimates from the USGS Soil-Water-Balance (SWB) model (Westenbroek and others, 2018) for the area surrounding the former BAAP (Nielsen, 2023). The SWB model uses climate data and readily available spatial datasets of landscape and soil properties to calculate a Thornthwaite-Mather soil water balance on a daily time step and produce an estimate of net infiltration (also referred to as potential groundwater recharge) of water beneath the rooting zone. Potential groundwater recharge grids were estimated by the SWB model on a daily time step and aggregated to total potential recharge for every time step in the model simulation. The SWB model was also used to estimate the runoff from upland watersheds on the Baraboo Hills that seep into the subsurface through losing stream segments. This runoff was applied as recharge from stream segments near the base of the Baraboo Hills using the WEL package. The average annual simulated potential recharge for the model area from 2000 to 2020 is shown in figure 12 to illustrate the general pattern of recharge cross the site. Additional details on the SWB model development and results for the area near the former BAAP are in appendix 2 of this report, and model files are available in a model archive (Nielsen, 2023).

The internal boundary conditions imposed on the Badger model (fig. 13) were represented using the following packages:

Permitted groundwater withdrawals for municipal and irrigation wells in the model area were minimal and distant from the plumes so they were inherited as the steady-state pumping rates from the Columbia County groundwater flow model (Gotkowitz and others, 2021). Transient pumping from the interim remedial measure/MIRM pump-and-treat wells (Joel Janssen, Project Manager, SpecPro Professional Services, LLC, written commun., February 2, 2021) on the former BAAP site was also included. Minor withdrawals from private residential wells are not included because most of the water pumped from these wells is returned to the system through onsite septic, and the rates are typically small enough to not greatly affect groundwater flow directions.

The SFR2 package allows the model to estimate the stream stage and flux between the aquifer and the stream for Otter Creek and other smaller streams. In this model, streamflow represents base-flow conditions reflecting exchanges with the groundwater system and does not include stormflow components. The interaction between the stream and the aquifer is determined by (1) the stage in the stream relative to the groundwater elevation (hydraulic head) in the aquifer and (2) a conductance term based on the streambed dimensions and the vertical hydraulic conductivity (Kv) of the streambed sediment. Losses from the stream to the aquifer are controlled by the streamflow in that reach. The package was set up using the SFRmaker tool (Leaf, 2018) with hydrography information for this area from the National Hydrography Dataset Plus (McKay and others, 2012) and the land surface elevations from the light detection and ranging digital elevation model (USGS, 2019; University of Wisconsin–Madison, 2019a,b).

The RIV package is a head-dependent flux boundary condition based on the hydraulic head in the aquifer relative to a specified river elevation and a conductance term. River elevations were assigned using average water surface elevations above and below the hydroelectric dam with an assumed linear gradient in river surface downstream from the dam (water elevations in dataset provided by Amanda Blank, Site Manager, Hydroelectric & Gas Generation, Alliant Energy, written commun., December 29, 2020 [These data are not publicly available. Contact Alliant Energy for further information.]). The linear gradient was calculated using the average water elevation below the dam (starting at model row 288) and the light detection and ranging elevation at the outlet of Otter Creek, as used in the SFR2 package. The conductance of the riverbed depends on the model cell area, the thickness of the riverbed, and the Kv of the riverbed material. To account for sedimentation above the dam, which is thought to reduce the Kv, the conductance of the riverbed was divided into five zones: three representing the area upstream from the dam and two for the area downstream from the dam (fig. 13). Representing the Wisconsin River with the RIV package is less sophisticated than the SFR2 package but is reasonable because the river stage does not change much during the simulation. Trial simulations demonstrated that using the RIV package with fixed stage values can accurately represent the river as a discharge zone for regional flow as well as localized flow affects from the dam.

Aquifer Properties

Aquifer properties, including hydraulic conductivity (Kh and Kv), specific storage (Ss), and specific yield (Sy), were defined by the hydrostratigraphic zones from the unconsolidated hydrofacies for the upper 10 model layers and the mapped bedrock units for the lower 4 layers (see app. 3). Zone categories for the unconsolidated included gravel, sand and gravel, sand, mixed coarse, mixed fine, silt and clay, and clay. Zone categories for the bedrock layers included the upper bedrock unit, the Eau Claire Formation (including separate weathered zone), the Wonewoc Formation (including separate weathered zone), the Mount Simon Formation, and a zone representing the weathered upper bedrock and the Baraboo transition zone for the thin wedge of unconsolidated and bedrock material along the southern base of the Baraboo Hills. The Eau Claire and Wonewoc weathered units represent saprolitic zones near the top surface of the bedrock where weathering can give rise to enhanced hydraulic conductivity. Figure 14 shows the hydrostratigraphic units for the bedrock and unconsolidated material in the 14 model layers.

A starting value for the hydraulic conductivity, and lower and upper bounds on the value were assigned for each zone in the bedrock and unconsolidated materials. During model calibration (described in the next section), the assigned value and bounds were updated to give the model simulation results that better match observed behavior of the system.

Model Target Groups

Using groundwater level observations and streamflow measurements, 14 target sets were created to calibrate the Badger model. The sets and observation group names used in the PESTPP–iES calibration are detailed in appendix 4. A high number of target sets was designed to explore the many aspects of the model parameters (for example, aquifer properties, recharge, and conductance) using the site historical data, mostly water levels, measured at various locations, depths, and times. Of particular interest were target sets that helped constrain the flow field near the plumes such as the MODPATH particles in set 10. The different target sets were weighted based on the uncertainty in the observations and to balance group representation relative to the total starting calibration phi. In PESTPP–iES nomenclature, phi is the sum of the squared weighted residuals between the observed and modeled values for all the target groups and is what the calibration process seeks to minimize. Figure 15 shows the relative contribution of each of the target groups (described in detail in app. 4, table 4.1; search using the set numbers) as a percentage of the initial model phi. Full details of the observation weighting scheme used for calibration can be found in the flow model archive (Reeves and Corson-Dosch, 2023)

Model Parameters

Model parameters include model inputs constrained to a range of values based on prior knowledge but are recognized as having some degree of uncertainty because of sparse measurement, measurement error, and structural errors inherent to the model (such as parameter simplification artifacts). Operationally, using multipliers as the adjustable parameters for the model has been efficient. The PESTPP–iES program is used to adjust the multipliers. In a preprocessing step, initial property values are multiplied by the multipliers to update model property values. The initial values were hand-adjusted before the PESTPP–iES calibration starting with values either inherited from the Columbia County model or from the literature. Most of the model parameters for the Badger model are pilot-point multipliers for hydraulic conductivity for the 14 hydrostratigraphic units (table 3). A pilot point is a specified location within a zone or unit where a value is adjusted; pilot points provide an estimate for the parameterized property that can change within the zone or unit. Other parameters include multipliers for recharge by zone, storage properties, and drain conductance. The only parameters that are not multipliers on initial values, but rather relate to values themselves, are parameters for the vertical anisotropy factor for hydraulic conductivity by zone, river conductance for five sections along the Wisconsin River, and porosity by model layer (table 3). The porosity parameters were not used in the flow model calibration.

To simulate groundwater flow, hydraulic conductivity must be assigned for every finite-difference cell in the model. The hydrostratigraphic units (fig. 14) for the unconsolidated material and bedrock at the site discussed previously and detailed in appendix 3 were used to guide the initial estimates for Kh and Kv for the model. Some units only cover part of the model area and, thus, have few observations to inform how their properties vary spatially. For these units, a single hydraulic conductivity multiplier was used. Other units, like the sand, cover a large part of the model area and are expected to have variation in hydraulic conductivity within the unit. For these units, pilot points were used. Pilot-point parameters are multiplier values applied to initial hydraulic conductivity value at the locations shown in figure 16. The starting pilot-point multiplier values were 1, which maintains the original hydraulic conductivity field. The PEST groundwater utility FAC2REAL was used to interpolate from updated hydraulic conductivity estimates at the pilot points to a horizonal two-dimensional array of values for each finite-difference cell using kriging (Watermark Numerical Computing, 2020). PARM3D (Watermark Numerical Computing, 2020) was used to interpolate between model layers, and, in this case, assigned the same value for cells in the same hydraulic conductivity zone for all model layers at a given x-y location. For the Badger model, a grid of 1,254 evenly spaced pilot points was established (fig. 16), and the distance between the pilot points was 1,500 ft.

The Kv was estimated by applying an anisotropy ratio (Kv:Kh) for the 14 hydrostratigraphic units. The anisotropy was uniform within each of the 14 units, so this parameter group had 14 parameters. The initial anisotropy assigned to most layers was 10 except for the Eau Claire unit, which was set to 100; and the weathered Wonewoc and Eau Claire units, which were set to 1.

The estimated recharge from the SWB model was considered uncertain and allowed to change by recharge zone. Appendix 2 describes the 16 recharge zones and the process used to determine an initial estimate for each zone. The updated recharge values are discussed in the Calibration Results section in this report.

Aquifer storage properties of specific storage and specific yield are needed for the transient model. However, data to constrain these properties are scarce, and the annual time spans for model stress periods also make it difficult to refine these properties in the calibration process. Therefore, only 6 aquifer storage parameters were used, and they are all multipliers against the initial values from Columbia County model: a single Ss and Sy for the unconsolidated units (z5–z11; table 3), a single Ss and Sy for the bedrock units (z1–z4; table 3), and a single Ss and Sy for the transitional and weathered bedrock units (z12–z14; table 3).

Wetlands and small ponds in the model area were simulated using the MODFLOW DRN package. A single multiplier for the drain conductance for each feature was assigned as a parameter. The Wisconsin River strongly affects the groundwater flow patterns in the model area. It was simulated in the model using the MODFLOW RIV package that requires river stage and riverbed conductance as inputs. The riverbed conductance was considered uncertain, and five zones (fig. 13) along the river were assigned different conductance parameters. The area immediately upstream from the Prairie du Sac hydroelectric dam was assumed to have the lowest riverbed conductance because the strong hydraulic gradient from the river towards the aquifer and lower water velocities could lead to more deposition of fines along the riverbed and, therefore, a lower hydraulic conductivity. Two other sections were imposed upstream from the dam: one section in the likely transition area where the hydraulic gradients switch from gaining to losing, and one section far upstream where the Wisconsin River is gaining and a higher initial value of riverbed conductance is assumed. One section was placed in the immediate tailwater of the dam where more turbulent flow likely flushed fine streambed material and the other was placed downstream from the tailwater section to represent normal streambed conditions. The parameters for these sections are direct riverbed conductance values passed to the groundwater flow model rather than multipliers.