Hydrogeologic Units

The hydrogeologic framework consists of the five hydrogeologic units described in the conceptual model: (1) fine-grained alluvium, (2) coarse-grained alluvium, (3) glacial drift, (4) buried-valley sediment, and (5) bedrock. Each of these hydrologic units was incorporated as a distinct layer of the hydrogeologic framework and numerical model. These units were identified, delineated, and interpolated in three dimensions for the framework using previous lithostratigraphic models and existing maps of bedrock altitude and surficial geology (Horick, 1973; Hallberg, 1980; Runkle, 1985; Miller and Burras, 2015; Iowa Geological Survey, 2020). Hydrologic parameter zones were assigned to each layer; the glacial drift and buried-valley sediment were assigned one zone each, and the fine-grained alluvium, coarse-grained alluvium, and bedrock layers were subdivided into multiple hydrogeologic zones based on compositional or landform variations within the layer.

Fine- and Coarse-Grained Alluvium

The Iowa River alluvial aquifer consists of fine- and coarse-grained alluvial units in the model area. Fine- and coarse-grained alluvium is present in the historical and present channels of the Iowa River and its tributaries (“river valleys”) and consists of sand and gravels that were deposited by these rivers and streams. The thicknesses of the alluvial units in the model area were estimated using land surface altitude and driller’s logs, where available, in the model area (Iowa Geological Survey, 2024). Land surface altitude for the hydrogeologic framework was determined using the mean land surface altitude over the 100-m model cell size from a hydro-conditioned and hydro-enforced digital elevation model (DEM; Iowa State University, 2020).

The extent of the alluvial units in the model area was estimated using a map of alluvial deposits in Iowa from the Iowa Department of Natural Resources (DNR; Iowa Department of Natural Resources, 2020). Thickness of the alluvial units was estimated using driller’s logs from the Iowa GeoSam database (Iowa Geological Survey, 2024), and logs from wells drilled near the Meskwaki Settlement were prioritized. Fine-grained alluvium, where present, was assigned a thickness in the active model area ranging from about 5 to 19.5 m, and coarse-grained alluvium, where present, was assigned a uniform thickness of 13.7 m (figs. 3 and 4). The fine- and coarse-grained alluvium were each divided into three zones based on their location: fine-grained alluvium in the Iowa River valley (zone 1); fine-grained alluvium in the valleys of the major tributaries of the Iowa River (zone 2); fine-grained alluvium underlying minor tributaries of the Iowa River which, in the model area, are generally third- or fourth-order streams with respect to the Iowa River (zone 3); coarse-grained alluvium in the Iowa River valley (zone 4); coarse-grained alluvium in the valleys of the major tributaries of the Iowa River (zone 5); and coarse-grained alluvium underlying minor tributaries of the Iowa River (zone 6; table 2).

Table 2.    

Hydrologic properties of each hydrogeologic unit in the numerical model.

[K_h, horizontal hydraulic conductivity; m/d, meter per day; K_z/K_h, vertical hydraulic conductivity; S_s, specific storage; m^-1, per meter; S_y, specific yield]

Table 2.    Hydrologic properties of each hydrogeologic unit in the numerical model.

Layer	Zone number	Description	Kh (m/d)	Kz/Kh ratio (percent)	Ss (m-1)	Sy (percent)	Number of pilot points
1	1	Fine-grained alluvium in Iowa River valley	20	10	1.00x10−5	5	92
	2	Fine-grained alluvium in valleys of major tributaries of the Iowa River	20	10			119
	3	Fine-grained alluvium underlying minor tributaries of the Iowa River	10	10			296
2	4	Coarse-grained alluvium in Iowa River valley	100	10	3.54x10−5	15	92
	5	Coarse-grained alluvium in valleys of major tributaries of the Iowa River	100	10			109
	6	Coarse-grained alluvium underlying minor tributaries of the Iowa River	50	10			296
3	7	Glacial drift	0.001	10	2.74x10−7	3	258
4	8	Buried valley sediment	50	10	5.25x10−5	10	99
5	9	Pennsylvanian shale	0.0003	10	2.74x10−7	1	11
	10	Mississippian limestone/dolomite, higher Kh	142	10			62
	11	Devonian shale	0.0003	10			136
	12	Silurian-Devonian limestone/dolomite	2.2	10			37
	13	Mississippian limestone/dolomite, lower Kh	0.5	10	5.25x10−6	5	23

Glacial Drift

In the valleys of the Iowa River and its tributaries, glacial drift has been eroded by streamflow, and alluvial sediments are present at the land surface. Outside of the river valleys, glacial drift typically overlies the bedrock unit and is exposed at the land surface. The thickness of the glacial drift unit not within the river valleys was estimated by taking the difference between land surface altitude, estimated from a DEM of the model area (Iowa State University, 2020), and the altitude of the top of the bedrock unit or, where present, the altitude of the top of the buried-valley sediment (fig. 5). In the river valleys, the thickness of the glacial drift unit was estimated by taking the difference between the bottom of the coarse-grained alluvial unit and the top of the bedrock unit or, where present, the altitude of the top of the buried-valley sediment. Where glacial drift is present within the model area, its thickness ranges from about 5 m in the river valleys to about 130 m outside the river valleys (fig. 5). Glacial drift in the model area was assigned to zone 7 (table 2).

Buried-Valley Sediment

Buried valleys are ancient river channels that carved into underlying bedrock and were subsequently filled in with alluvial and glacial outwash sand and gravels (buried-valley sediment). Buried-valley aquifers are documented in the model area, including the Poweshiek and Belle Plaine channels (Prior and others, 2003). The exact location and extent of buried-valley sediment can be difficult to determine since there is little to no evidence of its presence at the land surface. The location and thickness of buried-valley sediment within the model area was therefore estimated by geo-referencing mapped buried-valley aquifers in Iowa provided by the Iowa DNR (Prior and others, 2003) and outlining the major buried valleys (fig. 6). These surfaces were then intersected with bedrock altitude and interpolated using the nearest-neighbor method. This method was confirmed by the Iowa DNR as the same method used to produce their map of buried-valley aquifers (Stephanie Tassier-Surine, Iowa Department of Natural Resources, written commun., 2023). Where buried-valley sediment is present in the model area, its thickness ranges from about 5 to 75 m (fig. 6). Buried-valley sediment in the model area was assigned to zone 8 (table 2).

Bedrock

The altitude of the top of the bedrock unit in the model was estimated using Iowa Geological Survey bedrock altitude maps (rasters) and lithologic logs from the GeoSam database (Iowa Department of Natural Resources, 2017; Iowa Geological Survey, 2024). The top altitude of the bedrock unit was determined using the mean bedrock surface altitude over the 100-m cell size (fig. 7). The bottom of the bedrock unit was uniformly set to 200 m below the top of the bedrock unit. The bedrock hydrogeologic units slope toward the southwest and have been eroded across their upper leading edges (Prior and others, 2003), which results in younger layers being exposed at the bedrock surface in the west in the model area, with successively older bedrock units expressed at the bedrock surface toward the east (fig. 2). There are four upper bedrock units present within the study area, which were divided into five zones: (1) Pennsylvanian shale, which forms the Pennsylvanian confining unit, was assigned zone 9; (2) Mississippian limestone/dolomite, which forms the Mississippian aquifer, was assigned zones 10 and 13; (3) Devonian shale, which forms the Devonian confining unit, was assigned zone 11; and (4) Silurian-Devonian limestone/dolomite, which forms the Silurian-Devonian aquifer, was assigned zone 12. The Mississippian limestone/dolomite was divided into two zones, an area of higher horizontal hydraulic conductivity (zone 10) and an area of lower horizontal hydraulic conductivity (zone 13), based on results from initial numerical model calibration. The Pennsylvanian and Devonian shale act as confining units in the model area, meaning they do not readily transmit groundwater, while the Mississippian and Silurian-Devonian limestone/dolomite are considered aquifers in the model area as they contain sufficient saturated permeable material to yield substantial quantities of water to wells.

Precipitation Recharge and Evapotranspiration

In the model area, the Iowa River alluvial aquifer, glacial drift, and bedrock outcrops receive recharge from infiltration of precipitation that falls on the land surface. ET_g in the model area is evaporation from the saturated groundwater zone as well as from plant transpiration (Wilson and Moore, 1998). Recharge from precipitation and ET_g were estimated using a Precipitation-Runoff Modeling System (PRMS) model that was previously created for nine river basins in eastern Iowa, including the Iowa River Basin (Haj and others, 2015). The PRMS model of the Iowa River Basin was calibrated for the period of October 2002 through September 2012 to estimate daily streamflow at sites without streamgages within the basin. PRMS is a deterministic, distributed-parameter, physical-process-based modeling system developed to evaluate the response of various combinations of climate and land use on streamflow and general watershed hydrology (Markstrom and others, 2015). It simulates hydrologic processes including evaporation, transpiration, runoff, infiltration, and interflow using computer source code called “modules” and can simulate hydrologic water budgets at the watershed scale for temporal scales ranging from days to centuries (Markstrom and others, 2008).

In PRMS, basins are subdivided into a series of contiguous and internally homogenous spatial units called hydrologic response units (HRUs) based on hydrologic and physical characteristics such as land surface altitude, slope, aspect, plant type and cover, land use, soil morphology, geology, drainage boundaries, distribution of precipitation, temperature, solar radiation, and flow direction (Markstrom and others, 2008). The area of an HRU is a type of PRMS “parameter,” or an input value that does not change during the course of a simulation. The parameters required for a simulation depend on which modules are used. HRUs produce streamflow and transfer it to other HRUs and to the drainage network, represented by stream segments (Markstrom and others, 2015). The Iowa River Basin PRMS model consists of 2,340 HRUs and 1,174 stream segments; of these, 165 HRUs and 88 stream segments are within the active model area of this study.

The PRMS model of the Iowa River Basin was run for the HRUs and stream segments located within the active model area using a temporally expanded input dataset for the period of January 1, 1980, through August 30, 2022, to generate input datasets for the MODFLOW 6 model (figs. 8 and 9). In PRMS, a user specifies which variables (values that can change each time step during a simulation) to compute and output. The implemented variables for the PRMS HRUs were “recharge,” which is the simulated recharge to the groundwater system, and “potet,” which is the simulated potential evapotranspiration (PET) for each HRU. Each model cell was spatially joined to its associated HRU from the PRMS model using ArcGIS, and the recharge and PET values calculated for an HRU were assigned to each active model cell.

Groundwater Withdrawal

The Meskwaki Settlement has several drinking water production wells completed in the unconsolidated Iowa River alluvial aquifer. Monthly groundwater withdrawal data for these wells and screened well depths were provided by the Meskwaki Natural Resources Department (Leland Searles, Meskwaki Natural Resources Department, written commun., 2023). Supply wells were first installed in 1984 and currently serve a population of approximately 1,100 people on the Meskwaki Settlement. Pumping records were available on a monthly basis from 2018 to 2022 as a total volume withdrawn from all wells.

Groundwater withdrawals from an additional 39 production wells within the model area were included in the model. These volumes were obtained from the Iowa DNR Water Allocation Compliance and Online Permitting database (Iowa Department of Natural Resources, 2024). Monthly records of pumping were available for 2011–22 for all wells, and annual records were available for 1993–2022 for most wells. Groundwater withdrawals for irrigation, domestic, or livestock uses are relatively small compared to withdrawals for municipal and industrial uses in the model area and were not included in the numerical model.

Groundwater and Surface-Water Interactions

Water in the Iowa River alluvial aquifer is assumed to be in hydraulic connection with the Iowa River and its tributaries in the model area (Prior and others, 2003); therefore, groundwater discharge to, and groundwater recharge from, surface water likely occurs in the model area. Base flow is the portion of streamflow that is not attributable to direct runoff or melting snow, and it is sustained by groundwater discharge (Barlow and Leake, 2012). During dry periods, the base flow of the Iowa River is primarily water released from groundwater storage, and during periods of high streamflow, water from the river recharges the aquifer through the streambanks (Prior and others, 2003). Base flow can be calculated from streamflow recorded at streamgages using a process called hydrograph separation, in which a hydrograph is graphically separated into base flow and the portion of the streamflow hydrograph that exceeds base flow (Koch, 1970). The part of the streamflow hydrograph that exceeds base flow is referred to as runoff or residual flow (Koch, 1970). Base-flow gains and losses can then be estimated by subtracting base flow calculated using hydrograph separation at an upstream streamgage (or, if there are multiple upstream streamgages, the sum of base flow from all contributing upstream streamgages) from base flow calculated at a downstream streamgage (Koch, 1970). The calculated difference in base flow between a downstream streamgage and any upstream streamgage(s) that contributed to the downstream streamflow in the model area for a given month is the groundwater flux through the streambed between two or more streamgages and is hereafter referred to as “streambed flux.”

Streamflow data were available for USGS streamgages in the model area (U.S. Geological Survey, 2024; fig. 1; table 1). Base flow for each streamgage was calculated with the HYSEP sliding interval method of hydrograph separation (Sloto and Crouse, 1996), which uses an algorithm to draw connecting lines between low points of the streamflow hydrograph, using the software Groundwater Toolbox (Barlow and others, 2015). For the common period of record from 2015 to 2022, streambed flux was calculated in the model area by subtracting the base flow at the upstream USGS streamgages (Iowa River at Marshalltown, Iowa; Timber Creek near Marshalltown, Iowa; Iowa River at County Highway E49 near Tama, Iowa; Richland Creek near Haven, Iowa; Walnut Creek near Hartwick, Iowa; and Salt Creek near Elberon, Iowa) from the streamflow at downstream USGS streamgage Iowa River near Belle Plaine, Iowa, for the same period. For 1980–2015, the model period not within the common period of record for these streamgages, a percent difference ratio of median winter (October–March) base flows throughout the entire model period (1980–2022) to median winter base flows for the common period of record (2015–22) was calculated for the USGS streamgage Iowa River at Marshalltown, Iowa, which was the streamgage with the longest period of record in the model area. The calculated percent difference ratio of median winter base flows at USGS streamgage Iowa River at Marshalltown, Iowa, was then multiplied by the base flow recorded during the common period of record at the other streamgages in the model area to estimate base flow during the model period for which there was no recorded base flow. The streambed flux values tended to be low during lower streamflow, and higher streambed flux values typically corresponded to higher streamflow in the Iowa River.

The PRMS model of the Iowa River Basin was also used to calculate inflow to the model from the Iowa River and Linn Creek, a tributary of the Iowa River, where it enters the active model area. The PRMS variables implemented for the stream segments were “seg_inflow,” which is the total streamflow entering a segment, and “seg_outflow,” which is the total streamflow leaving a segment.

Spatial and Vertical Discretization

The active model area was designed to include the total area of the Meskwaki Settlement and is based on the HUC–10 boundary of the Iowa River and its major tributaries surrounding the Meskwaki Settlement. This area was chosen to include the USGS streamgage Iowa River at Marshalltown, Iowa, located upstream of the Meskwaki Settlement, and USGS streamgage Iowa River near Belle Plaine, Iowa, located downstream from the Meskwaki Settlement, for the purposes of streamflow routing and assigning calibration targets. The numerical model was spatially discretized into 100-m square cells of varying thickness and arranged in a grid with five layers, 496 rows, and 856 columns (fig. 10).

The model was vertically discretized into five layers based on the extent and thicknesses of the five hydrogeologic units described in the hydrogeologic framework: layer one consists of fine-grained alluvium, layer two consists of coarse-grained alluvium, layer three consists of glacial drift, layer four consists of buried-valley sediment, and layer five consists of bedrock. Because MODFLOW 6 allows units to “pinch out” where they are not present due to lateral tapering or thinning out, each layer corresponds to one hydrogeologic unit. If an inactive, “pinched out” cell is vertically between two active cells, the inactive cell was designated as a “vertical pass-through cell” using the IDOMAIN array in the DIS package (Langevin and others, 2017). MODFLOW 6 directly connects the active cells overlying and underlying the vertical pass-through cells, effectively removing them from the simulation (Langevin and others, 2017). The active model area and layering includes 808,265 active cells, excluding pass-through cells.

The altitude of the top of the model was estimated using a hydro-enforced and hydro-conditioned DEM (Iowa State University, 2020), as described in the “Hydrogeologic Units” section of this report. The bottom altitude of the model was set at an arbitrary 200 m below the top of the bedrock unit to give the bottom layer of the model a uniform thickness. Altitudes of the top and bottom of each cell for each layer used in the numerical model were estimated using the methods described in the “Hydrogeologic Units” section of this report and are available in the accompanying USGS data release (Davis and Goldstein, 2025).

Temporal Discretization

The numerical model for the Iowa River alluvial aquifer and underlying hydrogeologic units was temporally discretized using blocks of time called “stress periods,” within which all hydrologic stresses are constant and numerical model output applies to the end of each stress period. The groundwater flow model was temporally discretized into a single steady-state stress period followed by 512 monthly transient stress periods. The steady-state part of the model was constructed to approximate steady-state (mean) conditions and was assumed to generally represent low-flow conditions. Inputs to the steady-state part of the model were assumed to be representative of steady-state conditions and include recharge and ET_g, which was based on the mean daily rate from the PRMS model over the period 1980–2021, and mean streamflow data from 1979. Aquifer storage is considered negligible under steady-state conditions and is therefore not included in formulations of the groundwater flow equation for steady-state models (Langevin and others, 2017). The steady-state part of the model was used to provide starting conditions for the time-varying transient model.

The transient part of the model was used to simulate conditions for January 1, 1980–August 30, 2022, using monthly stress periods and the temporal unit of days. Aquifer stresses during the transient stress periods were representative of the mean of hydrologic stresses during each month. The rate of change in aquifer storage in a model cell during a transient stress period is equal to the sum of all flows into and out of the model cell and is not required to equal zero like in the steady-state part of the model. Aquifer storage was therefore considered during the transient formulations of the groundwater-flow equation calculated by MODFLOW 6.

Hydrologic Properties

The model utilized a series of discrete locations (pilot points) to which initial parameter values were applied and allowed to change during model calibration. Pilot points were evenly distributed within each hydrogeologic unit zone. The final distribution of hydrologic properties in each zone was interpolated to the model grid from pilot point values determined during model calibration using the kriging tool in ArcGIS (Esri, 2024).

Initial values of horizontal hydraulic conductivity were assigned at pilot points for the fine-grained (layer 1) and coarse-grained (layer 2) alluvial units within the model area based on the zones from the hydrogeologic framework: zone 1 was for fine-grained alluvium in the valleys of the Iowa River, zone 2 was for fine-grained alluvium in the valleys of the major tributaries of the Iowa River, zone 3 was for fine-grained alluvium underlying minor tributaries of the Iowa River, zone 4 was for coarse-grained alluvium in the valleys of the Iowa River, zone 5 was for coarse-grained alluvium in the valleys of the major tributaries of the Iowa River, and zone 6 was for coarse-grained alluvium underlying minor tributaries of the Iowa River (fig. 11A, 11B; table 2). Initial horizontal hydraulic conductivity values were based on values for similar hydrogeologic materials suggested in Morris and Johnson (1967), Freeze and Cherry (1979), and Bayless and others (2017). For the fine-grained alluvium in layer 1, pilot points in zones 1 and 2 were assigned hydraulic conductivities of 20 meters per day (m/d) and in zone 3 they were assigned a value of 10 m/d (table 2). For the coarse-grained alluvium in layer 2, pilot points in zones 4 and 5 were assigned a value of 100 m/d and in zone 6 they were assigned a value of 50 m/d (table 2). Horizontal hydraulic conductivity pilot points for the glacial drift (layer 3) were initially set to a uniform value of 0.001 m/d using one zone (zone 7; table 2), and horizontal hydraulic conductivity pilot points for the buried-valley alluvium (layer 4) were initially set to a uniform value of 50 m/d using one zone (zone 8; table 2). Bedrock in layer 5 was subdivided into five zones: zone 9 for Pennsylvanian shale, zone 10 for the higher conductivity Mississippian limestone/dolomite, zone 11 for Devonian shale, zone 12 for Silurian-Devonian limestone/dolomite, and zone 13 for the lower conductivity Mississippian limestone/dolomite (fig. 11C; table 2). The pilot points in these zones were initially assigned horizontal hydraulic conductivity values of 0.0003, 142, 0.0003, 2.2, and 0.5 m/d, respectively (table 2). Vertical hydraulic conductivity values for all layers were initially assigned using a vertical hydraulic conductivity anisotropy ratio that was 10 percent of the horizontal hydraulic conductivity at the same location (Todd, 1980).

Aquifer storage properties, used during the formulations of the groundwater-flow equation during the transient part of the simulation, consisted of specific yield and specific storage. These properties were assigned as constant for each model layer using pilot points and subsequently adjusted further during numerical model calibration. Initial values for specific yield and specific storage for all layers were from published ranges (Freeze and Cherry, 1979; Bayless and others, 2017) for similar hydrogeologic materials. Fine-grained alluvium (layer 1, zones 1–3) had an initial specific storage of 1.00x10⁻⁵ inverse meters (m⁻¹) and an initial specific yield of 5 percent; coarse-grained alluvium (layer 2, zones 4–6) had an initial specific storage of 3.54x10⁻⁵ m⁻¹ and an initial specific yield of 15 percent; glacial drift (layer 3, zone 7) had an initial specific storage of 2.74x10⁻⁷ m⁻¹ and an initial specific yield of 3 percent; buried valley sediment (layer 4, zone 8) had an initial specific storage of 5.25x10⁻⁵ m⁻¹ and an initial specific yield of 10 percent; bedrock other than the lower K_h Mississippian limestone/dolomite (layer 5, zones 9–12) had an initial specific storage of 2.74x10⁻⁷ m⁻¹ and an initial specific yield of 1 percent; and the lower K_h Mississippian limestone/dolomite (layer 5, zone 13) had an initial specific storage of 5.25x10^-6 m⁻¹ and an initial specific yield of 5 percent (table 2).

Boundary Conditions

Aquifer stresses in the numerical model were simulated using various MODFLOW 6 Stress and Advanced Stress packages. Lateral and lower numerical model boundaries were simulated as no-flow boundaries and, as such, no specific package was used. Recharge from infiltration of precipitation in the model area was simulated using the RCH package; direct ET_g was simulated using the EVT package; groundwater withdrawal was simulated using the WEL package; and the interaction of groundwater and surface water in the model area was simulated using the SFR package.

Surface-Water Features

Surface-water features in the model area contribute water to or drain water from the groundwater system depending on the hydraulic gradient between the surface-water feature and the groundwater head. Streams and rivers were the only surface-water features represented in the numerical model and were simulated using the SFR package (Langevin and others, 2017). The SFR package uses the continuity equation and assumption of piecewise steady (nonchanging in discrete periods), uniform (nonchanging in location), and constant-density streamflow (Langevin and others, 2017). This assumption ensures that during all times, volumetric inflow and outflow rates are 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 the report; Langevin and others, 2017). In addition to calculating streamflow, the SFR package also calculates seepage between the surface-water feature and the aquifer. The stream network represented by the SFR package was divided into reaches with each reach representing a section of a stream that corresponds to a specific model cell (Langevin and others, 2017). Reaches were numbered sequentially beginning at USGS streamgage Iowa River at Marshalltown, Iowa, and ending at USGS streamgage Iowa River near Belle Plaine, Iowa (fig. 12). The sequential numbering provided a reference for routing surface-water flow among the SFR reaches.

Attributes required for each SFR reach were length, width, gradient, altitude of the top of the streambed, streambed hydraulic conductivity, streambed thickness, and Manning’s roughness coefficient of the stream channel (Arcement and Schneider, 1989). Stream length and width and streambed top altitude and thickness were assigned in units of meters, river gradients were assigned in units of meters per meter, streambed hydraulic conductivity values were assigned in units of meters per day, and Manning’s roughness coefficients were unitless. Additional information used in the SFR package included surface-water inflow to SFR reaches at USGS streamgages Iowa River at Marshalltown, Iowa, and Iowa River near Belle Plaine, Iowa. SFR reaches were assigned to the highest active numerical model layer. The length of each SFR reach was calculated using ArcGIS geoprocessing tools (Esri, 2024) by intersecting the stream features with the model grid, and the width of each SFR reach was estimated from aerial imagery. The stream network in the model area was segmented at river and tributary confluences (hereafter referred to as “segments”; fig. 12). The length of each segment was calculated using ArcGIS geoprocessing tools (Esri, 2024), and the upstream and downstream altitudes of each segment were determined based on a DEM for the model area (Iowa State University, 2020). Downstream altitudes that were either equal to or higher than upstream altitudes, as determined from the DEM, were manually adjusted and set equal to 0.1 m below the upstream altitude. The upstream and downstream altitudes and length of each segment were used to calculate river gradients for each SFR reach. River gradients less than 1.0x10⁻⁴ meters per meter were assigned a minimum gradient of 1.0x10⁻⁴ meters per meter. The upstream and downstream altitudes, length of each segment, and length of each SFR reach were used to calculate the riverbed top for the midpoint of each SFR reach. Riverbed thicknesses for each SFR reach were assigned a nominal and uniform value of 1.0 m. If the riverbed bottom altitude for an individual SFR reach (calculated by subtracting riverbed thickness from top altitude) was below the bottom of the uppermost active layer, the layer bottom altitude for that cell was assigned to 1.0 m below riverbed bottom altitude. In some cases, this method produced reach altitudes that were above the top of the model layer to which the reach was assigned. The reach altitudes that were above the top of the numerical model were manually adjusted to ensure the altitude was below the top of the numerical model layer to which each reach was assigned. SFR reaches were initially assigned a uniform riverbed hydraulic conductivity of 0.3048 m/d and a uniform Manning’s roughness coefficient of 0.037; however, these values were adjusted during numerical model calibration, as described in the “Numerical Model Calibration Approach” section of this report.

Calibration Targets and Weighting

Model calibration is completed using a process called history matching. History matching refers to the process of attempting to closely match real-world hydrologic observations of interest (calibration targets) to model-simulated values for a defined period (Welter and others, 2015). In general, calibration targets are direct measurement values, such as hydraulic head. However, in some cases, calibration targets may be calculated from measurement values, such as changes in hydraulic head. History matching is accomplished by attempting to reduce the sum of squared weighted residuals, or the objective function, for all calibration targets (Welter and others, 2015). Each calibration target in the history matching process is assigned a weight, and a residual is the difference between a calibration target and model-simulated value. Each observation is assigned to a group of similar observation types. Calibration targets for hydraulic head and changes in hydraulic head were assigned to groups based on model layer, and calibration targets for monthly mean streambed flux and cumulative monthly streambed flux each were assigned to their own group (table 4). Weighting adjusts the contribution of each residual to the overall objective function and combining calibration targets into groups ensures that no single observation type or group becomes a dominant contributor in the calibration process (Welter and others, 2015). Highly weighted calibration targets will have greater importance in the calibration process compared to calibration targets with lower weights because smaller changes in a residual for highly weighted calibration targets will generate larger changes to the overall objective function (Welter and others, 2015).

Calibration Parameters

Model parameters adjusted during numerical model calibration include horizontal hydraulic conductivity, vertical anisotropy ratio, specific storage, specific yield, streambed hydraulic conductivity, recharge, evapotranspiration from groundwater, and well multipliers. Horizontal hydraulic conductivity, vertical anisotropy ratio, specific storage, and specific yield were estimated at pilot point locations distributed within the model area; recharge and ET_g multipliers were estimated for each of the PRMS HRUs; well pumping multipliers were estimated for each well; and streambed hydraulic conductivity values were estimated for each stream segment in the SFR package (fig. 12). The upper and lower bounds for each calibration parameter, except for the recharge, ET_g, and well multipliers, were set at one-tenth and ten times the initial value. The estimates of recharge, ET_g, and pumping from wells were adjusted independently during the calibration process using multipliers. The multipliers were applied uniformly to all numerical model stress periods (steady-state and transient stress periods) for the spatially varied PRMS estimates used in the RCH and EVT packages and for each production well represented in the model. The initial value of the recharge and ET_g multipliers was set at 1 and allowed to range from 0.33 to 1.66, and the initial value of the well multiplier was set at 1 and allowed to range from 0.85 to 1.20.

Horizontal and vertical hydraulic conductivity, specific storage, and specific yield were calibrated at pilot point locations (figs. 15 and 16). Pilot points are discrete locations distributed throughout the active model area that represent surrogate parameters from which hydrologic property values are interpolated to the numerical model grid using a kriging algorithm (Doherty and Hunt, 2010). The use of pilot points and the method of kriging allows for greater variation in hydrologic properties compared to using either a single value for each model layer or subdividing each layer into piecewise zones of uniform hydraulic conductivity (Doherty and Hunt, 2010). As a result, the process allows a smooth variation in hydrologic properties in the model area interpolated from estimated values at pilot points. Final parameter estimates were checked to ensure that they were appropriate based on literature values.

Optimal Parameter Estimates

The final calibrated parameter values of horizontal and vertical hydraulic conductivity, specific yield, specific storage, streambed hydraulic conductivity, recharge, and ET_g were considered reasonable for the hydrogeologic materials and conditions in the model area for January 1980–August 2022. In general, the mean calibrated parameter values were reasonably close to previously published values (table 5). The mean horizontal hydraulic conductivity for the fine-grained alluvium, estimated at pilot points during the calibration process, was about 36 m/d in the Iowa River valley (zone 1), 33 m/d in valleys of the major tributaries of the Iowa River (zone 2), and 17 m/d underlying minor tributaries of the Iowa River (zone 3; fig. 14A; table 5). Calibrated horizontal hydraulic conductivity ranged from about 2 to 179 m/d, 2 to 200 m/d, and 1 to 100 m/d at pilot point locations in zones 1, 2, and 3, respectively (fig. 14A; table 5). The mean calibrated vertical anisotropy ratio for the fine-grained alluvium was 25 percent at pilot points in zone 1, 15 percent in zone 2, and 17 percent in zone 3 (fig. 15A), which were generally similar to the initial value of 10 percent (table 5), and ranged from about 1 to 100 percent at pilot point locations in zones 1, 2, and 3 (table 5; fig. 15A). The mean calibrated specific storage was 1.58x10⁻⁵ m⁻¹ in zone 1, 2.08x10⁻⁵ m⁻¹ in zone 2, and 1.75x10⁻⁵ m⁻¹ in zone 3, which are similar to the initial specific storage for all three zones of 1.00x10⁻⁵ m⁻¹ (table 5). For specific yield of the fine-grained alluvium, the mean calibrated values were 8 percent in zones 1 and 2, and 9 percent in zone 3, which are similar to the initial specific yield for all three zones of 5 percent (table 5).

The mean horizontal hydraulic conductivity for the coarse-grained alluvium, estimated at pilot points during the calibration process, was 181 m/d in the Iowa River valley (zone 4), 190 m/d in valleys of major tributaries of the Iowa River (zone 5), and 95 m/d underlying minor tributaries of the Iowa River (zone 6; fig. 14B; table 5). Calibrated horizontal hydraulic conductivity ranged from 10 to 1,000 m/d in zones 4 and 5, and 5 to 500 m/d at pilot point locations in zone 6 (fig. 14B; table 5). The mean calibrated vertical anisotropy ratios for the coarse-grained alluvium were 16 percent at pilot points in zones 4, 5, and 6 (fig. 15B), which were generally similar to the initial value of 10 percent (table 5). Mean calibrated vertical anisotropy ratios for the coarse-grained alluvium ranged from 1 to 100 percent in zones 4, 5, and 6 (fig. 15B; table 5). The mean calibrated specific storage was 6.53x10⁻⁵ m⁻¹ in zone 4, 6.43x10⁻⁵ m⁻¹ in zone 5, and 5.56x10⁻⁵ m⁻¹ in zone 6, which were similar to the initial specific storage of 3.45x10⁻⁵ m⁻¹ for all three zones. For specific yield of the coarse-grained alluvium, the mean calibrated values were 17 percent in zone 4, 26 percent in zone 5, and 23 percent in zone 6, which are similar to the initial specific yield for all three zones of 15 percent (table 5).

The mean horizontal hydraulic conductivity for the glacial drift (zone 7), estimated at pilot points during the calibration process, was 0.53 m/d and ranged from 0.03 to 3 m/d (fig. 14C; table 5). The mean calibrated vertical anisotropy ratio for the glacial drift was 18 percent, which was similar to the initial value of 10 percent (table 5) and ranged from 1 to 100 percent (fig. 15C; table 5). The mean calibrated specific storage was 4.85x10⁻⁷ m⁻¹, which was similar to the initial specific storage of 2.74x10⁻⁷ m⁻¹, and the mean calibrated specific yield was 5 percent, which was similar to the initial specific yield of 3 percent (table 5).

The mean horizontal hydraulic conductivity for the buried valley sediment (zone 8), estimated at pilot points during the calibration process, was 92 m/d and ranged from 5 to 500 m/d (fig. 14D; table 5). The mean calibrated vertical anisotropy ratio for the buried valley sediment was 20 percent, which was similar to the initial value of 10 percent (table 5) and ranged from 1 to 100 percent (fig. 15D; table 5). The mean calibrated specific storage for the buried valley sediment was 9.48x10⁻⁵ m⁻¹, which was similar to the initial specific storage of 5.25x10⁻⁵ m⁻¹, and the mean calibrated specific yield was 17 percent, which was similar to the initial specific yield of 10 percent (table 5).

The mean horizontal hydraulic conductivity for the Pennsylvanian shale (zone 9), higher hydraulic conductivity Mississippian limestone/dolomite (zone 10), Devonian shale (zone 11), Silurian-Devonian limestone/dolomite (zone 12), and lower hydraulic conductivity Mississippian limestone/dolomite (zone 13), estimated at pilot points during the calibration process, was 0.62 m/d, 26 m/d, 0.42 m/d, 11 m/d, and 0.77 m/d, respectively (fig. 14E; table 5). Calibrated horizontal hydraulic conductivity ranged from 0.07 to 3 m/d, 1.5 to 143 m/d, 0.03 to 2.5 m/d, 0.25 to 101 m/d, and 0.09 to 2.8 m/d at pilot point locations in zones 9, 10, 11, 12, and 13, respectively (fig. 14E; table 5). The mean calibrated vertical anisotropy ratio for the bedrock unit was 23 percent at pilot points in zone 9, 21 percent in zone 10, 16 percent in zone 11, 13 percent in zone 12, and 24 percent in zone 13 (fig. 15E), which were generally similar to the initial value of 10 percent (table 5). The vertical anisotropy ratio ranged from 5 to 62 percent at pilot point locations in zone 9, 1 to 100 percent in zones 10, 11, and 13, and 1 to 70 percent in zone 12 (table 5). The mean calibrated specific storage of the bedrock units were 9.74x10⁻⁷ m⁻¹ in zone 9, 4.40x10⁻⁷ m⁻¹ in zone 10, 4.24x10⁻⁷ m⁻¹ in zone 11, 3.75x10⁻⁷ m⁻¹ in zone 12, and 9.52x10⁻⁶ m⁻¹ in zone 13, which are similar to the initial specific storage of 2.74x10⁻⁷ m⁻¹ for zones 9–12 and 5.25x10⁻⁶ for zone 13. For specific yield of the bedrock units, the mean calibrated values were 1 percent in zone 9, 20 percent in zone 10, 1 percent in zone 11, 22 percent in zone 12, and 11 percent in zone 13, which are similar to the initial specific yield of 1 percent for zones 9-12 and 5 percent for zone 13 (table 5).

The calibrated streambed hydraulic conductivity for stream segments in the SFR package ranged from 0.003 to 3.05 m/d with a mean of 0.84 m/d. The mean calibrated streambed hydraulic conductivity was greater than the initial value of 0.003 m/d. The mean calibrated recharge multiplier was 1.08 and ranged from 0.33 to 1.67, the mean calibrated evapotranspiration multiplier was 1.01 and ranged from 0.33 to 1.67, and the mean calibrated well multiplier was 1.00 and ranged from 0.85 to 1.20.

Comparison of Simulated and Observed Water Table Altitudes

Numerical model performance was assessed primarily on its ability to match calibration targets for groundwater levels (hydraulic head altitudes) and changes in hydraulic head. Hydraulic head altitude targets were assessed for 738 wells (fig. 16), and hydraulic head change targets were assessed for 33 wells that had more than one measurement.

The hydraulic head altitude targets used in numerical model calibration were compared with numerical model-simulated hydraulic head calculated at the same location and time as the observations. Observed hydraulic head altitudes were plotted with simulated hydraulic head altitudes and compared to a 1:1 perfect-fit line (fig. 16). The linear regression coefficient of determination (R²) for the simulated versus observed hydraulic head altitudes in layers 1 and 2, which make up the Iowa River alluvial aquifer, was 0.92, and the R² value for the simulated versus observed hydraulic head altitudes within all layers was 0.56 (fig. 16). In general, simulated hydraulic head altitudes in layers 1 and 2 tended to underestimate observed values, and as observed hydraulic head altitudes increased, so did the difference between simulated and observed hydraulic head altitudes (fig. 16). For all model layers, simulated hydraulic head altitudes tended to overestimate observed hydraulic head altitudes for observed values less than about 270 m above NAVD 88 and underestimate hydraulic head altitudes for observed values greater than about 270 m above NAVD 88 (fig. 16).

The differences between observed and simulated hydraulic head altitudes, otherwise known as residuals, were plotted by layer using a histogram (fig. 17) to assess the overall model-to-measurement fit. The differences were calculated as observed minus simulated values, so a negative residual indicates that the simulated hydraulic head altitude was greater than the observed hydraulic head altitude at the given location. The histogram shows that, in general, hydraulic head altitude residuals for wells in layer 1 tended to be positive, residuals for wells in layer 5 tended to be negative, and residuals for wells in layers 2 and 3 tended to be more evenly distributed between positive and negative (fig. 17).

The absolute hydraulic head altitude residuals (which do not take into account whether simulated values were over- or underestimations of observed values) for all model layers ranged from 0 to 62.5 m (fig. 17), with a mean of 8.8 m and a standard deviation of 9.6 m. Absolute hydraulic head altitude residuals were evenly distributed spatially across the active model area. For layers 1 and 2, absolute hydraulic head residuals ranged from 0 to 57.5 m, with a mean of 5.7 m and standard deviation of 7.3 m. For all model layers, 46.5 percent of simulated groundwater levels were within plus or minus 5 m of the observed values, 68.6 percent were within plus or minus 10 m (with 11.1 percent being positive residuals and 11.1 percent being negative), and 93.4 percent were within plus or minus 25 m (with 11.1 percent being positive residuals and 13.7 percent being negative; fig. 17).

Observed and simulated hydraulic head hydrographs were created for select observation wells in layers 1 and 2 with more than three water-level observations for their period of record during the model time period (fig. 18). The hydrographs are a visual representation of hydraulic head altitude and hydraulic head change observation targets. The first value in a hydrograph represents the hydraulic head altitude at the time of the first measurement at an observation well. A visual assessment of the hydrographs indicated that in general, for wells located in layers 1 and 2, the numerical model underestimated observed hydraulic head altitudes and tended to underestimate the magnitude of the change in hydraulic head observed at wells. The timing of simulated lowest hydraulic head altitudes at wells generally matched the timing of observed lowest hydraulic head altitudes, and the timing of simulated highest hydraulic head altitudes generally matched the timing of observed highest hydraulic head altitudes (fig. 18).

Comparison of Streambed Fluxes

Numerical model performance was also assessed based on the differences (or change) in observed and simulated base flow. Streambed fluxes were calculated for each cell of the SFR package in the active model area during the period of record for each streamgage and were plotted as the change in streambed flux between streamgages in the active model area (fig. 19). When compared to observed streambed fluxes (fig. 19), the simulated streambed fluxes are underestimated, and in some cases simulated streambed flux accounted for one-fifth of observed streambed flux, typically during periods with higher streambed flux (fig. 19). In particular, the model accurately simulates the lower streambed flux values but does not accurately represent periods where there is an observed peak in streambed flux. This underestimation may be due, in part, to the fact that base flow was calibrated using only values for the winter months (October–March), when streamflow is generally lowest, although base flow tends to remain relatively steady during the year compared to streamflow. Additionally, inflows to the model area where the Iowa River and Linn Creek cross the active model boundary were determined from base flow estimated by the PRMS model, as described in the “Groundwater and Surface-Water Interactions” section of this report. These streambed flux values, calculated using the inflows to the model from the Iowa River or Linn Creek, may be underestimated and could lead to larger differences in observed and simulated streambed flux when compared to streambed fluxes calculated from the differences in observed and simulated base flow calculated using upstream and downstream streamgages.

Simulated Groundwater Budget

The simulated groundwater budget characterizes inflows and outflows to the Iowa River alluvial aquifer and underlying hydrogeologic units. Steady-state and transient groundwater budgets were calculated using the program ZONEBUDGET distributed with MODFLOW 6 (Langevin and others, 2017) using cell-by-cell flow data from the numerical model output. Groundwater withdrawals were not simulated during the steady-state stress period; therefore, the steady-state groundwater budget does not include a groundwater withdrawal component. The simulated steady-state groundwater budget for the model consisted of inflows from precipitation recharge (RCH package) and stream infiltration (SFR package) and outflows from ET_g (EVT package) and groundwater discharge to streams (SFR package). During the steady-state stress period, simulated inflow for the active model area was about 817,500 m/d from precipitation recharge and about 63,000 m/d from the infiltration of streamflow (table 6). Simulated outflow for the steady-state period was about 51,300 m/d from ET_g and 829,900 m/d from groundwater discharge to streams (table 6).

The groundwater budget for the transient period consisted of the same components from the steady-state period and also included inflows from and outflows to storage (STO package) and outflows from well withdrawals (WEL package). The groundwater budget for the transient period presented in this report was calculated as the mean of each monthly budget component for the entire transient period (January 1980–August 2022) and was expressed as percentages of the mean total monthly inflow and outflow. The percentage of total inflow or total outflow presented in this paragraph equal slightly more than 100 percent (100.1 percent) because of rounding. Based on the mean budget components calculated for the transient period, precipitation recharge accounts for about 60.6 percent of the total inflow, inflow from storage accounts for about 34.5 percent, and stream infiltration accounts for about 5.0 percent (table 7). For the mean outflows for the transient period, groundwater discharge to streams accounts for about 60.3 percent of the total outflow, outflow to storage accounts for about 34.5 percent, ET_g accounts for about 3.8 percent, and well withdrawals account for about 1.5 percent of total outflows (table 7).