Hydrogeologic Units as Model Layers

The hydrogeologic framework consists of the six hydrogeologic units described in the conceptual model: (1) surficial materials, (2) Des Moines River alluvium, (3) glacial outwash, (4) glacial till, (5) buried valley alluvium, and (6) bedrock. These hydrogeologic units were incorporated as distinct layers in the hydrogeologic framework and numerical model. Hydrologic parameter zones were assigned to each layer; the glacial outwash, glacial till, and bedrock layers contain only one zone each, whereas surficial materials, alluvium, and buried valley alluvium layers were subdivided into multiple hydrogeologic zones based on compositional or landform variations within the layer.

Layer 1: Land-Surface Material

The top model layer (layer 1) represents the hydrogeologic materials that exist at the land surface in the model area (fig. 4). The top elevation of this layer was determined using the average land surface elevation over the 60-m model cell size from a digital elevation model (DEM; Iowa Geospatial Data, 2022) The bottom elevation for layer 1 was set to the top of the underlying units. Median thickness for the surface layer is 2.0 m, with a range from 1.0 m to 53.0 m. In most areas of the model, surficial deposits are conceptualized as 2.0 m thick. In general, the cells in layer 1 with high thickness are cells representing Saylorville Dam (where the full height of the dam down to bedrock is included in layer 1) and in some surface stream cells (where layer 1 was thickened to accommodate a reasonable stream gradient along streamflow paths). Layer 1 was delineated into four hydrogeologic zones based on general soil permeability classifications: alluvial and outwash-based soils (zone 1), till-based soils (zone 2), areas of fill (zone 3), and quarries (zone 4). Although quarries were delineated as a separate hydrogeologic zone within layer 1, they were ultimately represented in the model using a separate boundary condition package and so zone 4 was not included in model parameterization.

Layer 2: Des Moines River Alluvium

Layer 2 is the top layer of alluvium in the Des Moines River alluvial plain and consists primarily of Holocene-age sediments that were deposited by the Des Moines River (fig. 5). This unit was subdivided into three hydrogeologic zones: low (zone 5), intermediate (zone 6), and high (zone 7) alluvial terraces identified as the DeForest Formation in Quade and others (2003). DeForest Formation alluvial terraces are described as 1ؘ–5 or 7 m of silty clay loam, clay loam, or loam overlying the Noah Creek formation (Quade and others, 2003; glacial outwash sand and gravel). This material was delineated within the Des Moines River alluvial plain only. Beaver Creek also has associated alluvium deposits, but in this model, these are included in the surface material layer because well lithology logs indicate that the Beaver Creek alluvium is relatively thinner than the alluvium associated with the Des Moines River. This unit is part of the Des Moines River alluvial aquifer, along with underlying older sands and gravels.

The surface elevation of this unit was set to 3 m beneath a smoothed land surface DEM sampled to a 60-m raster grid, with a median elevation of 242 m. In the upstream part of the Des Moines River alluvial plain within the model, where the valley is relatively wide, alluvium overlies older glacial outwash sand and gravel and, in some areas, buried valley alluvium. Elevation data for these lower units are sparse, so the alluvium unit bottom elevation was fixed to 234 m in most cells and higher only where necessary to accommodate higher bedrock elevation.

Layer 3: Glacial Outwash

Layer 3 in the model represents glacial outwash deposits in the Beaver Creek valley and buried valley outside of the Des Moines River valley. Glacial outwash in this model is identified as the Noah Creek formation, described as 3–15 m of coarse-grained sand and gravel overlying Dows formation glacial till or bedrock (fig. 6; Quade and others, 2003). In the glacial outwash channel and Beaver Creek valley, this unit overlies a layer of fine-grained massive glacial till. Glacial outwash also is present in the river valley, where outwash overlies either bedrock or buried valley alluvium. For this model, glacial outwash in the Des Moines River valley has been included with buried valley alluvium in layer 5. The glacial outwash and buried valley alluvium are distinct units with separate depositional histories, and outside of the river valley they are separated by a layer of glacial till; however, within the Des Moines River Valley these units are in direct contact and form a single aquifer (along with the overlying more modern alluvium). Additionally, glacial outwash in the Des Moines River valley has greater thickness and depth than the glacial outwash assigned to the Beaver Creek valley and the buried channel outside the river valley, and the glacial outwash in the river valley is horizontally and vertically disconnected from these other outwash deposits. Glacial outwash within the Des Moines River valley was therefore modeled as a single unit with underlying buried valley alluvium (layer 5) to better represent the layer connectivity in the numerical model.

Outside of the Des Moines River valley, the outwash surface elevation was set to 2 m below a smoothed land surface elevation grid, and the base elevation was set to the surface elevation of the underlying till unit. The outwash was pinched out in areas where shallow depth to bedrock prevented at least 1 m of unit thickness. Layer 3 contains only one hydrogeologic zone (zone 8).

Layer 4: Glacial Till

Layer 4 in the model represents Wisconsin-age glacial till that underlies uplands surrounding the Des Moines River valley and parts of lowland areas (Quade and others, 2003). Quade and others (2003) identified this as the Dows formation and describes different Dows formation subunits within the bounds of this model, each consisting of surface sediments overlying massive, dense loam diamicton (till). For the purposes of this model, these units are combined into one homogenous glacial till material, and the associated surficial sediments form part of the surface layer of the model (fig. 7).

Glacial till is present throughout the model except in the Des Moines River alluvial plain and in places where shallow depth to bedrock precludes adequate accommodation space for till. In the upland areas, glacial till top elevation was set to 2 m below a smoothed surface DEM and adjusted to maintain at least 2 m for the overlying surficial materials model layer. Till bottom elevation in uplands was set to the bedrock surface, and the unit was set to zero thickness in any cell where it was not possible to have at least 1 m of thickness. In the Beaver Creek valley and the buried valley areas, the surface elevation of the till unit was initially fixed at 247 m then adjusted to be no more than 8 m or less than 4 m lower than the smoothed land surface. As with upland cells, till thickness was set to zero in cells without accommodation space for at least 1 m of till. Layer 4 contains only one hydrogeologic zone (zone 9).

Layer 5: Buried Valley Alluvium

Layer 5 represents a layer of pre-Wisconsinan alluvial sand and gravel at the bottom of the Quaternary sequence, overlying the shale bedrock of the buried channel that underlies Beaver Creek valley in the study area and extends into the southeast part of the model area (fig. 8). This buried channel and its sand and gravel fill have been noted in a geologic study (Bettis and Hoyer, 1986) and in well lithologic logs from the area (IGS, 2022). In most places, this buried channel deposit is overlain by glacial till and another sequence of sand and gravel outwash (the Dows and Noah Creek formations); however, where the buried channel deposit crosses the Des Moines River valley, the glacial till layer is absent and this buried channel alluvium directly contacts overlying outwash sand and gravel or more modern alluvium.

In the buried valley where this alluvium is present, its top surface elevation was fixed at 237 m. This is consistent with lithologic log data where the alluvium has been noted (IGS, 2022). Bottom elevation for this unit was set to the bedrock surface elevation. The thickness of this unit in the model ranges from 1 to 17 m, and the layer is precluded in areas within the valleys where bedrock top elevation is higher than 236 m. This material also underlies part of the Des Moines River alluvial plain, where it forms an aquifer with overlying glacial outwash and more modern alluvium. Within the river valley, the surface elevation of the combined glacial outwash and buried valley alluvium unit was fixed at 234 m, its base elevation was set to the bedrock surface elevation, and the unit was pinched out in areas at the valley margin with shallow depth to bedrock. These values were all generalized from existing well logs, which are sparse in some areas. Layer 5 was subdivided into four zones: buried valley alluvium (zone 10); combined buried valley alluvium and glacial outwash in areas of the river valley where they are both present (zone 11); glacial outwash in the Des Moines River valley (zone 12); and a small region of buried valley fill underlying Rock Creek (zone 13, fig. 1).

Layer 6: Bedrock

Layer 6 in the model represents shale of the Cherokee Group that underlies the model area (fig. 9). The elevation of the bedrock surface was determined using a published bedrock elevation map (Witzke and others, 2010), lithologic logs from the GeoSam database (IGS, 2022), the published bedrock surface from the 2021 Des Moines River groundwater model (FitzGerald and others, 2022), and geophysical data collected in the model area by the USGS (Gruhn and others, 2021). For this model, bedrock is assumed to be a single homogenous unit with low hydraulic conductivity. Bottom elevation of this unit was set to 100 m above North American Vertical Datum of 1988 in every active cell. Layer 6 contains just one hydrogeologic parameterization zone (zone 14).

Water Budget Components

Water budget components for this model include sources of recharge to groundwater and sources of groundwater loss. In this model, sources of recharge to groundwater include seepage from streams, infiltration of precipitation, and groundwater inflow from outside the model boundary. Sources of groundwater loss include discharge to streams and other surface water bodies including quarries, ET_g, withdrawals from wells, and groundwater flow outward through the model area boundary. Recharge and evapotranspiration estimates were generated using a previously published Precipitation Modelling Runoff System (PRMS) model of the Des Moines River Basin (Haj and others, 2015). Streamflow data used to estimate streambed seepage (streambed flux) were assembled from records of three USGS streamgages in the model area with continuous data for all or part of the model period (USGS, 2023a, b, c; fig. 2). Groundwater withdrawal rates were compiled from records furnished by DMWW (Dylan Hart, Black & Veatch Corporation, written commun., 2023; Davis and Bristow, 2024) and records retrieved from the State of Iowa Water Allocation Compliance and Online Permitting system (Iowa Department of Natural Resources, 2023).

Streamflow and Its Interaction with the Groundwater System

Streams in the model area are hydrologically connected to the underlying aquifers and contribute water to and drain water from the groundwater system. Streamflow gains and losses were calculated by subtracting upstream discharge at the Beaver Creek and Saylorville gages from the downstream discharge at the 2nd Avenue gage (fig. 10). The stretch of Beaver Creek and the Des Moines River between these three gages gained an average of 306,000 cubic meters (m³) per month in combined groundwater discharge and overland runoff from October 1986 through July 2022. Streamflow gain happened in 225 of the 310 months of this period and streamflow loss, because of recharge to the aquifer from river seepage or evapotranspiration, happened in 85 of 310 months (fig. 10A). The average monthly streamflow gain is highest in April through July (fig. 10B). In winter months (October through March), average monthly streamflow gain is 132,000 m³.

The Saylorville gage is approximately 3 km downstream from the Saylorville Dam outlet structure. One important ungaged creek, Rock Creek, flows into the Des Moines River between the dam outlet and the gage (fig. 1). Input streamflow for the MODFLOW model Streamflow Routing (SFR) Package for the Des Moines River below the Saylorville Dam was calculated by subtracting the PRMS-calculated streamflow for Rock Creek from the observed Des Moines River streamflow at the Saylorville gage.

Lateral and Lower Boundaries

Most of the active model area boundary is a no-flow boundary following topographic highs underlain by glacial till or bedrock (model layers 4 and 6), with the assumption that the groundwater drainage basin is coincident with the surface water drainage basin in the model area. The base of the model also is represented as a no-flow boundary. Lateral no-flow model boundaries are established by deactivating cells outside the boundary.

Groundwater Inflow and Outflow Boundaries

Aquifers in the model intersect the model boundary in three locations: the Des Moines River alluvial aquifer intersects the model boundary at the southeast edge of the model where the river exits the model area; and the buried valley aquifer intersects the model boundary in the northwest and southeast. Groundwater movement in and out of the model for these three boundaries was represented in the model using the MODFLOW 6 GHB Package (Langevin and others, 2017). The downstream alluvial aquifer GHB boundary was assigned to model layer 2 and head for each cell was assigned a constant value based on the average Des Moines River stage at the 2nd Avenue gage plus an elevation offset ranging from 0 to 0.79 m to approximate the floodplain topography and maintain a water table gradient toward the river. The GHB boundaries representing the buried valley alluvium were assigned to layer 5, the lowest aquifer in those cells. Hydraulic head for the northwest boundary was assigned a constant value equal to the approximate surface elevation of Beaver Creek in that area (255 m), plus an elevation offset ranging from 0 to 0.6 m to make a slight gradient in the water table toward Beaver Creek. Hydraulic head for the southeast buried valley boundary was fixed at 246.5 m for all cells associated with the southeast buried valley boundary.

Recharge from Precipitation

Recharge from precipitation infiltration into the aquifers was simulated using the MODFLOW 6 RCH Package. The RCH Package is designed to simulate aerially distributed recharge to the groundwater system (Langevin and others, 2017). Recharge was specified for each cell for each monthly stress period in units of meters per day and was applied to the top cell throughout the model. Recharge rates were calculated using the PRMS model (as described in the “Water Budget Components” section) for each PRMS HRU polygon and then each grid cell was assigned the recharge rate of the nearest PRMS HRU. An additional multiplier based on surficial geology was applied for each cell, as described in the “Water Budget Components” section. These combined spatially distributed rates were used as initial estimates of recharge to the numerical model and included as parameters for calibration.

Surface Water Interaction

Streams in the model were simulated using the MODFLOW 6 SFR Package, which calculates streamflow and streambed seepage (streambed flux) using the continuity equation. SFR calculations assume piecewise steady (nonchanging in discrete time periods), uniform (nonchanging in location), and constant-density streamflow, such that volumetric inflow and outflow rates are always equal, and no water is added to or removed from storage in the surface channels. The numerical model routes streamflow through a network of rectangular channels (which may include rivers, streams, canals, and ditches and are referred to collectively as reaches in the remainder of this report; Langevin and others, 2017). Reaches represent a stream associated with a particular cell in the numerical groundwater flow model. A stream on the landscape is represented in the numerical model by reaches delineated within individual model cells, numbered sequentially from upstream to downstream, and assigned uniform physical and hydrologic properties. Rates of overland flow, precipitation, and evaporation also can be assigned for each reach in the SFR Package; however, these rates are not assigned for this model. In this report, segments are conceptually defined as groups of one or more individual reaches between stream confluences. Streambed hydraulic conductivity values were assigned uniformly to the reaches that make up stream segments, except for segments of the Des Moines River, which had values of streambed hydraulic conductivity uniquely assigned to reaches based on CRP surveys, discussed later in this section.

Streams including the Des Moines River, Beaver Creek, and perennial tributary creeks at least 1-m wide were included in this model as sequentially numbered SFR reaches, with each reach corresponding to a single model cell as described in the previous paragraph. All SFR reaches were assigned to layer 1, the top-most layer in each cell. Each reach was assigned physical characteristics of length, width, gradient, top elevation of streambed, streambed thickness, Manning’s roughness for the channel, and streambed hydraulic conductivity. Length for each reach was calculated by intersecting each stream channel center line with model grid cells, and width was approximated from aerial imagery (Google Earth, 2022). Streambed top elevations for stream segment endpoints were determined using a 3-m resolution DEM and interpolated to individual reaches along each stream segment (Iowa Geospatial Data, 2022). Bed thickness for each reach was nominally fixed at 1 m. In cases where the bottom of the streambed was below the bottom of layer 1, the elevation of the bottom of layer 1 was adjusted to 1 m below the elevation of the bottom of the streambed to accommodate the reach. Manning’s roughness coefficient (n) was fixed at 0.037 for all reaches (Arcement and Schneider, 1989), and streambed hydraulic conductivity was initially set to 0.0154 m/d for all stream segments in the SFR Package. Streambed hydraulic conductivity for each stream segment was included as a parameter for calibration.

Streambed hydraulic conductivity in streams can differ based on differences in deposition and accumulation of sediment, and these natural variations in streambed hydraulic conductivity can affect seepage rates (Hatch and others, 2010). To better characterize variation in the Des Moines River streambed, apparent resistivity of the top 1 m of the streambed was used to generate a multiplier to linearly scale streambed hydraulic conductivity for individual reaches in areas where CRP data were collected (fig. 12). This approach assumes a linear correlation of unconsolidated sediment grain size to apparent electrical resistivity and hydraulic conductivity. The hydraulic conductivity multiplier for each SFR reach was calculated as the ratio of average streambed resistivity within a 60-m radius around the reach’s center point to average streambed resistivity for all reaches. This means that SFR reaches with lower-than-average apparent resistivity have a multiplier of less than 1.0, resulting in a hydraulic conductivity value lower than the uniform segment hydraulic conductivity value, and reaches with greater-than-average apparent resistivity have a multiplier greater than 1.0, resulting in a hydraulic conductivity value greater than the uniform segment hydraulic conductivity. Apparent resistivity values less than 50 ohm-meter (ohm-m) were identified as water column points, so the streambed was identified at each reach midpoint as the topmost data point with resistivity greater than 50 ohm-m, and apparent resistivity of the streambed was averaged from the top meter of the streambed. An additional scaling factor was used to define the slope of the linear relation between hydraulic conductivity and apparent resistivity of the streambed. The slope scaling factor was initially set to 1.0 and included as a parameter used in model calibration. If this relation resulted in a streambed hydraulic conductivity less than 1.0x10⁻³ m/d, streambed hydraulic conductivity was assigned to a minimum of 1.0x10⁻³ m/d.

Several water-filled quarries on the surface of the alluvial plain were represented using the MODFLOW 6 GHB Package, which simulates flow into or out of a cell from an external source as a proportion of the difference between the hydraulic head in the cell and the hydraulic head assigned to the external source (Langevin and others, 2017). GHB cells representing the quarries were assigned to layer 2 (river alluvium), and layer 1 was inactive in these areas. The water level for these quarry boundary cells was fixed at land surface elevation as determined by the DEM (Iowa Geospatial Data, 2022).

Groundwater Evapotranspiration (ET_g)

ET_g was simulated with the MODFLOW 6 EVT Package and applied to the uppermost cell in the model (at the land surface) in units of meters per day (Langevin and others, 2017). The EVT Package simulates the effects of plant transpiration and direct evaporation by removing water from the saturated groundwater regime (Langevin and others, 2017). Simulation of ET_g was based on the following assumptions (Langevin and others, 2017): (1) when the water table is at or above a specified elevation (the evapotranspiration surface), ET_g loss is at a specified fixed rate (the maximum ET_g rate); (2) when the depth of the water table below the evapotranspiration surface exceeds a specified value (called the “extinction depth” in this report), ET_g ceases; and (3) when between the evapotranspiration surface and the extinction depth, ET_g varies in a piecewise-linear fashion with water table elevation. Extinction depths were assigned on the basis of land cover, as defined using a simplified classification of the 2020 USDA Cropland Data Layer (USDA, 2020): the extinction depth of grassland, barren land, and developed areas was set to 1 m; farmland to 1.5 m; woodland to 3 m; and areas of open water were set to 0 m. Initial maximum monthly ET_g rates were assigned uniformly to the PRMS HRUs in the model area and were calculated using the PRMS model outputs that were sampled to each model cell from the nearest HRU polygon.

Groundwater Withdrawal

Groundwater withdrawal was represented with the MODFLOW 6 WEL Package (Langevin and others, 2017). RCWs were assigned to the cell closest to the center of the collector well, and other wells were assigned to the cell in which they were located. Monthly average groundwater withdrawal rates were specified for each well in units of cubic meters per day. The “auto flow reduction” option for the WEL Package reduces the groundwater withdrawal rate as the water level in a well approaches the bottom of the cell. The auto flow reduction option was set to 0.2 for all WEL Packages in the model, meaning the withdrawal rates could be reduced when the water level in a WEL Package cell was simulated at or less than 20 percent of the vertical cell thickness.

Calibration Targets

Observed values for calibration were determined for the steady-state and transient model stress periods. The first stress period of the model representing steady-state conditions was calibrated for hydrologic conditions before the 1980s and before substantial development of groundwater resources in the model area. The remaining monthly transient stress periods simulated monthly conditions for January 1, 1980, through July 31, 2022. Transient observations were consolidated into monthly values.

Observation data used for model calibration targets included groundwater levels, changes in groundwater levels, pumping rates for the RCWs, streamflow, and streambed flux (table 4). Groundwater level observations included water level elevations, water level changes, and water level inequality observations. Water level data were acquired from 13 USGS monitoring wells, including 4 with single measurements and 11 with multiple months of data (table 1; fig. 2). Additional water level data were from RCW records furnished by DMWW (Dylan Hart, Black & Veatch Corporation, written commun., 2023; Davis and Bristow, 2024), 10 monitoring wells in the alluvial plain used as part of a wetland habitat monitoring project between 2010 and 2015 (Vern Rash, DMWW, written commun., 2022 and Davis and Bristow, 2024), and 97 wells from the State of Iowa’s GeoSam database that had single water level measurements, generally collected at the time the well was drilled (Davis and Bristow, 2024). Water level elevation observations were assigned from wells with single water level measurements or from the first measurement for wells with a series of measurements. Monthly water level change observations were calculated for wells with multiple measurements. Daily pumping rates at the RCWs, averaged into monthly observations, were also used as calibration targets. Minimum design water levels for the collector wells, representing threshold groundwater levels below which pumping is reduced, were used as inequality observations. Inequality observations are a unique observation type implemented in PESTPP–IES (White and others, 2020) in which the residual is zero unless model output is greater than or less than the calibration target; the user specifies which inequality condition applies. In this case, the minimum water level for the collector well operations were used as a “greater than” observation, so simulated water levels greater than the minimum design water level were not penalized in calibration, but those lower than the minimum water level generated a residual.

Streambed flux (streamflow gain or loss) calibration targets were calculated as change in streamflow between the upstream and downstream USGS gages in the model area, and cumulative monthly streambed flux was calculated from the monthly streambed flux values. The Saylorville and Beaver Creek gages have monthly streamflow data available for the entire model period, January 1980 through July 2022, while the 2nd Avenue gage has streamflow data for October 1996 through July 2022; gage streamflow differences were calculated for the period from 1996 to 2022 with data for all three gages. Additionally, only measurements for October–March of each year were used during calibration to reduce expected effect of overland runoff on streamflow measured at gages.

Because this model is designed to simulate water levels in the alluvial aquifer in the area of the RCWs, groundwater observations from the alluvial aquifer materials were assigned the highest weights in calibration (observation groups wl_alt_obs(1), wl_alt_obs(3), and wl_alt_obs(4), table 4). Water level elevation observations for wells in layer 6 (wl_alt_obs(2); bedrock) were assigned a lower weight than those from the other layers because these are likely representative of bedrock aquifer units that underlie the shale bedrock confining layer at the bottom of this model and therefore are not well-connected to the Des Moines River alluvial aquifer. Streamflow targets were given relatively low weight in calibration because the model’s objective was to simulate groundwater levels rather than streamflow targets.

Calibration Parameters

Parameters adjusted during model calibration include hydrogeologic properties (horizontal hydraulic conductivity, vertical anisotropy, specific storage, and specific yield), streambed hydraulic conductivity, recharge multipliers, evapotranspiration multipliers, well withdrawal multipliers, quarry bed hydraulic conductivity, and groundwater flux boundary hydraulic conductivity (table 5). Horizontal hydraulic conductivity, vertical anisotropy, specific storage, and specific yield were distributed among the layers using a combination of hydrogeologic unit zones and pilot points (shown in figs. 13 through 18 in “Calibration and Numerical Model Results”). Distribution of pilot points and starting parameter values are described in the “Hydrogeologic Properties” section. Streambed hydraulic conductivity for the SFR Package was assigned uniformly between stream confluences for all streams other than the Des Moines River, which used variable streambed hydraulic conductivity based on the variation of cell apparent electrical resistivity using a multiplier. An additional calibrated parameter was used to adjust the slope of the resistivity multiplier relation. Recharge and evapotranspiration rate multipliers were assigned by HRU; recharge rate multipliers for soil type zones were also used. Well pumping rate multipliers were used for individual wells. Calibration parameters also included quarry-bed hydraulic conductivity for each quarry, and GHB hydraulic conductivity for groundwater flux boundaries was assigned uniquely for each boundary. All calibration parameters, weights, and target observations are included in the associated groundwater model data release (Davis and Bristow, 2024).

Optimal Parameter Estimates

Optimal model parameters were selected from the “base” realization of the second PESTPP-IES iteration. In general, calibration of the model resulted in parameter values similar to initial estimates. Calibrated parameter values are summarized in tables 6 and 7, and complete calibration results are included in the data release accompanying this report (Davis and Bristow, 2024). Calibrated horizontal conductivity values for each model layer (table 6) are illustrated along with pilot points and hydrologic parameter zones in figures 13 through 18. In layers 2 and 5, areas of the highest hydraulic conductivity were near the river channel and RCWs.

Calibrated recharge multipliers for the three soil type zones in layer 1—(1) manmade fill, (2) till and bedrock-based surficial material, and (3) alluvial, eolian, and outwash-based surficial material—were 0.05, 0.208, and 2.09, respectively. These zones are the same as the hydraulic conductivity parameterization zones in layer 1, delineated in figure 13. The calibrated recharge multiplier for till- and bedrock-based surficial material (0.208) is much smaller than the initial estimate of 0.5, while the calibrated value for alluvial, eolian, and outwash-based surficial material (2.09) is higher than the initial estimate of 1 (table 7).

Comparison of Observed and Simulated Water levels

Model performance was assessed on its ability to match hydraulic observations of interest, including water level elevations, monthly average changes in water levels (water level changes), minimum design water levels for the RCWs, and monthly average streambed flux. Plots comparing observed and simulated water levels at the RCWs (fig. 19) and other observation wells (fig. 20) demonstrate the ability of the model to simulate observed water levels at the RCWs and to satisfy the RCW minimum design water level (figs. 19–20). In addition, observed and simulated water table elevations were plotted and compared to a 1:1 perfect-fit line for 119 wells with single water level measurements used as calibration data (figs. 21A, B; table 4). The model accurately simulated water levels and seasonal water level fluctuations at the RCWs, and it accurately simulated water table levels above the RCW design minimum level for all stress periods with RCW water level observation data (figs. 19A, B). Average absolute value water level residuals for the model period with available observations are 1.20 m at RCW1 and 1.49 m at RCW2 (Davis and Bristow, 2024). The water level scatter plot demonstrates a slight positive bias in simulated water levels (figs. 21A, B). Comparison of model output water levels to observed values by layer is instructive. The average absolute value of residuals in layers 2 and 5 (the model layers representing the alluvial aquifer) is 1.57 m and 2.47 m, respectively, while that in layer 6 (bedrock) is 24.22 m (table 8). Layer 6 represents a homogenous shale bedrock confining unit for this model, which is appropriate for the model’s focus on the alluvial aquifer. Water level observations assigned to layer 6 are likely representative of water levels in deeper bedrock aquifers that are not connected to the alluvial aquifer targeted in this model and are thus poorly simulated by this model.

Simulated Groundwater Budgets

Simulated groundwater budgets were calculated using Zone Budget version 6, a program packaged with MODFLOW 6 that characterizes model inflows and outflows specified in the model output files. Water budgets with model inflows and outflows were calculated for each monthly stress period of the transient model (table 9; fig. 22). The average water budget, in cubic meters per day, for the entire transient model period (1980 through 2022) is provided in table 9, along with the average water budget before and after the introduction of RCWs. Water budgets were also averaged into annual total flux volumes (fig. 22). All zone budget output files are available in the data release accompanying this report (Davis and Bristow, 2024).

Simulated inflows to the transient model include recharge from precipitation, infiltration from streams, infiltration from quarries, groundwater inflow through flux boundaries at the model edge, and gain from storage. Outflows include evapotranspiration, discharge to streams, discharge to quarries, outflow through flux boundaries, groundwater withdrawal, and loss to storage. The largest average model inflows for the entire transient period were recharge from precipitation (35,210 m³/d or 5.8 centimeters per year [cm/y] across the active model area) and infiltration from streams (20,573 m³/d; table 9). The largest average model outflows for the entire transient period were evapotranspiration (38,242 m³/d or 6.3 cm/y) and discharge to streams (12,176 m³/d). RCWs became operational in late 2011 and after that date, withdrawals average 21,177 m³/d (table 9). The largest magnitude of change from the period before RCW initial operation to after is in storage. Before October 2011, average net gain from storage (gain into the aquifer from storage minus loss from the aquifer to storage) was negative in every layer (−524 m³/d), and after that date average net gain is positive in every layer (a combined average of 3,437 m³/d). Other notable changes in groundwater fluxes after the initiation of RCW operation include an increase in infiltration from quarries and decrease in discharge to quarries and increase in infiltration through streambeds. The water budget also indicates substantial seasonal and annual fluctuation in net storage and other budget components (fig. 22) before and after the initiation of RCW groundwater withdrawals. Loss from the aquifer to storage is highest in the late spring and early summer months, and withdrawals from storage highest in later summer and autumn months.

Comparison of observed streambed flux with simulated net discharge through streambeds (discharge minus infiltration) indicates a substantial mismatch in the magnitude of the values (fig. 23; note the different Y-axis scales for each dataset). Because this model was designed to simulate water levels at the RCWs in the Des Moines River alluvial aquifer and not streamflow, streamflow calibration targets were weighted much lower than groundwater level targets in calibration (table 4). Observation data for seepage to streams were calculated as the difference between upstream and downstream streamflow measurements, ignoring the contribution of overland surface runoff to streamflow. Observation data were calculated for late fall and winter months only (October through March) in an effort to minimize the contribution of runoff to this net streamflow calculation, but this may have been insufficient to control for overland runoff.

Sensitivity of Model Parameters

Model parameter sensitivity is a measure of the effect that a change in model input, such as changes to parameter values, has on outputs, such as calibration targets. During calibration, PESTPP–IES runs the model repeatedly while slightly adjusting parameter values for each run and then compares model output results to output from the previous model run. Parameter sensitivity is a measure of how much output values change in proportion to how much a parameter has changed; high sensitivity parameters have a relatively higher effect on model results than low sensitivity parameters. Parameter sensitivity tests for this model indicated that horizontal hydraulic conductivity for pilot points in three zones of layer 5 (zones 10, 11, and 13) were the most sensitive calibration parameters (fig. 24). These zones represent buried valley alluvium, buried valley and alluvium and outwash combined, and buried valley fill in the Rock Creek valley. The RCWs targeted by this study withdraw water from zone 11, which represents the bottom part of the alluvial aquifer around where the wells are located. SFR streambed hydraulic conductivity also had relatively high sensitivity. Parameter sensitivity was analyzed by parameter group averages; for parameter groups with more than one individual value, such as hydraulic conductivities distributed to pilot points or recharge multipliers for each HRU, the sensitivity value presented here is the average of all group members. All parameter sensitivity results are included in the accompanying model archive data release (Davis and Bristow, 2024).