Hydrologic Budget

Much can be learned about a hydrologic system through analysis of the soil-water balance. Thornthwaite (1948) recognized the value of observing the “march of actual evapotranspiration” over the course of a year to quantify the growing-season periods where natural precipitation inputs are insufficient to balance evapotranspiration needs of crops and other vegetation. The components of the soil-water-balance are used here to estimate net infiltration—the amount of soil moisture that passes through the root zone and into the unsaturated soil zone that lies beneath it in most places throughout the Central Sands region. With time, most of this net infiltration eventually flows to the water table, becoming groundwater recharge. Net infiltration estimates are input directly into the groundwater-flow models described later in the report.

Net infiltration is one of the most difficult components of the hydrologic budget to estimate; limited direct observations of net infiltration are available. The use of a soil-water-balance approach depends on properly defining the remaining components of the water budget; net infiltration may then be calculated as the difference between the inputs and outputs of water to the root zone soil layer. A simplified description of the root zone soil-water balance is $net infiltration=rainfall+snowmelt+irrigation-evapotranspiration-runoff- ∆soil moisture$, where $∆soil moisture$ is the change in soil moisture from one day to the next (Westenbroek and others, 2010, 2018).

Annual precipitation in the Central Sands region ranges from less than 30 inches [in.] in 2012 to more than 114.3 cm (45 in. in 2018 [fig. 2]). For reference, the 30-year mean annual precipitation (calculated for the years 1981 through 2010) ranges between 81.3 and 86.4 cm (32 and 34 in., respectively) across the Central Sands region (PRISM Climate Group, 2020). The mean annual precipitation amount received over the Central Sands region during 2016–18 was well above the 30-year mean annual precipitation amount (fig. 2).

Actual evapotranspiration (ET), the sum of bare soil evaporation and plant transpiration, varies widely with land use, crop or vegetative cover type, and the time of year. Actual ET over the region is generally low—less than 1 in.—for January, February, and March (Reitz and others, 2017a). Actual ET increases steadily with springtime vegetation and crop growth, peaking in about July at about 5 in. per month. Actual ET continues to decrease through August, September, and October, and is generally less than 1 in. per month for November and December.

Runoff in the Central Sands region is generally low but does affect the water budget. Annual mean runoff amounts for the study area range between 2 and 4 in., respectively; Reitz and others, 2017b); the years with larger runoff amounts than other years are the same years for which gross precipitation is higher than the 30-year “normal” amounts.

Irrigation amounts, as estimated by means of the Food and Agriculture Organization of the United Nations irrigation and drainage paper 56 (FAO-56) methodology (Allen and others, 1998), vary appreciably from year to year and from one crop type to another. Use of reported or metered irrigation amounts in the water-budget calculations would be preferable to use over estimated values, but reported values were only available for the history-matching period of 2012–18. For the longer range of climatic conditions required for the Lake Ecosystem Response Assessment (LERA) scenarios (1981–2018) simulated in the lake study, irrigation amounts computed by Soil Water Balance (SWB) (see SWB subsection of the “Numerical Models” section for a detailed description) were used. Irrigation amounts represent between 0.8 and 1.3 cm (0.3 and 0.5 inch per year, respectively) of water when averaged over the entire study area during 2012–18.

Finally, the net annual infiltration amount can be estimated once all the other water-budget components have been accounted for. Net infiltration averaged over the study area ranges from 6 in. in 2012 to 13 in. in 2018. April and October tend to have the highest estimated net infiltration amounts, with June, July, and August having the lowest amounts.

Aquifer Properties

The surficial aquifer ranges in thickness from thin to absent near topographic high points in the bedrock surface, to more than 100 m thick along the Almond Moraine (fig. 1) east of Plainfield, Wisconsin. Horizontal hydraulic conductivity in the surficial aquifer ranges from less than 1 m per day locally to 150 m per day or more, and values of 10–80 m per day are common (Weeks and Stangland, 1971; Hart and others, 2015; Bradbury and others, 2017). Horizontal hydraulic conductivity is generally higher in the well-sorted, undisturbed outwash west of the moraines, especially in the northwestern part of the study area. Vertical to horizontal anisotropy is commonly between 10 and 100. Specific yield in the surficial aquifer ranges from about 0.12 to 0.33, with a mean of about 0.17 (Weeks and Stangland, 1971; Hart and others, 2015; Bradbury and others, 2017).

Plainfield Tunnel Channel Lakes

The Plainfield Tunnel Channel Lakes area is focused on two of the study lakes, Plainfield and Long Lakes. Both lakes are shallow (no more than 15 ft deep) and were formed in a collapsed tunnel channel (for example, WDNR, 2021, app. A). The tunnel channel area near Plainfield contains several lakes in addition to the two (Plainfield and Long Lakes) considered in this area for this study. Hydraulically, the Plainfield Tunnel Channel Lakes are near the groundwater divide, in an area with no streams. Regional groundwater flow in this area is from north to south and towards the stream networks to the east and west (Central Sands model: Mechenich, 2012). Throughflow is the typical flow condition in the Plainfield Tunnel Channel Lake. Groundwater discharges into the north side of both lakes, and leakage flows out of the southern side of both lakes into the surficial aquifer (as described later in the “Regional Model Results” section). However, the lack of streams in this area results in greater fluctuations in groundwater levels in response to changes in recharge and pumping than in areas with more streams, where surface water buffers water-table fluctuations. This result is consistent with observations in seepage lakes in a similar glacial environment in Cape Cod, Massachusetts (McCobb and others, 1999), where the faster response of the surrounding groundwater system relative to the lake results in reversal of flow into or out of the lake depending on whether surrounding groundwater levels are rising or falling. In the spring of 2019, after an extended rise in the water table (WDNR, 2021, app. A), a transition to nearly all groundwater discharge around the perimeter of Plainfield and Long Lakes was observed, making the lakes act as sinks to the groundwater system. At other times, when water levels are declining, the lakes may leak around their perimeters and act as sources of water to the groundwater system. Variable groundwater levels and shallow bathymetries make the Plainfield Tunnel Channel Lakes prone to large changes in exposed shoreline, volume and surface area, and even periodic drying.

Pleasant Lake

The third study lake—Pleasant Lake—is relatively deep (nearly 25 feet) compared to the Plainfield Tunnel Channel Lakes and situated in a collapse feature just east of the moraine, near two headwater streams (Chaffee Creek and Tagatz Creek) (fig. 1). Regional groundwater flow in this area is to the east-southeast, away from the undrained area toward Chaffee Creek to the east and Tagatz Creek to the south as described later in the “Regional Model Results” section. The typical flow condition of Pleasant Lake is, therefore, one of groundwater discharge into the western and northern sides and leakage out to the groundwater system on the southern and eastern sides. Hydraulic gradients around Pleasant Lake are generally stable because of nearby streams and because the lake is downgradient from the groundwater divide. The steeper sides and greater depth of the Pleasant Lake Basin also make the lake shoreline area less sensitive to seasonal water-level fluctuations. Unlike the Plainfield Tunnel Channel Lakes, where depth to bedrock is 164 ft or more, bedrock depth around Pleasant Lake is much shallower, approximately 30 ft, near the northern and southeastern edges of the lake, and cropping out in nearby areas south of the lake as discussed further by WDNR (WDNR, 2021, app. A). Shallow bedrock and steeper topography around Pleasant Lake (than around the other lakes considered) result in more surface runoff than in the Plainfield Tunnel Channel Lakes. Shallow bedrock also results in decreased groundwater exchange with Pleasant Lake.

Model Domain

The regional model boundary is shown in figure 1. The active model extent was selected to include natural flow boundaries, where available, and all pumping stresses that may affect the hydraulic heads at the specified boundaries of the two lake inset models. Regional model flow boundaries are the Wisconsin and Plover Rivers to the west; the Tomorrow, Waupaca, and Wolf Rivers to the north; and the Fox River to the east. Areas beyond the regional boundaries were inactive in the model.

Model Grid and Layering

The regional model consists of a uniform grid of 572 rows and 533 columns of 200-m cells. The model contains four layers to represent the hydrostratigraphic units discussed in the “Conceptualization of Groundwater Flow in the Central Sands Region” section and includes the following:

The top elevation of layer 1 is resampled from a 10-m digital elevation model (DEM) of the model domain (WDNR, 2019). The bottom elevations of layers 1 and 2 are defined by the bottom of the upper coarse layer and middle fine layer discussed in WDNR (2021, app. A). The bottom of layer 3 is the top of the sandstone bedrock unit, and the bottom of the layer 4 is the top of the Precambrian bedrock. Additional geologic information concerning the model layers is available from WDNR (WDNR, 2021, app. A).

Boundary Conditions

Boundary conditions that affected groundwater movement in the regional model included infiltration originating from precipitation and snowmelt, groundwater exchanges through the model perimeter and the stream network, and groundwater withdrawals through pumping. Boundary conditions used in the regional model are shown in figure 3.

Spatially averaged net infiltration (aquifer recharge) for the 2012–18 period ranged from about 3.9–23.6 in/y in the regional model domain (fig. 4). This net infiltration was specified using the MODFLOW Recharge (RCH) package (Harbaugh, 2005).

High-capacity pumping wells operating during any part of the 2012–18 period were included in the regional model (fig. 3). Pumping was specified using the MODFLOW Well (WEL) package (Harbaugh, 2005). Well construction and reported monthly pumping information for high-capacity wells were provided by the WDNR (https://dnr.wisconsin.gov/topic/WaterUse/data.html). Well locations were assigned to the model layer with highest transmissivity within the well’s open interval. Wells without open-interval information were assigned to the model layer with the highest transmissivity at their location.

Lateral flow boundaries were represented in the regional model using the MODFLOW General Head Boundary (GHB) package (Harbaugh, 2005). These flow boundaries included the Wisconsin and Plover Rivers to the west, the Tomorrow, Waupaca, and Wolf Rivers to the north, and the Fox River to the east (figs. 1 and 3). Groundwater exchanges with GHB cells are computed using the difference between the simulated hydraulic head in the aquifer and the assigned GHB elevation, and a conductance term. The elevation of the GHB cells was set to the minimum DEM elevation (WDNR, 2019) within the model-cell area. The conductance term is a function of the cell area, the assumed thickness of the riverbed, and the vertical hydraulic conductivity of the riverbed material. GHB cells were assigned a conductance of 0.5 meter squared per day (m²/d), which, assuming a 1-m thick riverbed, is equivalent to a 1.25×10⁻⁵ meter per day (m/d) vertical hydraulic conductivity. This conductance was adjusted during history matching to 3.9 m²/day, which is equivalent to a 9.75×10⁻⁵ m/d (3.20×10⁻⁴ feet per day [ft/d]) vertical hydraulic conductivity. This value is similar to the calibrated vertical hydraulic conductivity of 4×10⁻⁴ ft/d for GHB cells representing Lake Superior in Leaf and others (2015).

Lateral flow boundaries form the perimeter of the active area of the regional model except along the southern boundary. The southern model boundary has no major river nearby and is far enough from the study lakes as not to affect groundwater flow near the lakes. As a result, the southern model boundary was set as a no-flow boundary, consistent with a previous two-layer groundwater-flow model (Mechenich, 2012) of the Central Sands region (referred to here as the Mechenich Model), located parallel to the assumed groundwater-flow direction east to the Fox River or west to the Wisconsin River.

Streams were represented in the regional model using the MODFLOW-NWT Streamflow Routing (SFR2) package (fig. 3) (Niswonger and Prudic, 2005). SFR2 calculates exchanges between the aquifer and stream system while accounting for total streamflow in the channel. SFR2 input was developed from National Hydrography Dataset version 2 hydrography using the SFRmaker software (Leaf and others, 2021). National Hydrography Dataset flowline features were joined spatially to the model grids, with the resulting line/grid cell intersections and associated attributes forming the basis for SFR2 reaches. Streambed-top elevations were sampled from the light detection and ranging (lidar)-based DEM (WDNR, 2019) as the minimum elevation within a 100-m buffer around each SFR2 reach. Similar to the GHB, SFR2 has a conductance term that represents the thickness of the streambed, the model-cell area of the SFR2 cell, and the vertical hydraulic conductivity of the streambed. The vertical hydraulic conductivity of the streambed was estimated for each stream reach during history matching. Vertical hydraulic conductivity values of the streambed started (initial model-simulation run) at 1 m/d with final calibrated (history matching) values ranging from 0.03 to 80.5 m/d.

Aquifer Properties

Aquifer properties represented in the regional model include vertical and horizontal hydraulic conductivity, specific storage, and specific yield. Initial values were assigned based on literature values described throughout the report and then adjusted during history matching.

The initial horizontal hydraulic conductivity for the eastern half of the unconsolidated model layers (1–3) was estimated using the coarse/fine fractions from WDNR (WDNR, 2021, app. A) and a power-law approach similar to Feinstein and others (2010). The power-law parameters were optimized so that the mean of the initial hydraulic conductivity values produced by the power law agreed with the mean values from the Mechenich Model. Hydraulic conductivity for the western half of the unconsolidated layers, west of the terminal moraine, where coarse/fine fraction estimates were not made in WDNR (2021, app. A), used information from the upper model layer, representing unconsolidated materials, of the two-layer Mechenich Model of the Central Sands region. The initial hydraulic conductivity of the sandstone bedrock in the regional model (layer 4) came from layer 2, representing bedrock, in the Mechenich Model. Lakes across the model domain were included as high hydraulic conductivity zones in layer 1 with a horizontal hydraulic conductivity of 10,000 m/d.

Initial vertical hydraulic conductivity of the fine-grained New Rome Member simulated in layer 2 was set to 0.01 m/d based on vertical hydraulic conductivity estimates for the New Rome Member of 0.2 and 8×10⁻⁵ m/d from Hart and others (2015). The initial horizontal hydraulic conductivity of the New Rome Member was set at 0.1 m/d, which is 10 times the vertical hydraulic conductivity. The specific yield for the New Rome Member was set to 0.20 based on a range of literature values for silt, clay, and fine sand specific yields of 0.06–0.33 (Duffield, 2019). The specific storage was set to 9.2×10⁻⁴m⁻¹, slightly less than the Hart and others (2015) estimate of 9.8×10⁻⁴ m⁻¹ for the transition zone of the New Rome and at the lower end of the literature range of 9.2×10⁻⁴–2.0×10⁻² m⁻¹ (2.8×10⁻⁴–6.2×10⁻³ft⁻¹) for clays (Duffield, 2019). The aquifer properties of the New Rome Member after history matching were a horizontal hydraulic conductivity of 9.2×10⁻² m/d, vertical hydraulic conductivity of 2.0×10⁻³ m/d, a specific yield of 0.2, and a specific storage of 9.2×10⁻⁴ m⁻¹. The horizontal hydraulic aquifer properties of the regional model after history matching are shown in figure 5.

Horizontal hydraulic conductivities ranged from 0.09 to 152 m/d in the three unconsolidated model layers (layers 1–3) (fig. 5). The mean horizontal hydraulic conductivity was highest in layers 1 and 3 (46 and 18 m/d, respectively) and lowest in layer 2 (13 m/d). These values are consistent with layers 1 and 3 representing coarser glacial material and layer 2 representing finer material to the east and the New Rome Member to the west. The mean horizontal hydraulic conductivity was 9 m/d (30 ft/d) for the sandstone bedrock, represented by model layer 4.

Vertical hydraulic conductivities of the three unconsolidated model layers ranged from 0.002 to 2.93 m/d, with a mean that was lowest in layer 2 (0.37 m/d) and highest in layer 1 (1.01 m/d) (fig. 6). The vertical hydraulic conductivity of the sandstone bedrock (layer 4) had a mean of 0.13 m/d (0.43 ft/d).

The specific yield means for the three unconsolidated layers (1–3) ranged from 0.19 to 0.37 and was 0.14 for the bedrock sandstone (layer 4) (fig. 7). These values are within the literature ranges for till, clay, sand, and gravel specific yields of 0.06–0.33 and within the literature ranges for the sandstone bedrock of 0.06–0.27 (Duffield, 2019). The mean specific storage for the three unconsolidated layers ranged from about 3.3×10⁻⁴ to 4.8×10^-4m⁻¹ (fig. 8). Model values for the unconsolidated material (layers 1–3) fell within the specific storage literature ranges (Duffield, 2019) for loose to dense sand and gravel (from 1.0×10⁻⁴ to 4.3×10⁻⁴ m⁻¹) and clays (from 9.2×10⁻⁴ to 2.0×10⁻²m⁻¹) and close to the Hart and others (2015) estimate of specific storage in the sandy aquifer underlying the New Rome Member of 3.0×10⁻⁴m⁻¹ (9×10⁻⁵ft⁻¹). The specific storage was about 9.8×10⁻⁵m⁻¹(3.0×10⁻⁵m⁻¹) for the sandstone bedrock (layer 4) and is just above the literature ranges (Duffield, 2019) for fissured rock from 3.3×10⁻⁶ to 6.9×10⁻⁵m⁻¹.

History-Matching Approach and Results

History matching was performed with the regional model. History matching is a process referred to as “parameter estimation” and, sometimes, “model calibration,” although “calibration” implies a precise unique fit of the model to data. As described herein, calibration is not an accurate representation of the history-matrching process to systematically adjust parameter values in the model such that associated model outputs are consistent with historical observations including hydraulic-head and base-flow values. The general history matching approach follows Bayes’ theorem (Tarantola, 2005) for parameter estimation. In this approach, a prior estimate of model parameters and their credible ranges (uncertainty) is first used. These prior parameter values and uncertainty are selected by expert knowledge, literature values, and available direct measurements. Through a systematic conditioning step, the parameter values are updated to be consistent with observations that correspond to model outputs. This step is referred to as an “update” and results in a posterior set of parameter values (“posterior” meaning “after the update”). Many algorithms can be used to perform this conditioning step. For the regional model, the parameter estimation package for high-performance computing ([PEST++]; White and others, 2021) software was used. The choice of software, in part, reflects rapid development of PEST-related tools throughout the study timeframe. The iterative ensemble smoother (iES) history matching algorithm in White and others (2021) was used for the inset models. The history matching files are provided in Fienen and others (2021b).

The measured values of groundwater hydraulic heads and streamflow used for history matching (referred as “targets”) came from various data sources that are summarized in table 1. Data from a total of 177 streamflow and 464 well and lake-level observation locations were used during the history matching. Streamflow measurements were processed to focus on base flow using the techniques of Sloto and Crouse (1996) and Wahl and Wahl (1988). Simulated groundwater elevations, lake-level elevations, and streamflow from the steady-state period at the start of the transient model were compared to data targets averaged over the 2012–18 period. Simulated groundwater elevations, lake-level elevations, and streamflow from the transient stress periods were compared with targets closest to the end of the month represented by each stress period, for stress periods where measured data were available.

A post-history-matching comparison between the measured and simulated values for lake and groundwater-elevation targets are shown on figure 9, and the streamflow targets are shown on figure 10. The closer these values plot along the 1:1 line, the closer the match between the simulated and measured values. Groundwater-elevation targets mostly plotted on or near the 1:1 line with some outliers scattered above and below the line. The streamflow targets generally plotted on or near the 1:1 line with a slight under-estimated bias (simulated lower than measured values), particularly in comparison to the WDNR base-flow data (fig. 10).

The residuals (differences between simulated and measured values) for the steady-state hydraulic-head and streamflow targets are shown spatially in figures 11 and 12, respectively. Most hydraulic-head targets are from well-construction reports, with accurate locations limited to the central eastern part of the model domain, where location-corrected values were provided by the Wisconsin Geological and Natural History Survey (WDNR, 2021, app. A), focused close to the study lakes. Hydraulic-head residuals were generally well distributed throughout the model domain with some locally biased over-estimated values (simulated higher than measured values) in the central part of the domain. A comparison of flux (streamflow) target residuals indicates a general pattern of under-estimated streamflows, although the bias is less near the inset model areas (for example, Upper Pine River, Chaffee Creek, and Tagatz Creek; fig. 12).

Transient hydraulic-head and streamflow targets for a select subset of wells and streams are shown in figures 13 and 14, respectively. Well locations with with long-term records (covering the entire calibration period) (fig. 13) were selected to represent wells across the model domain and included a range of good to poorly matched simulated and measured water levels. Overall, the simulated transient hydraulic heads at most wells were in good agreement with measured long-term groundwater-level data.

Few stream locations had daily streamflow data in the model domain. Tenmile Creek (fig. 12) has the most complete daily record for the history-matching period, and the simulated fluxes generally were in good agreement with the trends and timing of peak flows for the period of record. The simulated fluxes generally were lower than the measured streamflow highs but were in good agreement with the mid and lower base flows. Simulated fluxes in the headwater streams near the focus area lakes reasonably matched the magnitude and trends measured in these streams, except for Tagatz Creek, where the simulated base flows in the headwaters were appreciably lower than measured base flows. This resulted, in part, because of bedrock anomalies in the area around Tagatz Creek. These anomalies were addressed in the inset models by locally adjusting the bedrock elevation near the creek to simulate incision of the stream.

Model parameters (totalling 1,775) were adjusted during the history matching and included:

The aquifer properties included in the history matching were horizontal and vertical hydraulic conductivity, specific storage, and specific yield. These properties had initial values assigned as discussed in the “Aquifer Properties” section. These initial parameter arrays were then adjusted using a network of pilot points representing multipliers and zone multipliers for each layer. Zones and pilot points are shown in figure 15 for each of the four regional model layers and include:

The aquifer properties after history matching are presented and discussed in the “Aquifer Properties” section. The GHB conductance is discussed in the “Boundary Conditions” section. The final multipliers on the SWB-estimated net infiltration ranged from 0.8 to 1.2 with a mean of 1.02, meaning that, on average, the adjustment to net infiltration was a 2 percent increase over what was estimated with SWB. The final multipliers on the reported well-pumping rates ranged from 0.75 to 1.25 with a mean of 0.95, meaning that, on average, the adjustment to reported pumping rates was a 5 percent decrease in the reported rates.

Regional Model Results

The water-table elevation contours and cumulative simulated streamflow for the regional model are presented in figure 16. The water table is mounded along the topographic high formed by the moraine that runs north-south through the center of the model domain. Groundwater flows from this water table high to the east and west towards the major rivers that form the major hydrologic boundaries in the region. Groundwater-flow directions are assumed to be generally perpendicular to water-table elevation contours.

A simulated water budget of the major changes in inflows and outflows from 2012 to 2018 is shown in figure 17. Water enters the groundwater system through recharge and stream leakage under losing streamflow conditions and exits through streams under gaining streamflow conditions, pumping from high-capacity wells, and discharge to the major rivers at the lateral boundaries of the groundwater system. Excess inflows and outflows are balanced throughout the system by groundwater storage replenishment and depletion. During periods of high recharge (often in the spring), groundwater storage has a gain, which can be thought of as excess water leaving the groundwater system and entering into storage (Storage_in on fig. 17). During the growing season, when pumping is higher and recharge is lower, water is being removed from groundwater storage and entering the groundwater system (Storage_out on fig. 17).

Inset Model Domains and Horizontal Discretization

The inset model domains were designed to encompass the nearest headwater streams to the east and west of the study lakes and most of the high-capacity wells that may affect the lake water balances, while maintaining reasonable model runtimes. A uniform grid with 20-m horizontal resolution in an area surrounding the lakes was selected to adequately represent the detailed bathymetry and shoreline geometry of the lake basins. Initial testing of runtimes with the MODFLOW-NWT code and a 20-m grid spacing over the entire inset model domains resulted in unacceptably long runtimes (on the order of hours).

Because of the long runtimes, an LGR approach was adopted using the multiple model capabilities of MODFLOW 6 (Langevin and others, 2017). The models focused on the lakes each consist of two submodels—the “inset” submodel aligned with the regional model grid at the same 200-m (656.2-ft) resolution and a locally refined “LGR” submodel with a uniform 20-m resolution encompassing a rectangular area around the lake(s) of interest (fig. 18 for Plainfield Tunnel Channel Lakes and fig. 19 for Pleasant Lake). The Plainfield Tunnel Channel Lakes inset submodel contains 90 rows and 110 columns, and the LGR submodel contains 120 rows and 250 columns. The Pleasant Lake inset submodel contains 100 rows and 100 columns, and the LGR submodel contains 100 rows and 120 columns. Within MODFLOW 6, the two submodels for each lake(s) are coupled within the same solution matrix and solved simultaneously as a “simulation” (Langevin and others, 2017). This coupling greatly reduced the overall number of model cells, cutting the runtimes for the base-case simulation (2012–18 period) to approximately 10 minutes.

Inset Model Layering

Layering in the inset models was developed from the same data sources as the regional model; however, the exception being that layer 1 was typically split into two layers for the inset models and was adjusted for lake bathymetry where the lakes are simulated. The model top (top of layer 1) is based on mean elevations sampled for each model cell from a lidar-based DEM (WDNR, 2019), except within the basins of Plainfield, Long, and Pleasant Lakes. Bathymetry for these lakes, developed by WDNR (2019), was subtracted from the DEM elevations to develop the model top. The layer bottom surfaces in the inset models, where lakes are not present, are congruent with the regional model except that layers 1 and 2 are evenly subdivided from layer 1 in the regional model to better represent hydraulic gradients near surface-water features.

Unlike MODFLOW-NWT, MODFLOW 6 allows for discontinuous layering, meaning that model cells can be removed from the model solution in places where a hydrogeologic unit is absent. In the Plainfield Tunnel Channel Lakes inset model, this feature was used to remove cells from layer 3 in areas west of the Hancock Moraine, where “layer 2” from WDNR (2021, app. A) was not defined (fig. 20). Similarly, in the Pleasant Lake inset model, cells were removed from “layer 3” in places where the New Rome Member or “layer 2” from WDNR (2021, app. A) were absent. Cells were also removed from the models in places where the layer bottoms were within a meter of the land surface, for example, near high points in bedrock surface or along lake bottoms. Removing unneeded model cells, which would otherwise need to be carried throughout the model at some minimum thickness (for example, 1 m), can further increase the stability and speed of the model-simulation run.

Time Discretization

Time discretization for the base versions of the inset models is similar to the regional parent model with the exception that the initial steady-state period represents mean conditions for the period from 2012 through 2015. Subsequent transient monthly stress periods represent mean conditions for each month through 2018. Each transient stress period was subdivided into five time steps using a time-step multiplier of 1.2. This time-step multiplier differs from the regional model multiplier of 1.5 because inset model convergence is more difficult with the lake package, and a lower multiplier can sometimes improve model convergence.

Boundary Conditions

Boundary conditions within the inset model domains include regional groundwater flow across the model perimeters; terrestrial recharge originating from precipitation, snowmelt and irrigation; and groundwater/surface-water interactions with lakes and streams. The model-perimeter boundaries were simulated as specified-head values obtained from the regional model.

Aquifer Properties

Horizontal and vertical hydraulic conductivity were initially set at the values estimated for the regional model by history matching. To avoid having parameters adjusted to extreme values in the local-scale models’ history matching process, when applying multipliers to localized values that are already approaching extreme values, specific storage (Ss) was initially set to 1×10⁻⁶ m^−1.. Specific yield (Sy) was initially set uniformly to 0.15 (dimensionless) throughout the model domain, based on previous investigations for similar deposits.

Parameter Estimation and Uncertainty Analysis

The history matching approach used in this analysis is described in detail in Corson-Dosch and others (2022). Similar to history matching with the regional model, history matching was performed to refine parameter estimates for the inset and LGR models with the same observation data used for the regional model analysis, with additional observations in the area around the lakes. For the inset models, the iES implementation (White, 2018) of PEST++ (version 5.0.0; White and others, 2021) was applied. The goal of iES is to provide parameter estimates and to quantify the uncertainty of those estimates. This uncertainty is characterized by a range of estimated parameters that each reproduce simulated hydraulic heads and streamflows within an acceptable range of measured values. Each observation is assigned an observation weight that, in principle, corresponds to the inverse of the standard deviation assuming a normal distribution of uncertainty around the measured observation value. As a result, in comparing the model output with measured values, the uncertainty of the measured and simulated values warrant consideration and, ideally, these values would overlap.

Observation Data

Observation data were compiled for the inset models from the same dataset used for history matching in the regional parent model. Some additional observations focused around the study lakes were added for the inset models and some spatial water-level difference observations (measures of the horizontal hydraulic gradient between shallow piezometers (usually less than 10 m deep) installed near the lakes and the simulated lake elevations) were added along with the temporal difference observations. Observation weights were initially assigned based on assumptions regarding a level of fit between model outputs and nearby observations. These weights were adjusted to balance the objective function that drives the iES regression including assigning a weight of 0.0 to some classes of observations. These adjustments are all intended to steer the iES results toward a model design that is focused on the complex groundwater/surface-water interactions in the study area, and, in particular, matching observed lake elevations and streamflow values—a key objective of this analysis. The observations used, by observation group, and values and observation weights for the Pleasant Lake and Plainfield Tunnel Channel Lakes inset models are summarized in tables 2 and 3, respectively.

Table 2.    

Observation data groups, values, and descriptions for Pleasant Lake inset model, Central Sands region, central Wisconsin.

[PEST, parameter estimation package; WGNHS, Wisconsin Geological and Natural History Survey; NA, not applicable; USGS, U.S. Geological Survey; NWIS, National Water Information System database; WDNR, Wisconsin Department of Natural Resources]

Table 2.    Observation data groups, values, and descriptions for Pleasant Lake inset model, Central Sands region, central Wisconsin.

Group name in the PEST files	Observed values	Target type	Description	Total number of observations	Number of weighted observations	Number of zero-weighted observations	Weight	Weight-informed standard deviation	Data source
hds_wgnhs_tr, heads_wgnhs	280.311 to 327.17	Hydraulic head	Well-construction report that includes groundwater elevation measured after a well was drilled. Locations were determined by the WGNHS.	100	0	100	0	NA	WDNR, 2022
nwis_dvs	298.00 to 314.72	Hydraulic head	Groundwater elevations at locations with daily data collected by USGS.	14	14	0	3.78 to 7.13	0.14 to 0.26	USGS, 2021
nwis_fm	264.09 to 300.84	Hydraulic head	Groundwater elevations at miscellaneous locations measured by USGS.	2	0	2	0	NA	USGS, 2021
nwisdvs_tr	297.75 to 315.47	Hydraulic head	Groundwater elevations at locations with daily data collected by USGS.	133	133	0	1.26 to 2.04	0.49 to 0.79	USGS, 2021
nwisdvs_tr_tdiff	−0.25 to 0.54	Head temporal difference	Groundwater elevations at locations with daily data collected by USGS.	119	119	0	7.32	0.14	Derived from data elsewhere in this table
nwisfm_tr	263.12 to 300.84	Hydraulic head	Groundwater elevations at miscellaneous locations measured by USGS.	5	0	5	0	NA	USGS, 2021
nwisfm_tr_tdiff	0.12 to 0.59	Hydraulic head temporal difference	Groundwater elevations at miscellaneous locations measured by USGS.	3	0	3	0	NA	Derived from data elsewhere in this table
usgs_stages	299.09	Lake level	Lake elevations measured by USGS.	1	1	0	63.67	0.02	USGS, 2021
usgs_stages_tr	298.94 to 299.25	Lake level	Transient lake elevations measured by USGS.	8	8	0	50.46	0.02	USGS, 2021
nwisdvflx_tr	19,963 to 34,446	Streamflow	Streamflow at locations with daily data measured by USGS.	16	16	0	0.0012 to 0.00215134	464.83 to 802.05	USGS, 2021
wdnrflx_tr	0 to 14, 8776	Streamflow	WDNR streamflow measurements made during base-flow conditions. No adjustments made.	171	127	44	0 to 0.0040	247.22 to 6552.83	WDNR, 2022
wdnr_lakes	298.98	Lake level	Lake elevations measured by WNDR.	1	1	0	93.99	0.01	WDNR, 2022
wdnr_lb	298.98 to 299.12	Lakebed head	Lakebed piezometer groundwater elevation measured by WDNR.	10	10	0	6.04	0.17	WDNR, 2022
wdnr_lb_sdiff	–0.11 to 0.05	Lakebed head spatial difference	Lakebed piezometer groundwater elevation measured by WDNR.	10	10	0	8.19	0.12	WDNR, 2022
wdnrlks_ann_tr	298.42 to 299.16	Lake level	Lake elevations measured by WNDR—long term.	6	6	0	61.019	0.016	WDNR, 2022
wdnrlks_mon_tr	298.92 to 299.18	Lake level	Lake elevations measured by WNDR—recent (2017–18)	9	9	0	54.75	0.018	WDNR, 2022
wdnrwells_tr	297.29 to 299.88	Hydraulic head	Wells installed for this study and measured by WGNHS and WDNR.	68	68	0	1.25 to 2.45	0.4084 to 0.80	WDNR, 2022
wdnrwells_tr_sdiff	–0.61 to 1.80	Hydraulic head spatial difference	Wells installed for this study and measured by WGNHS and WDNR.	56	56	0	2.31 to 4.52	0.22 to 0.43	Derived from data elsewhere in this table
wdnrwells_tr_tdiff	–0.062 to 0.29	Hydraulic head temporal difference	Wells installed for this study and measured by WGNHS and WDNR.	56	56	0	12	0.083	Derived from data elsewhere in this table

Table 3.    

Observation data groups, values, and descriptions for Plainfield Tunnel Channel Lakes inset model, Central Sands region, central Wisconsin.

[PEST, parameter estimation package; WGNHS, Wisconsin Geological and Natural History Survey; NA, not applicable; USGS, U.S. Geological Survey; NWIS, National Water Information System database; WDNR, Wisconsin Department of Natural Resources]

Table 3.    Observation data groups, values, and descriptions for Plainfield Tunnel Channel Lakes inset model, Central Sands region, central Wisconsin.

Group name in the PEST files	Observed values	Target type	Description	Total number of observations	Number of weighted observations	Number of zero-weighted observations	Weight	Weight-informed standard deviation	Data source
hds_wgnhs_tr, heads_wgnhs	299.61 to 335.15	Hydraulic head	Well-construction report that includes groundwater elevation measured after a well was drilled. Locations were determined by the WGNHS.	106	0	106	0	NA	WDNR, 2022
nwis_dvs	330.79 to 334.91	Hydraulic head	Groundwater elevations at locations with daily data collected by USGS.	4	4	0	16.38 to 26.16	0.038 to 0.06	USGS, 2021
nwis_fm	310.75 to 327.11	Hydraulic head	Groundwater elevations at miscellaneous locations measured by USGS.	2	0	2	0	NA	USGS, 2021
nwisdvs_tr	329.90 to 335.35	Hydraulic head	Groundwater elevations at locations with daily data collected by USGS.	82	82	0	2.32 to 3.76	0.27 to 0.43	USGS, 2021
nwisdvs_tr_tdiff	−0.26 to 2.01	Hydraulic head temporal difference	Groundwater elevations at locations with daily data collected by USGS.	78	78	0	0.79 to 9.50	0.11 to 1.26	NA
nwisfm_tr	310.44 to 327.46	Hydraulic head	Groundwater elevations at miscellaneous locations measured by USGS.	64	0	64	0	NA	USGS, 2021
nwisfm_tr_tdiff	−0.57 to 0.80	Hydraulic head temporal difference	Groundwater elevations at miscellaneous locations measured by USGS.	62	0	62	0	NA	Derived from data elsewhere in this table.
usgs_stages	330.69 to 335.22	Lake level	Lake elevations measured by USGS.	3	3	0	17.70	0.056	USGS, 2021
wdnrflx_tr	0 to 14434.8	Streamflow	WDNR streamflow measurements made during base flow conditions. No adjustments made.	88	73	15	0 to 0.00808321	123.71 to 1149.53	WDNR, 2022
wdnrlks_ann_tr	333.65 to 334.97	Lake level	Lake elevations measured by WNDR—long term	10	10	0	25.15	0.040	WDNR, 2022
wdnrlks_mon_tr	334.74 to 335.47	Lake level	Lake elevations measured by WNDR—recent (2017–18).	16	16	0	16.62	0.060	WDNR, 2022
wdnrwells_tr	334.82 to 335.58	Hydraulic head	Wells installed for this study and measured by WGNHS and WDNR.	121	121	0	0.10 to 1.95	0.51 to 1.00	WDNR, 2022
wdnrwells_tr_sdiff	−0.16 to 0.36	Hydraulic head spatial difference	Wells installed for this study and measured by WGNHS and WDNR.	105	105	0	3.99 to 7.80	0.13 to 0.25	Derived from data elsewhere in this table.
wdnrwells_tr_tdiff	−0.065 to 0.21	Hydraulic head temporal difference	Wells installed for this study and measured by WGNHS and WDNR.	97	97	0	10.32	0.097	Derived from data elsewhere in this table.

Water level and streamflow observations were assigned their values at the end of the month in which they were measured. In this way, intermediate results at a timescale shorter than the stress periods are not considered. Intermediate results were not considered to be consistent with the way stresses (including pumping and recharge) are implemented in model simulation. With monthly stress periods, all stresses that occur at any time in the month must be applied uniformly at the start of the stress period. As a result, observations at intermediate times within the month are not truly representative of their submonthly timescale location, so the end of the month is considered most representative. For streamflow, total streamflow was used to assign observation value with the assumptions that there may be some stormflow in the observations that are not simulated in the model as only base flow is simulated. However, base flow is assumed to compose most of the streamflow in this area.

The locations of water-level observations for the Pleasant Lake and Plainfield Tunnel Channel Lakes inset models are shown in figures 21 and 22, respectively. The labels correspond to observation names that were used in the PEST++ input files.

Parameterization

The parameterization strategy for the inset models was to use multipliers for most parameters to preserve, to the extent possible, the prior conceptualization of the regional parent model. In this approach, the initial parameter values used in the inset models resulted from the previous history matching in the regional model. The multipliers all were set to an initial value of 1.0 and were given narrow upper and lower bounds for adjustment. In this way, particularly with pilot points that apply to a small spatial footprint in the property fields, only minor parameter adjustments are allowed. This approach avoids overfitting with iES. For specific storage and specific yield, however, a multiscale approach to multipliers was applied following the techniques of McKenna and others (2019) and White and others (2020). In this approach, multipliers to entire zones (or, in this case, model layers) can be combined with multipliers at a finer scale (in this case, pilot points). Streambed conductance and lakebed conductivity were both parameterized using values in the properties’ actual measurement units rather than multipliers. The parameterization for the Pleasant Lake and Plainfield Tunnel Channel Lakes inset models are summarized in tables 4 and 5, respectively. Pilot-point locations for horizontal and vertical hydraulic conductivity, specific storage, and specific yield are described and shown in the “Parameter Estimates” section for the Plainfield Tunnel Channel inset model and Pleasant Lake inset model.

Table 4.    

Parameter groups, values, and descriptions for Pleasant Lake inset model, Central Sands region, central Wisconsin.

[PEST, parameter estimation package; K, hydraulic conductivity; SFR, Streamflow Routing package]

Table 4.    Parameter groups, values, and descriptions for Pleasant Lake inset model, Central Sands region, central Wisconsin.

Group name in the PEST files	Group description	Transform	Count	Initial value	Lower bound	Upper bound
_k33_pp_inset:0	Inset vertical K pilot point multipliers: layer 1	log	99	1.00	0.01	10.00
_k33_pp_inset:1	Inset vertical K pilot point multipliers: layer 2	log	99	1.00	0.01	10.00
_k33_pp_inset:2	Inset vertical K pilot point multipliers: layer 3	log	99	1.00	0.01	10.00
_k33_pp_inset:3	Inset vertical K pilot point multipliers: layer 4	log	99	1.00	0.01	10.00
_k33_pp_inset:4	Inset vertical K pilot point multipliers: layer 5	log	99	1.00	0.01	10.00
_k33_pp_parent:0	Parent vertical K pilot point multipliers: layer 1	log	81	1.00	0.01	10.00
_k33_pp_parent:1	Parent vertical K pilot point multipliers: layer 2	log	81	1.00	0.01	10.00
_k33_pp_parent:2	Parent vertical K pilot point multipliers: layer 3	log	81	1.00	0.01	10.00
_k33_pp_parent:3	Parent vertical K pilot point multipliers: layer 4	log	81	1.00	0.01	10.00
_k33_pp_parent:4	Parent vertical K pilot point multipliers: layer 5	log	81	1.00	0.01	10.00
_k_pp_inset:0	Inset horizontal K pilot point multipliers: layer 1	log	99	1.00	0.01	10.00
_k_pp_inset:1	Inset horizontal K pilot point multipliers: layer 2	log	99	1.00	0.01	10.00
_k_pp_inset:2	Inset horizontal K pilot point multipliers: layer 3	log	99	1.00	0.01	10.00
_k_pp_inset:3	Inset horizontal K pilot point multipliers: layer 4	log	99	1.00	0.01	10.00
_k_pp_inset:4	Inset horizontal K pilot point multipliers: layer 5	log	99	1.00	0.01	10.00
_k_pp_parent:0	Parent horizontal K pilot point multipliers: layer 1	log	81	1.00	0.01	10.00
_k_pp_parent:1	Parent horizontal K pilot point multipliers: layer 2	log	81	1.00	0.01	10.00
_k_pp_parent:2	Parent horizontal K pilot point multipliers: layer 3	log	81	1.00	0.01	10.00
_k_pp_parent:3	Parent horizontal K pilot point multipliers: layer 4	log	81	1.00	0.01	10.00
_k_pp_parent:4	Parent horizontal K pilot point multipliers: layer 5	log	81	1.00	0.01	10.00
_ss_con_inset__multiplier	Inset specific storage multipliers by layer	log	5	1.00	0.70	100.00
_ss_con_parent__multiplier	Parent specific storage multipliers by layer	log	5	1.00	0.70	100.00
_ss_pp_inset:0	Inset specific storage pilot point multipliers: layer 1	log	99	1.00	0.70	100.00
_ss_pp_inset:1	Inset specific storage pilot point multipliers: layer 2	log	99	1.00	0.70	100.00
_ss_pp_inset:2	Inset specific storage pilot point multipliers: layer 3	log	99	1.00	0.70	100.00
_ss_pp_inset:3	Inset specific storage pilot point multipliers: layer 4	log	99	1.00	0.70	100.00
_ss_pp_inset:4	Inset specific storage pilot point multipliers: layer 5	log	99	1.00	0.70	100.00
_ss_pp_parent:0	Parent specific storage pilot point multipliers: layer 1	log	81	1.00	0.70	100.00
_ss_pp_parent:1	Parent specific storage pilot point multipliers: layer 2	log	81	1.00	0.70	100.00
_ss_pp_parent:2	Parent specific storage pilot point multipliers: layer 3	log	81	1.00	0.70	100.00
_ss_pp_parent:3	Parent specific storage pilot point multipliers: layer 4	log	81	1.00	0.70	100.00
_ss_pp_parent:4	Parent specific storage pilot point multipliers: layer 5	log	81	1.00	0.70	100.00
_sy_con_inset__multiplier	Inset specific yield multipliers by layer	log	5	1.00	0.75	2.00
_sy_con_parent__multiplier	Parent specific yield multipliers by layer	log	5	1.00	0.75	2.00
_sy_pp_inset:0	Inset specific yield pilot point multipliers: layer 1	log	99	1.00	0.75	2.00
_sy_pp_inset:1	Inset specific yield pilot point multipliers: layer 2	log	99	1.00	0.75	2.00
_sy_pp_inset:2	Inset specific yield pilot point multipliers: layer 3	log	99	1.00	0.75	2.00
_sy_pp_inset:3	Inset specific yield pilot point multipliers: layer 4	log	99	1.00	0.75	2.00
_sy_pp_inset:4	Inset specific yield pilot point multipliers: layer 5	log	99	1.00	0.75	2.00
_sy_pp_parent:0	Parent specific yield pilot point multipliers: layer 1	log	81	1.00	0.75	2.00
_sy_pp_parent:1	Parent specific yield pilot point multipliers: layer 2	log	81	1.00	0.75	2.00
_sy_pp_parent:2	Parent specific yield pilot point multipliers: layer 3	log	81	1.00	0.75	2.00
_sy_pp_parent:3	Parent specific yield pilot point multipliers: layer 4	log	81	1.00	0.75	2.00
_sy_pp_parent:4	Parent specific yield pilot point multipliers: layer 5	log	81	1.00	0.75	2.00
lak_mults	Lake Evaporation/precipitation multipliers by stress period	log	170	1.00	0.90	1.10
lakebedk	Lakebed conductivity values by littoral or profundal zone	log	2	0.02 to 0.04	1.00e-04 to 1.00e-02	0.10 to 5.00
rch_child__multiplier	Inset recharge multipliers by stress period	log	85	1.00	0.80	1.20
rch_parent__multiplier	Parent recharge multipliers by stress period	log	85	1.00	0.80	1.20
sfrk	SFR conductance values by reach	log	48	1.00	0.50	2.00
wel__parent_cnst	Well pumping multipliers by stress period	log	85	1.00	0.80	1.20

Table 5.    

Parameter groups, values, and descriptions for Plainfield Tunnel Channel Lakes inset model, Central Sands region, central Wisconsin.

[PEST, parameter estimation package; K, hydraulic conductivity; SFR, Streamflow Routing package]

Table 5.    Parameter groups, values, and descriptions for Plainfield Tunnel Channel Lakes inset model, Central Sands region, central Wisconsin.

Group name in the PEST files	Group description	Transform	Count	Initial value	Lower bound	Upper bound
_k33_pp_inset:0	Inset vertical K pilot point multipliers: layer 1	log	264	1.00	0.01	10.00
_k33_pp_inset:1	Inset vertical K pilot point multipliers: layer 2	log	264	1.00	0.01	10.00
_k33_pp_inset:2	Inset vertical K pilot point multipliers: layer 3	log	264	1.00	0.01	10.00
_k33_pp_inset:3	Inset vertical K pilot point multipliers: layer 4	log	264	1.00	0.01	10.00
_k33_pp_inset:4	Inset vertical K pilot point multipliers: layer 5	log	264	1.00	0.01	10.00
_k33_pp_parent:0	Parent vertical K pilot point multipliers: layer 1	log	80	1.00	0.01	10.00
_k33_pp_parent:1	Parent vertical K pilot point multipliers: layer 2	log	80	1.00	0.01	10.00
_k33_pp_parent:2	Parent vertical K pilot point multipliers: layer 3	log	80	1.00	0.01	10.00
_k33_pp_parent:3	Parent vertical K pilot point multipliers: layer 4	log	80	1.00	0.01	10.00
_k33_pp_parent:4	Parent vertical K pilot point multipliers: layer 5	log	80	1.00	0.01	10.00
_k_pp_inset:0	Inset horizontal K pilot point multipliers: layer 1	log	264	1.00	0.01	10.00
_k_pp_inset:1	Inset horizontal K pilot point multipliers: layer 2	log	264	1.00	0.01	10.00
_k_pp_inset:2	Inset horizontal K pilot point multipliers: layer 3	log	264	1.00	0.01	10.00
_k_pp_inset:3	Inset horizontal K pilot point multipliers: layer 4	log	264	1.00	0.01	10.00
_k_pp_inset:4	Inset horizontal K pilot point multipliers: layer 5	log	264	1.00	0.01	10.00
_k_pp_parent:0	Parent horizontal K pilot point multipliers: layer 1	log	80	1.00	0.01	10.00
_k_pp_parent:1	Parent horizontal K pilot point multipliers: layer 2	log	80	1.00	0.01	10.00
_k_pp_parent:2	Parent horizontal K pilot point multipliers: layer 3	log	80	1.00	0.01	10.00
_k_pp_parent:3	Parent horizontal K pilot point multipliers: layer 4	log	80	1.00	0.01	10.00
_k_pp_parent:4	Parent horizontal K pilot point multipliers: layer 5	log	80	1.00	0.01	10.00
_ss_con_inset__multiplier	Inset specific storage multipliers by layer	log	5	1.00	0.70	100.00
_ss_con_parent__multiplier	Parent specific storage multipliers by layer	log	5	1.00	0.70	100.00
_ss_pp_inset:0	Inset specific storage pilot point multipliers: layer 1	log	264	1.00	0.70	100.00
_ss_pp_inset:1	Inset specific storage pilot point multipliers: layer 2	log	264	1.00	0.70	100.00
_ss_pp_inset:2	Inset specific storage pilot point multipliers: layer 3	log	264	1.00	0.70	100.00
_ss_pp_inset:3	Inset specific storage pilot point multipliers: layer 4	log	264	1.00	0.70	100.00
_ss_pp_inset:4	Inset specific storage pilot point multipliers: layer 5	log	264	1.00	0.70	100.00
_ss_pp_parent:0	Parent specific storage pilot point multipliers: layer 1	log	80	1.00	0.70	100.00
_ss_pp_parent:1	Parent specific storage pilot point multipliers: layer 2	log	80	1.00	0.70	100.00
_ss_pp_parent:2	Parent specific storage pilot point multipliers: layer 3	log	80	1.00	0.70	100.00
_ss_pp_parent:3	Parent specific storage pilot point multipliers: layer 4	log	80	1.00	0.70	100.00
_ss_pp_parent:4	Parent specific storage pilot point multipliers: layer 5	log	80	1.00	0.70	100.00
_sy_con_inset__multiplier	Inset specific yield multipliers by layer	log	5	1.00	0.75	2.00
_sy_con_parent__multiplier	Parent specific yield multipliers by layer	log	5	1.00	0.75	2.00
_sy_pp_inset:0	Inset Specific Yield pilot point multipliers: layer 1	log	264	1.00	0.75	2.00
_sy_pp_inset:1	Inset specific yield pilot point multipliers: layer 2	log	264	1.00	0.75	2.00
_sy_pp_inset:2	Inset specific yield pilot point multipliers: layer 3	log	264	1.00	0.75	2.00
_sy_pp_inset:3	Inset specific yield pilot point multipliers: layer 4	log	264	1.00	0.75	2.00
_sy_pp_inset:4	Inset specific yield pilot point multipliers: layer 5	log	264	1.00	0.75	2.00
_sy_pp_parent:0	Parent specific yield pilot point multipliers: layer 1	log	80	1.00	0.75	2.00
_sy_pp_parent:1	Parent specific yield pilot point multipliers: layer 2	log	80	1.00	0.75	2.00
_sy_pp_parent:2	Parent specific yield pilot point multipliers: layer 3	log	80	1.00	0.75	2.00
_sy_pp_parent:3	Parent specific yield pilot point multipliers: layer 4	log	80	1.00	0.75	2.00
_sy_pp_parent:4	Parent specific yield pilot point multipliers: layer 5	log	80	1.00	0.75	2.00
lak_mults	Lake evaporation/precipitation multipliers by stress period	log	680	1.00	0.90	1.10
lakebedk	Lakebed conductivity values by littoral or profundal zone	log	8	0.02 to 0.04	1.00e-04 to 1.00e-02	0.10 to 5.00
rch_child__multiplier	Inset recharge multipliers by stress period	log	85	1.00	0.80	1.20
rch_parent__multiplier	Parent recharge multipliers by stress period	log	85	1.00	0.80	1.20
sfrk	SFR conductance values by reach	log	15	1.00	0.50	2.00
wel__parent_cnst	Well pumping multipliers by stress period	log	85	1.00	0.80	1.20
wel_inset_cnst	Well pumping multipliers by stress period	log	85	1.00	0.80	1.20

History-Matching Results

The ensemble approach in iES results in an ensemble of objective function values from iteration to iteration for the Pleasant Lake and Plainfield Tunnel Channel Lakes inset models (figs. 23 and 24, respectively). The blue curve represents the base ensemble realization, and the lighter gray lines indicate the trajectories of all the other ensemble realizations (figs. 23 and 24). A wide range of parameter combinations can be examined with this approach, and some of those combinations result in MODFLOW-6 models that fail to run to completion. Others result in poor correspondence between model outputs and observations, expressed as high objective function values. As a result, two subjective steps considered here are required to settle on a final posterior distribution of parameters and model outputs. First, an iteration is chosen from the objective function progress. This selection is made based on a heuristic balance between level of fit and variability in the ensemble. At later iterations, overfitting and collapse of the ensemble are common; therefore, an early iteration typically is chosen. In both inset models, the third iteration was chosen as an acceptable level of fit. Second, a histogram of objective function values is evaluated at the chosen iteration and a cutoff of the objective function is chosen to remove outlier ensemble realizations that result in an unacceptably high objective function value—this process is referred to as “rejection sampling.” The histograms, cutoff objective function values, and effects on ensemble size for the Pleasant and Plainfield Tunnel Channel Lakes inset models are shown on figures 25 and 26, respectively.

Comparison of Measured and Simulated Observations—By Observation Group

The match between model-simulated and measured observation values can be evaluated for each observation group as a 1:1 line plot and as residuals. Residuals are calculated by subtracting the model-simulated value from the measured value for each observation. As described previously for the regional model, values plotting closer to the 1:1 line indicate a better match between model-simulated and measured values. The residuals plots provide a means to assess prediction bias in the history matching results over the entire range of observation values. The uncertainty of both the observations and the ensembles of model outputs for those observations is depicted as +/−2xσ, where σ is calculated as the inverse of the observation weights and is calculated empirically from the ensembles for model outputs. This is a general metric of uncertainty in how much information is expected to be provided by each observation and how variable the resulting parameter estimate ensembles are. The fit and residuals plots for the Pleasant Lake inset model are shown in figures 27–34, and the fit and residuals plots for the Plainfield Tunnel Channel Lakes inset model are shown in figures 35–39. The observation group names are defined in tables 2 and 3.

Fit Between Measured and Model-Simulated Observations—Spatial Distributions

The spatial patterns of weighted residuals for observations used for history matching in the Pleasant Lake and Plainfield Tunnel Channel Lakes inset models are shown in figures 40 and 41, respectively. These results were derived only from the base-case parameter set obtained from the iES ensembles, and, therefore, do not reflect the range of simulation results available in the ensembles. Most of the observations represent transient conditions; therefore, the mean of the base realization over the entire simulation time is presented. Only a partial understanding of the history matching process can be determined from figures 40 and 41, but these figures are useful for examining spatial patterns. Residuals showed spatial bias; for example, water levels near the lake were simulated greater than measured values, but the magnitude of residuals was small. In general, simulated values that were greater than measured values were closer in fit than when simulated values were less than measured values. The streamflow observations showed little spatial bias with a mix of simulated values greater and less than measured values. Subsequent sections describe the relation of simulated and measured values at specific locations in time series are discussed in detail below.

Fit Between Measured and Simulated Observations–Time Series by Observation Location

Individual observation locations with transient data (hydraulic head and streamflow) provide multiple targets for history matching yet may result in multiple fit values for each location, spanning over model-simulation time. The uncertainty of the observed values, as determined for the iES process as the inverse of their observation weight, can be expressed as a 95-percent credible interval using bars that extend two times σ above and below the observed value at a given time. Unlike the 1:1 plots (figs. 38 and 39), an accompanying bar for the simulated values is not in figure 42. Instead, the time series for each ensemble realization is displayed separately as a thin gray line. This time series indicates not only the uncertainty of both observations and simulated values but also the continuity of values over time in each ensemble realization.

Observation points labeled as prior data conflict (PDC) indicate observations that were under prior data conflict, meaning that, in these cases, the initial ensemble evaluated by iES identifies observations where simulated values are statistically distant (White and others, 2021) from the measured values. Because iES evaluates a wide range of parameters in this initial run, the history matching process is unlikely to result in a reasonable set of parameters that will make the model fit close enough to be considered valid. In other words, the model is being used to flag observations that cannot be fit with the indicated model parameters. This can result from observation errors or from deficiencies in the model conceptualization or validity of model parameters, or both. In PEST++, these observations may be systematically assigned zero observation weight to prevent them from affecting the history matching process. Alternatively, PDC can serve as a diagnostic tool to point out areas of concern that would benefit from additional study to improve the quality of the observations, model conceptualization, and/or model parameters. In this study, the PDC labels were assigned to observations in prior data conflict to indicate that observations were unlikely to result in low residuals. These observations were retained with their initially assigned observation weight. Results for time series in each of the two inset models are described below.

Pleasant Lake Inset Model

The general patterns of observed lake levels are reproduced in the Pleasant Lake inset model and, at most times, the simulated results from the model realizations overlap with the measured values (fig. 42). The bars indicating uncertainty of the measured values are not visible in these plots because the observation weight was assigned high enough to these important observations that the implied uncertainty is small. The model results at the beginning of the simulation period may be affected by a need for a “spin-up” or “warm start” for the simulation of the lake level to reach a dynamic equilibrium with respect to changes in groundwater storage. To account for this potential discrepancy between simulated and measured lake levels, recent observations were assigned higher observation weight than values for time steps early in the model-simulation run that include data from early in the simulation period. The other period of large differences, between simulated and measured values, occured in the spring of 2018. This result may indicate that the lake level was anomalously high because of above-average surface runoff associated with the spring melt that year. Surface runoff is not simulated in this model, so under-simulation of the lake level under that above-average runoff condition is expected.

The patterns and magnitudes of simulated streamflow were generally similar to the measured values throughout the Pleasant Lake inset model (figs. 43–46). Various observations are labeled as PDC. These PDC observations are typically locations where either the measured streamflow was already low and dropped to near zero or represented as high flow because of runoff. Whereas the streamflows in the study area generally are dominated by base flow, some storm events do occur and the SFR package in MODFLOW 6 does not include explicit consideration of storm runoff.

The time-series results for the piezometers near Pleasant Lake in the wdnrwells_tr group—wells that represent hydraulic head conditions around the lake perimeter—are shown in figures 46–47. These model results are nearly all within the 95-percent credible range for the observations. The observation group names are defined in tables 2 and 3.

The groundwater hydraulic-head values (water levels) in the wells of the observation group nwisdvs_tr for which daily values are available in the National Water Information System database (U.S. Geological Survey, 2021) are shown in figures 48 and 49. Water levels measured in these wells reflect the overall pattern of groundwater levels in recent time. Where the greatest discrepancies are observed (particularly in fig. 49), the simulated values typically are too high rather than too low. This bias may be, at least in part, an unintended consequence caused by the assignment of high observation weight on streamflow and lake-level observations that results in a parameter set yielding higher water levels in the aquifer to better match those observations. For streamflow, in particular, this parameter could be compensating for water that enters the streams through processes not explicitly simulated in these models.

Plainfield Tunnel Channel Lakes Inset Model

The lake-level time series for the Plainfield Tunnel Channel Lakes inset model are shown in figure 50. The recent (2018) and long-term patterns (2017–18) for Plainfield Lake are shown in panels A and B, respectively. The recent (2018) and long-term (2013–18) patterns for Long Lake are shown in panels C and D respectively. The general patterns of lake levels are reproduced, and at most times the simulated realizations overlap with the observed (measured) values. Unlike Pleasant Lake, where some appreciable differences likely resulted from spring runoff, in Plainfield and Long Lakes, the overall dynamics are visible in the most recent monthly data as shown in panels A and C, respectively.

The time series of streamflow in the two streams (North Branch Tenmile Creek and Upper Pine River) with history matching data in the Plainfield Tunnel Channel Lakes inset model are shown in figure 51. The patterns and magnitudes of simulated streamflow were generally similar to the measured values. Nearly all ensemble simulated values fall within the 95-percent credible intervals for the streams. Streamflow in North Branch Tenmile Creek is simulated as higher values than those measured at the extremely low observed streamflow conditions between 2012 and 2015. Such low values are difficult to represent accurately with the model and difficult to quantify when measured in the field, as indicated by the relatively large credible interval around the observations. Similar to some of the observations in the Pleasant Lake model, at least one storm event in late 2017 simulated too low in both streams, likely because of runoff not being explicitly represented in the model. The simulated results in the Upper Pine River are closer to measured values with the exception of the late 2017 storm event.

The time-series results for the piezometers near Long Lake in the wdnrwells_tr observation group are shown in figures 52 and 53; the time-series results for the piezometers near Plainfield Lake in the wdnrwells_tr observation group are shown in figures 54 and 55. These wells represent hydraulic-head conditions around the lake perimeter, and the simulated results are nearly all within the 95-percent credible range for the observations.

Time series for groundwater hydraulic-head values in the wells of the nwisdvs_tr group, for which daily values are available in the USGS National Water Information System database, are shown in figure 56. Site 441146089250301, which includes observations from the beginning of the simulation, is shown in panel A. Time series for the wells shown in panels B, C, and D represent recent time (2018–19). Similar to the Pleasant Lake inset model, where the fit is poorest, simulated values generally are higher than measured values. An exception is the winter of 2017–18, where simulated values briefly decreased while measured values were increasing at site 441146089250301 (fig. 56A); however, the overall water-level response is similar in increasing water-level elevations at the end of the simulation period.

Focus Area Inset Model Results

After history matching, both inset models were updated with the “base” parameter set from iES iteration number 3, which represents the single best or “minimum error variance” estimate of parameter values (see the “Parameter Estimation and Uncertainty Analysis” section). The updated models were then run for the history matching period from 2012 through 2018 and used for the scenarios described in the “Lake Ecosystem Response Assessment (LERA)section. Results for the 2012 through 2018 history matching period are described below.

Parameter Estimates

The base parameter set spatial distributions for the inset models, after history matching, are shown in figures 57–64. The use of multipliers against the parameter values estimated at the regional scale resulted in parameter fields similar to those values in the regional model. In the Pleasant Lake model, in particular, some greater changes in parameter values were estimated in the inset models focused around the lake relative to the parent model region. The high observation weight assigned to lake elevation observations and the presence of headwater streams that are highly sensitive to small changes in parameters make the parameters in the inset models particularly critical for accurate model simulation and subject to greater changes than the surrounding area.

Groundwater Flow and Boundary-Condition Fluxes for Mean (2012–15) Conditions

Simulated steady-state water-table elevations and boundary-condition fluxes for both inset models, under mean conditions for the period of 2012–15, are shown in figures 65 and 66. Discharge from the groundwater-flow system to streams, wells, or across the specified-head perimeter boundaries out of the model domain is represented with a red color scale. A blue color scale represents discharge to the simulated groundwater-flow system from stream leakage or inflow into the model domain through the specified-head boundaries. Note that for streams, this color scheme is opposite of that shown in figure 12 for the regional parent model, where blue colors indicate groundwater discharge to streams (inflows from the perspective of the streams). Lake/groundwater interactions are displayed separately from the perspective of the lakes, with blue indicating groundwater discharge to the lakes and red indicating lake leakage to groundwater.

The overall flow patterns indicated by the simulated water-table elevation contours are consistent with the regional model (fig. 16). In both inset models, groundwater discharge along the perimeter is simulated mostly as leaving the model domains, a result of the model domains being located along the groundwater divide that runs through the streamless portion of the study area. A small amount of groundwater flow enters the domains from the north, consistent with regional flow patterns. Similarly, streams within the inset models are mostly gaining groundwater, with the exception of the Upper Pine River (fig. 12) and Tagatz Creek (fig. 66), which are in the hummocky collapsed terrain east of the moraines. Simulated groundwater discharge to and from the study lakes reflects the simulated water-table gradients, with the lakes gaining groundwater on the upgradient side and losing water to the groundwater system on the downgradient side.

Groundwater Flow and Boundary-Condition Fluxes for Mean (2012–15) Conditions

The inset submodel budgets are dominated by flows exchanged with the parent submodels (“FLOW-JA-FACE” term). This result is acceptable because the two submodels are fully coupled at the inner iteration level, meaning continuity in the groundwater-flow model final solution is preserved throughout the inset model domain. Otherwise, the inset submodel water budgets show some of the same seasonal trends as the parent submodel. Also, the inset submodel budget indicates a small magnitude of flux through the lakes relative to other sources and sinks.

The simulated water budgets for Plainfield and Long Lakes (fig. 71) show a seasonal pattern of mostly net groundwater outflow (out of the model domain) in the fall and winter months, with net groundwater inflow (into the model domain) in the spring, when groundwater levels are high. In both lakes, there is a pattern towards predominantly gaining conditions later in the simulation, with net groundwater inflow for most of 2018, consistent with the field observations by WDNR (WDNR, 2021, app. A). In contrast, Pleasant Lake has a relatively stationary seasonal cycle of predominantly positive net groundwater outflow, presumably because of the stabilizing effects of nearby Chaffee and Tagatz Creeks on the water table.

Lake Ecosystem Response Assessment (LERA)

One of the goals of this study was to evaluate the potential local-scale effects of agricultural water and land-use practices on lakes as well as on the regional groundwater-flow system. To accomplish this goal, WDNR generated two potential land-use scenarios to be assessed with the groundwater-flow models (WDNR, 2021, app. F) for details). One land-use scenario, referred to as the “current irrigated agriculture” scenario, is based on current (2018) land-use and irrigation practices (fig. 72). This scenario can be considered a “control” or “business as usual” case that reflects the best estimate of hydrologic responses given current conditions. A second land-use scenario, referred to as the “no irrigated agriculture” scenario, is one in which present (2018) irrigated agriculture areas are replaced with a land cover consistent with slope, soil-type, and other factors of the surrounding areas (fig. 73) (WDNR, 2021, app. F). This alternative land use is meant to represent a land-use pattern if irrigated agriculture had not been introduced to the Central Sands region. As a result, unlike land use in irrigated agricultural fields, where crop type varies from year to year, the land use for this alternative no irrigated agriculture scenario is static over time. In areas of the model domain where irrigated agricultural land use is not simulated for 2018 conditions, the two scenarios are the same. The no irrigated agriculture scenario is meant to provide a comparison point, which compared with the current irrigation scenario, helps identify potential effects from irrigated agriculture. These two land-use scenarios were run for a 38-year simulation time with the SWB code at the regional model scale to provide recharge and water-use values as needed for input into the regional model. The results from the regional model simulations then provided specified-head boundary conditions for the two focused lake inset models. The recharge and water use from the SWB simulations (Westenbroek and Fienen, 2022) were also assigned to the two inset models to provide lake-level and base-flow model outputs at the same locations as the history matching data. This analysis enabled a comparison of long-term responses in lakes and streams to land/water use with and without irrigated agriculture.

The 38-year period (climate years 1981–2018) is the timeframe for which daily PRISM precipitation and air-temperature data were available, both of which are necessary inputs to SWB and the MODFLOW-6 LAK package. This is referred to as “climate years” to highlight that only the climate signal is used for that time period—land use is defined separately, and the scenarios are not meant to simulate overall conditions in the specific climate years. Furthermore, this period included multiple cycles of the typical range of weather variability seen in the region (see WDNR 2021, app. B for additional discussion).

The land-use code for a given SWB model cell defines the kind of vegetation in the land parcel the cell is associated with. This code controls how evapotranspiration is calculated in SWB. For the no irrigated agriculture scenario, the vegetation assemblage was assumed to remain static in each cell over the 38-year LERA simulations.

The current (2018) irrigated agriculture scenario is dynamic through time, as most agriculture in the region is, and involves rotations of multiple crops. To account for this land-use variability in SWB, each parcel was assigned a crop rotation (WDNR, 2021, app. F) based on four categories that describe prevailing agricultural practices. The distribution of the main crop rotations that include “potato/vegetable,” “cash grain,” “dairy,” and “pasture/hay” is shown in figure 74. Non-agricultural land use is quantified as land that is owned by farms but not used for growing crops. Within each crop rotation, the current distribution (2018) of land-use codes was evaluated resulting in a discrete distribution of land-use codes for all parcels in the model domain assigned a rotation. The rotation determines which assemblage of crops are rotated through that location over time, and the discrete distributions within each rotation are used to assign specific crops for each year. To assign land use in each climate year, the full complement of parcels assigned a rotation was filled with a proportion of land-use codes corresponding to the current distributions as shown in figure 75.

The spatial distributions of the crop rotations and the alternative no irrigated agriculture scenario land use assumed in the same parcels for the regional model domain are shown in figures 72 and 73, respectively. The spatial distributions of crop rotations and the alternative no irrigated agriculture scenario land use assumed in the same parcels for the Plainfield Tunnel Channel Lakes and Pleasant Lake inset models are shown in figures 76 and 77, respectively. The intermorainal areas, where the study lakes are located, are areas where the “potato/vegetable” rotation is more prevalent relative to the other rotations more commonly in the outlying parts of the study area (figs. 76 and 77). For the no irrigated agriculture scenario, the alternate land use is assumed to be dominated by forest land use (WDNR, 2021, app. F).

For the purpose of this analysis, special consideration was given to potatoes within the “potato/vegetable” rotation parcels because potatoes (land-use code 43) account for 25 percent of total land in the “potato/vegetable” rotation (fig. 75). To account for typical growing practices, a 4-year recurrence of potatoes was enforced so that every parcel in the “potato/vegetable” rotation was simulated as growing potatoes once every 4 years, while still maintaining that 25 percent of all the “potato/vegetable” rotation be simulated as growing potatoes. After this requirement was met, the remaining parcels were randomly assigned crops following the distribution shown in figure 75.

The land-use files for the no irrigated agriculture scenario and the current irrigated agriculture scenario were provided as input to SWB, along with representative climate data from PRISM resulting in grids of agricultural demand and net infiltration, aggregated to monthly stress periods. Using parcel mapping from WDNR (WDNR, 2021, app. F), water estimated by the agricultural demand was aggregated to permitted irrigation wells. The resulting recharge and well information was then input to the regional model for simulation over the climate years to provide specified-head boundary conditions for the inset models. The inset models were then run using the boundaries available from the updated regional model, the recharge and well pumping data from the SWB runs, and evapotranspiration and precipitation input for the LAK package, updated using the representative PRISM data. The other parameters for the three MODFLOW models were obtained from history matching and using the optimal parameters for the regional model and the base realizations for the two inset models. To focus on differences between the no irrigated agriculture and current (2018) irrigated agriculture scenarios, the same specified-head boundaries from the regional model (the no irrigated agriculture results) were used for all inset model LERA runs.

Three-hundred fifty SWB parameter sets were generated using a Monte Carlo approach. These parameters are assumed to vary about the values assigned during the base SWB run. The parameters that were allowed to vary include many parameters related to the FAO-56 calculations: irrigation starting and ending dates, planting dates, crop-coefficient curve midpoint (Kcb_mid), the maximum allowable depletion (responsible for triggering virtual irrigation events), and the depletion fraction defining soil-moisture conditions leading to plant stress. The plant rooting depth was a non-FAQ-56 related parameter allowed to vary.

Each of the parameter sets discussed above was simulated with SWB, and the outputs were then used as input to the inset models to examine the responses and associated uncertainties of lake levels and streamflows derived from varying SWB input parameters. The simulated lake levels for all the study lakes are shown in figure 78. The darker blue and orange lines indicate the results for the base run parameter set for the Monte Carlo SWB runs for the no irrigation and current irrigation scenarios, respectively. Lighter blue lines show the ensemble of results resulting from sampling the uncertainty of the SWB parameters. Note that in panel B, short periods indicate that the lake level was simulated at the bottom of the lake bathymetry for some realizations. Similar plots for streamflow in two selected streams from the Plainfield Tunnel Channel Lakes inset model are shown in figure 79. Streamflow results for two sets of selected streams in the Pleasant Lake inset model are shown in figures 80 and 81.

A final LERA analysis run was a potential irrigated agriculture scenario where a land-use scenario was generated (WDNR, 2021, app. F). In this scenario, all land considered suitable for irrigated agriculture was simulated as irrigated agriculture. The land-use patterns were assigned rotations and annual crop distributions following the same methodology as the current irrigation scenario discussed above. Only the base SWB parameters were simulated for this sensitivity run and the lake stages for the base SWB parameters for the no irrigated agriculture, current irrigated agriculture, and potential irrigated agriculture scenarios are shown in figure 82.

Two final evaluations following the same conceptualization as the LERA runs were performed to more precisely evaluate the connection between irrigation and land use at individual parcels and the responses in lake levels. These were the “distance” and the “lag-time” evaluations, both of which were limited to Long Lake as a proof-of-concept. The distance from the geometric centroid of Long Lake to each permitted high-capacity well was used to rank wells from 0 (closest) to 319 (farthest). Using this information, the two scenarios evaluated (1) converted each well, cumulatively, from wells 0 to 319, in sequence, from current irrigated agriculture to no irrigated agriculture (the distance scenario) and (2) divided the 320 wells into 20 groups of 16 wells and converted all wells in each group from current irrigated agriculture to no irrigated agriculture (the lag-time scenario). The distance scenario illustrates the cumulative effects of pumping on lake levels as a function of distance. The lag-time scenario illustrates the magnitude and timing of effects from smaller well groups assigned based on similar distance on lake levels. The conversion from current irrigated agriculture to no irrigated agriculture was simulated for both scenarios by (1) removing the high-capacity well(s) from the simulation and (2) substituting the recharge calculated for the no irrigated agriculture in place of the recharge calculated for the current irrigated agriculture scenario for the parcel(s) associated with the well(s). This is the same conceptualization as the LERA runs above, but with subsets of wells selected for evaluation.

Results of the distance scenario are presented and discussed further by the WDNR (WDNR, 2021, app. B). For the lag-time scenario, the 20 groups of 16 wells each, grouped as quantiles of distance rank from the center of Long Lake, are shown in figure 83. The colors indicate the total amount of water simulated as pumped for irrigation in each group. The differences, over time, in the level of Long Lake over the 38-year LERA simulation time, simulating the wells in each group using the no irrigated agriculture recharge and pumping scenario relative to the current irrigated agriculture recharge and pumping scenario, are shown in figure 84. Positive values indicate higher lake levels in the no irrigated agriculture scenario relative to the current irrigated agriculture scenario.