Mining Features and Well Withdrawals

Two mining features, tailings basins and mine pits, use hydrologic stresses in addition to hydraulic-conductivity zones to model their effects on groundwater (fig. 9). Tailings basins are modeled using the MODFLOW General Head Boundary package (GHB) stresses (light blue, fig. 9; table 2) and using low hydraulic-conductivity zones in the uppermost model layer. The low hydraulic-conductivity zones account for the low hydraulic conductivity of the fine tailings sediments that coat the bottom of the tailings basins. While the low hydraulic-conductivity zone covered the entire area of each tailings basin, GHB hydraulic stresses (represented by a fixed basin stage [based on lidar elevation data] and a basin-lining conductance) are only in cells of the open-water part of each tailings basin. These GHB stresses represent the fact that mining operations constantly keep these pools full of water, potentially recharging groundwater.

Mine pits are modeled in detail using several techniques. Active mine pits are usually pumped to prevent flooding and allow mining, artificially lowering the water table to the mine-pit bottom (pumped mine pits, table 2). The pumping in active mine pits, which were generally pumped dry during the mining-model period, is represented using the MOFLOW Drain package (DRN; shades of pink to maroon, fig. 9). The elevation of the top of layer 1 is set to the mine-pit bottom, and drain cells allow all water discharging into mine pits to be removed from the model, keeping the mine pits dry. Drain stage is set to the minimum lidar elevation within a mine pit cell, representing the bottom of dry mine pits.

Inactive mine pits have been allowed to naturally fill completely with water or have their water level maintained somewhat lower by pumping or through a constructed or natural surface-water outlet (flooded mine pits, table 2). Inactive, flooded mine pits with no pumping or outlet have the top of model layer 1 set to the pool elevation (mean lidar elevation of the pool area), and model cells representing the water-filled volume of the mine pits are assigned a hydraulic conductivity of 1,121.044 ft/d in layer 1, 1,000 ft/d in layers 2–6, and 16.847 ft/d in layer 7. This choice allows any water entering a flooded mine pit to easily flow to any other cell in the same mine pit. Because depth data are unavailable, mine-pit bottom elevation is assumed to be 101 ft deeper (to the bottom of layer 6) than the pool elevation. No pool water level was imposed on flooded mine pits. Instead, MODFLOW calculated the water level in flooded mine pits based on the groundwater levels adjacent to the mine pits and water flows into, through, and out of the mine pits. The actual, measured pool elevation was used as a guide in calibration.

Inactive mine pits with pumped or surface-water outflows (lowered pools) are modeled with the same general approach used for mine pits with no pumping or surface-water outlets. In mines with lowered pools, mine-pit outflows are controlled by the stage of DRN hydrologic stresses placed throughout the area of the pit and in the appropriate model layer determined from the lidar elevation of the pool. The amount of flow from these DRN stresses are aggregated for each mine pit. These modeled flows are compared to flux targets in model calibration where measured mine outflows were reported in 2011 (A. Guertin, MN–DNR, written commun., 2019; R. Clark, Minnesota Pollution Control Agency, written commun., 2019), the same year that the lidar data used to determine the mine-pit pool elevations were collected.

Estimates of outflows from inactive mine pits during 2011 are available from the MN–DNR and Minnesota Pollution Control Agency (A. Guertin, MN–DNR, written commun., 2019; R. Clark, Minnesota Pollution Control Agency, written commun., 2019) for seven sites that do not discharge into streams with measured flow in the IRMA. Because these inactive mine-pit outflows cannot be used to help calibrate the model, the flows are represented using the MODFLOW Multi-Node Well 2 package (MNW2) (brown dots, fig. 9). Annual 2011 groundwater withdrawals from municipal, commercial, and industrial, nonmine-pit wells that require reporting to the MN–DNR (https://www.dnr.state.mn.us/waters/watermgmt_section/appropriations/wateruse.html, accessed April 2019) were also represented by the MNW2 package in the mining model (green-shaded dots, fig. 9). The MNW2 package was used for these features so that these features could be easily separated from lakes, which were modeled using the WEL package, in model mass-balance results.

All DRN hydrologic stresses (pumped mine pits, flooded mine pits with lowered pools, or mine pits with measured flows) were placed in cells of bedrock layers (5–8) because the mine pits are always constructed in the Biwabik Iron Formation. For each mine pit, DRN hydraulic stresses were assigned to cells in the bedrock layer containing the minimum lidar elevation of that mine pit. Because glacial sediments have been removed from mine pits, the thickness of these upper four layers (1–4) over mine pits is set to 0.25 ft each. Any model layer representing glacial material over a mine pit was given a very high hydraulic conductivity (see description of flooded mine pits, mentioned earlier) to ensure that water discharging to mine pits laterally from those layers can easily flow to DRN hydrologic stresses in the bedrock toward the bottom of the mine pit.

Waste-rock piles and mining-disturbed areas are represented as separate zones of moderately high hydraulic conductivity with no additional hydrologic stresses. These piles sit above the general land surface and consist of large pieces of broken hard bedrock, which justifies the moderately high hydraulic conductivity.

Observations

PEST calibrated the mining model to measured groundwater levels, stream base flow, mine-pit dewatering flows, and the degree to which the model correctly simulated wetland and nonwetland areas (discussed later). Of the 2,766 measurements forming the water level calibration dataset (fig. 11), 2,629 groundwater levels came from the Minnesota Well Index (T. Wahl, MGS, written commun., 2018, publicly available data at https://www.health.state.mn.us/communities/environment/water/mwi/index.html), which is maintained by the Minnesota Department of Health; 55 were supplied by the MN–DNR (2016); and 8 were supplied by the Minnesota Pollution Control Agency (A. Streitz, written commun., 2018). All calibration water levels are available in the model archive (Cowdery and others, 2023). The remaining 74 groundwater levels were from mine pits and mine-pit lakes determined from 2011 lidar elevation data. Most groundwater levels were single measurements (table 3) collected during drilling or development of wells for domestic, municipal water supply, and mining exploration uses. Average water levels were calculated at wells with more than one water level. Only levels from wells where the aquifer was known were included. Levels used for calibration were measured during 1900 to 2017, with 97 percent of those measurements being collected after 1970. Most measurements are from the 1990s and 2000s (25 percent and 30 percent, respectively). Fifty-seven percent of the water level observations were measured during the steady-state modeling period of 1995–2015, but during 2011, the year corresponding to lidar measurements used in determining lake and mine-pit water levels, only 3 percent (83 water levels) were available.

Base flow in streams (fig. 11) measured by the USGS and MN–DNR falls into three groups when used for calibration: base flow calculated at continuous-flow streamgages (7 sites), miscellaneous single streamflow measurements (6 sites), and streamflow from a group of sites (34 sites) at which one-time, concurrent synoptic flow was measured (table 4). All flow measurements were modified to model units (cubic feet per day [ft³/d]) and to better represent the steady-state period that the calibrated mining model simulates. A base flow estimate was calculated from continuous streamflow records at seven streamgages by hydrograph separation using the Institute of Hydrology Base Flow Index method (Institute of Hydrology, 1980; Wahl and Wahl, 1988). During 2012–15, four of these sites had complete daily average flow records. Three other sites had substantial daily flow records, but these records were incomplete during 2012–15. The incomplete sections of the records were estimated by calculating an average flow factor for the streamgage with missing records during periods of overlapping flow record with the USGS streamgage (USGS site 04015438) in the St. Louis River headwaters at Skibo, Minnesota (average daily flow at the Skibo streamgage during the period of overlapping record divided by the average daily flow at the Skibo streamgage during 2012–15). This average flow factor, multiplied by the average streamflow at the Skibo streamgage during missing streamflow periods, produced a complete average daily streamflow estimate for the streamgage with missing flow record for the 2012–15 period. Average daily base flow discharge, which was assumed to be groundwater discharge, was calculated at these seven sites during 2012–15.

Six miscellaneous, single flow measurements stored in the USGS National Water Information System database (USGS, 2022) and collected during 2012–15 were also used as flow observations in calibration. These measurements were corrected to produce 2012–15 representative flows using a factor of the flow at the Skibo streamgage on the measurement day divided by the average base flow at the Skibo streamgage during 2012–15. One set of synoptic streamflow measurements was collected for this study at 34 locations on August 15–17, 2018. The 34 synoptic measurements were also adjusted to produce 2012–15 representative flows with the same method used to adjust the miscellaneous flow measurements.

The flows of rejected recharge and surface seepage calculated by the UZF package are important components of the water budget in the IRMA, and some part of these flows likely contributes to base flow measured and used as flow calibration targets in the mining model. The rejected recharge and surface seepage that do not find their way to streams are eventually lost to ET. To accurately calibrate the model, the flows not lost to ET were added to model-produced stream base flows using a post-MODFLOW processing utility program because they are not routed within the implementation of the UZF package used. As noted previously, 20 percent of the rejected recharge was assumed to not leave the Basin through ET and was added to the appropriate stream base flow. In contrast, 80 percent of surface seepage was assumed to not leave the Basin through ET and was added to the appropriate stream base flow

Reported dewatering flow at 16 mine pits was also used as calibration flow targets. The aggregated outflow of DRN hydrologic stresses for each mine pit, which is used to maintain mine-pit water levels in the model, was compared to the average reported mine-pit dewatering flow. Aggregated dewatering flows were routed (as part of the postprocessing) to the appropriate streamgages to augment simulated stream base flows. Two mine pits on the western end of the IRMA were within the Mississippi River drainage area where streamflows were not measured. The discharge from these mine pits was not routed.

Total area of mapped permanent wetland in the IRMA was used as an observation (brown areas, fig. 10) and compared to total permanent wetland area produced by the model (area of groundwater-surface seepage. See definition of model-produced wetland area in “Wetland” section presented previously). The difference in these areas multiplied by a weighting factor of 0.004 (table 4) was added to the objective function, which was minimized during PEST calibration. Two types of wetland-area errors are possible. False-negative wetland areas are those where mapped wetlands exist, but the model does not produce surface seepage. False-positive wetland areas are those where there are no mapped wetlands, but the model produces surface seepage.

Weights

Weights were applied to the water levels to control their importance in parameter estimation. Calibration water levels were placed in the model layer corresponding to the depth of the well open interval. For wells where the geologic material of the open interval did not match that of the containing model layer, observation locations were moved vertically by as much as 80 ft to an adjacent cell that had matching geologic material. Observations that were moved received lower weight in calibration; water levels from wells with an open interval in geologic material that matched that of the corresponding model cell have more weight (number of water levels matching or not matching for each type of material is shown in table 5). Water levels from wells that are within three model cells of a boundary condition (including hydrologic stresses) also have less weight because specified fluxes at the lateral boundaries can produce modeled groundwater heads that are substantially different from measured heads. The weights used for calibration water levels are shown in table 4.

The weights applied to base flow observations favored values calculated from continuous hydrographs and those with higher flows (table 4). Weights were calculated by taking the appropriate weight dividend from table 4 and dividing it by the base flow to produce the weight used for that base flow. Average base flows calculated from continuous hydrographs have high weight dividends that produce relatively high weights because they better represent the average flow conditions that a steady-state model produces. Single base flow measurements have relatively low weight dividends that produce relatively low weights because these measurements are idiosyncratic of the flow conditions at the time of measurements and may not represent the long-term average conditions produced by a steady-state model. High base flow had higher weight dividends because higher flows were measured at locations that have a larger drainage area, thus integrating a larger part of the model. The choice of these weights ensures that the overall groundwater flows through the model have more effect on parameter values than do flows from a small area of the model.

Mine-pit discharge measurements were all given a relatively low weight of 20 because they have a large uncertainty and are all about the same magnitude. The weight for wetland false-positive area and false-negative area observations were 0.004, which is small, because the method used to determine the wetland area in a steady-state condition in the IRMA is uncertain, as described previously. Further, the method used to determine what areas of the steady-state model are wetlands also contains many assumptions and uncertainties.

Parameters

Model hydraulic conductivity was initially calibrated using zones (layers 1, 7, and 8) and pilot points (layers 2–6; Doherty, 2003). Pilot points is a parameterization method where values of the parameter are estimated at modeler-specified locations in the model. The model hydraulic-conductivity arrays were interpolated from the hydraulic conductivity at pilot-points using modeler-specified estimation choices. This approach offers a compromise between the computationally burdensome option of specifying a parameter in every model cell and the rigid representation of specifying a parameter in a small number of constant hydraulic-conductivity zones (Hunt and others, 2007; Anderson and others, 2015). Singular-value decomposition (Kalman, 1996; Aster and others, 2013) was used to ensure stability in the parameter-estimation process using settings suggested by Doherty and Hunt (2010). Tikhonov regularization (Tikhonov, 1963a; Tikhonov, 1963b; Doherty, 2003; Doherty, 2010) was used during parameter estimation to prevent overfitting (that is, unrealistic parameter values) following the approach of Doherty and Hunt (2010). See Doherty and Hunt (2010), Anderson and others (2015), and Leaf and others (2015) for a detailed description of the application of these methods.

The PEST calibration adjusted 306 parameters: 21 hydraulic-conductivity zones, 252 pilot points in 12 hydraulic-conductivity zones, 29 vertical hydraulic conductivity anisotropy zones, recharge multipliers in 2 zones, and streambed and drain conductance (tables 6 and 7). Pilot points were not used in layer 1 because the available detailed surficial geology mapping was able to provide appropriately heterogeneous hydraulic-conductivity distributions using zones of mapped sediment texture. The glacial geology detail produced for layers 2–4 was much lower than in layer 1 because it was interpolated from a few hundred well logs that had the range of sediment types aggregated into six material groups. Two of these materials (St. Louis Sublobe fine-grained and course-grained sediments) do not appear in layers 2–4. The remaining four glacial materials were treated as single material zones regardless of layer, and their hydraulic conductivity was calibrated using pilot points within those zones. The four bedrock materials in layers 5–8 were given separate hydraulic-conductivity zones in each layer to account for decreased fracturing with depth. Pilot points were placed within the zones in layers 5 and 6 (the uppermost two bedrock layers) because there were enough groundwater head measurements in these layers to inform hydraulic-conductivity variability. Pilot points were not placed in layers 7 and 8, because so few groundwater head measurements were available to constrain pilot point parameters. Piecewise-constant hydraulic-conductivity zones represented the four bedrock materials in layers 7 and 8.

The initial hydraulic conductivity for each zone was set within ranges reported from the literature to maintain an appropriate relative ranking of hydraulic conductivity of aquifer materials to one another (table 6). The hydraulic conductivity of Des Moines glacial fine-grained sediments is less than that of the Rainy Lobe, which is less than that of the Superior Lobe, based on the general texture of those tills (Wright, 1972; Johnson and others, 2016). Relations in the variability of hydraulic conductivity in bedrock were incorporated into the model by differing ranges in parameter constraints (table 5). The order of decreasing bedrock hydraulic conductivity is the Biwabik Iron Formation (which also is more highly variable), the Duluth Complex, the Virginia formation, and crystalline bedrock (Barr Engineering Company, 2014). Most upper and lower bounds for pilot-point hydraulic-conductivity parameters spanned two orders of magnitude for each zone, where the initial value was at or near the mean of its range. Zones were given wider bounds to reflect the range of reported literature values and expected upscaling artifacts that may result from using literature values in a groundwater-flow model (table 6). Hydraulic conductivity was initially free to vary in each unit independently during parameter estimation. Optimum hydraulic conductivities were evaluated for reasonableness at the end of the parameter estimation. Vertical anisotropy for each initial hydraulic-conductivity zone was calibrated as a separate parameter (table 7) regardless of whether that zone had pilot points. Using an anisotropy parameter with a lower bound of 1.0 ensured that optimal parameters always had the expected relative relation of horizontal hydraulic-conductivity values being equal or greater to the vertical hydraulic-conductivity value.

The potential areal recharge array was divided into two areas, each with a calibration multiplier parameter: “wetland” recharge areas of low slope and “upland” recharge areas of higher slope (table 7, fig. 10). For the purposes of calibration, wetland areas are those defined as peat accumulating permanent wetlands (see “Wetland” section). All other areas are considered uplands. The two wetland recharge area multipliers retain the relative potential-recharge distribution produced by the SWB model across the IRMA but allow PEST to calibrate the mining model to more closely simulate the correct total flux through the aquifers in the IRMA. The multiplier for upland areas was allowed to vary over a wider range during calibration because the uplands are more variable geologically, hydrologically, and biologically than are the wetland areas.

Water Balance

Total groundwater flow through the mining model is nearly 171 million ft³/d. In a steady-state groundwater model, inflows should balance outflows. In this model, the balance is 0.012 percent. The components of the mining-model water balance (table 8) that have the largest flows are areal recharge (in, 78 percent), areal groundwater seepage to the surface (surface seepage out, 45 percent), seepage to streams (out, 43 percent), and stream seepage to groundwater (in, 12 percent). These four flows comprise 90 percent of the model inflows and 88 percent of the model outflows. Flows in other components (for example, inflows and outflows from the model perimeter, lakes, and mine pits) are each less than 5 percent of the total groundwater flow. Simulated model surface seepage occurs throughout the IRMA but is concentrated in the lowlands. Surface seepage rates in areas of extensive wetlands were small (less than 0.01 ft/d per cell), while surface seepage rates occurring near streams were larger.

Total groundwater flow out to all mining features and pumping wells (man-made features) is 4.1 percent of all groundwater flow. This is partly counterbalanced by groundwater recharge (flow in) from tailings basins of 2.7 percent, leaving net groundwater flow out to man-made features at 1.4 percent. Most of this outflow is mine-pit (dewatering) pumping (2.7 percent), other (nondewatering) mine pumping (0.6 percent), and groundwater discharge to tailings basins (0.4 percent). Coincidentally, the amount of groundwater leaving through mine pits is equal to the amount of groundwater entering from tailings basins (table 8), although some groundwater also flows out to tailings basins as seepage (0.4 percent). The model was only able to produce 55 percent of the specified amount of municipal and industrial groundwater withdrawals without model cells becoming dry. Similarly, only 68 percent of the lateral-boundary flows out of the model could be produced. Together, these unproducible flows are 2.1 percent of total groundwater flow through the model.

Groundwater flow decreases with depth in the model (tables 9–11). More than 80 percent of groundwater flow in the IRMA model occurs in the upper four model layers of surficial and glacial materials; flow between layers (downward, upward, or net flows) decreases by between one-half and one order of magnitude per model layer downward. This pattern reflects the fact that most of the hydrologic stresses (water sources and sinks) and aquifer transmissivity in the model occur in the uppermost layers (recharge, streams, and lakes for example, typically occur in the uppermost four layers). The distribution of this flow in the unconsolidated layers (1–4) likely is not evenly distributed spatially, however. In areas like the northwestern and northeastern part of the IRMA, where these layers are thin (white areas, fig. 2), considerably more groundwater may be flowing in the upper bedrock layers than in other parts of the IRMA. In the bedrock layers (5–8), mine-pit drainage and well pumping create hydraulic gradients that likely result in higher flows than would otherwise occur in the absence of such hydrologic stresses.

Head Calibration Results

Calibration of the mining model produced a good fit for heads (coefficient of determination, R²=0.93; fig. 15). The mean head residual (measured minus simulated) for all heads was −1.89 ft, indicating that modeled heads are biased slightly low. The mean absolute error residual for all heads was 11.03 ft, with a standard deviation of 14.68 ft. (root mean squared error:18.36 ft.). The mean absolute error of heads measured in 6 of 74 wells in mines and mine pits was more than twice as high as those from other head groups, reflecting the inability of the model’s discretization to capture all hydrologic variability and heterogeneity near mining activities.

The difference between measured and simulated heads (modeled head errors) average near zero (−1.89 ft), indicating that the model produces heads that are not substantially biased. Most of large head errors cluster near the mining features of the Iron Range, however, and these errors tend to be positive, indicating that the calibrated model tends to produce heads in this area that are too low. The sign and magnitude of the head errors within each of the seven head groups shown in figures 15 and 16 also appear random, indicating that the calibrated heads are not biased among these groups. Some simulated heads differed substantially (28 heads, 1 percent, greater than 71 ft, absolute value) from the measured heads, however, but are not in any particular group.

Several features of the calibration heads are responsible for the disagreement between measured and modeled heads. The measured heads used to calibrate the model were collected over the last 120 years during times that do not represent a steady state, especially during droughts and wet periods. Further, many of the heads were collected by drillers and reported as water levels below land surface even though measurements are often made from a location higher than land surface. The location of some wells and, therefore, the land-surface elevation derived from DEMs is poorly known, which can affect the elevation of water levels. Because the discretization of the model grid and the geology it contains are crude approximations at the scale of a well, the model geologic material of a cell containing the open intervals of some wells is not the same material as that on the stratigraphic record of the well in the real world. If a model well is finished in a clayey till but the real well is finished in a buried sand, a substantial difference between modeled and measured water level can be expected. Finally, the steady state produced by the model cannot exist in the natural world, so modeled steady-state heads will never perfectly match measured heads.

Although the model generally simulated water levels well, water levels from mine pits and water levels near other mining features have larger residuals than water levels in areas with geologic materials. Part of the error in simulated mine-pit water levels likely is because of the heterogeneous hydraulic conductivity of the Biwabik Iron Formation and the overlying coarse-grained glacial deposits that supply the mine pits with water. The density of measured water levels and, therefore, pilot points used in the model are inadequate to capture all the geologic variability in the IRMA. Further, mine-pit water levels were determined from the 2011 lidar DEM. This year was drier than average and did not represent the long-term steady-state heads that the model produces. The large head residuals near nonmine-pit mining features likely come from their proximity to mine pits, which induce large head gradients that change constantly as mine pits are developed and dewatered. Finally, the locations of bedrock wells (and their water levels) are likely skewed toward parts of bedrock that locally have higher hydraulic conductivity. Well drillers likely select these locations for wells because they produce more water. Because the model calibration is dependent on measured water levels from these higher hydraulic-conductivity locations in the bedrock, the final bedrock hydraulic conductivities in the model likely are closer to the real average bedrock hydraulic conductivity in the IRMA, not the higher hydraulic-conductivity locations of wells. This would result in higher head errors near productive wells.

Head errors in bedrock are also generally very large (more than 100 ft) near mines, and nearly all of the residuals greater than 30 ft occur along the Iron Range (fig.16). A disproportionate number of these large residuals occur in the Biwabik Iron Formation for the reasons outlined previously: the formation’s hydraulic conductivity is heterogeneous, and many water levels come from mine pits, which are themselves heterogeneous. While 6.5 percent of heads errors were greater than 30 ft, 57 percent of these are in the Biwabik Iron Formation. Other bedrock zones had fewer than 5 percent of these large head errors.

Flow Calibration Results

Simulated base flow at streamgages also match measured streamflow well (R²=0.97 or greater of all but miscellaneous base-flow measurements; figs. 17 and 18). The model’s ability to reproduce flow is affected by several factors. Average flow at continuous streamgages is based on long periods of many measurements and better characterizes steady-state base flow from the modeled period (1995–2015) than do single synoptic or miscellaneous measurements of streamflow. Flow errors are biased particularly high (above the 1-to-1 line in fig. 17) at streams with smaller flows (less than 25 cubic feet per second). The large departure from the 1-to-1 line is an artifact of the logarithmic scale of the plot. A given error in terms of flow (for example, 10 cubic feet per second, rather than percent) will appear large for measurements with low flows.

Simulated mine-pit flows have negative flow errors (below the 1-to-1 line in fig. 17). This is likely because measured mine-pit flows were given relatively low weight during calibration to reflect uncertainties in these observations. Mine-pit flows used in calibration were selected from 2011 to better correspond to the mine-pit water elevations, which were estimated from the lidar elevation data collected in 2011. But 2011 was abnormally dry, so these mine pit flows likely represent average conditions poorly. The low calibration weight of mine-pit flows means that they had little effect over the calibrated values of parameters.

Wetland Presence/Absence

The calibrated mining model was able to reproduce 71 percent of the mapped permanent wetland areas in the IRMA (white area divided by white plus orange area on fig. 19). The model produced 997 mi² of wetlands (white plus green areas on fig. 19, defined as areas where the water table is within 1 ft of land surface), 71 percent of which overlap with mapped permanent wetlands (fig. 19 white areas equaling 658 mi²). This amount is 108 percent of actual wetland in the IRMA. The model produced 339 mi² of permanent wetlands where they are not mapped (green areas on fig. 19; false-positive areas). Most of the areas where the model incorrectly produced wetlands are in the southeast part of the St. Louis River Basin in areas between the mapped wetlands and rivers. This suggests that the water table is too high near rivers, perhaps indicating that the hydraulic conductivity of surficial materials is too high or perhaps that these riparian wetland areas are not correctly mapped as wetlands, possibly because they are seasonal or because they do not support the plants used to identify wetlands. Large areas where the model was unable to produce wetlands occur in the northern part of the IRMA (orange area on fig. 19; false-negative areas). Surficial sediments are very thin (less than a few feet thick) over much of this area. The model likely has difficulty producing heads high enough to produce wetlands in bedrock near the land surface because the hydraulic conductivities of these rocks are calibrated in areas where they are buried beneath tens to hundreds of feet of glacial sediments. The PEST parameter for wetland area was effective for total wetland area produced by the model but not where those wetlands were produced.

Well Withdrawals

Municipal (including commercial) well withdrawals account for 0.2 percent and mining process (other mining) well withdrawals account for 0.4 percent of the total modeled groundwater flow in the IRMA. However, the model simulated only 38 percent of specified municipal pumping and 60 percent of specified other-mining pumping. Three wells simulated no water pumped at all because they were open to model cells that were dry. The MODFLOW version used for these simulations reduces pumping rates specified by the user if the aquifer cannot support the specified rate. Because the specified pumping rates are self-reported by water users, this discrepancy could be a result of inaccurately reported pumping rates. Alternatively, errors in model properties controlling flow in the vicinity of these wells could also result in simulated pumping being lower than reported pumping. Requiring that the model produce specified rates of withdrawal could be used as a calibration target for future modeling (for example, Parsen and others, 2016), if pumping were an important model objective for the model. Municipal, commercial, and industrial pumping is not the focus of this study and accounts for less than 1 percent of the total groundwater flow in the model.

Calibrated Parameters

Hydraulic-conductivity zones calibrated using pilot points produced zones of varying heterogeneity. The most heterogeneous bedrock zones were the Biwabik Iron Formation (hydraulic-conductivity range: 2.6–25 ft/d) and crystalline bedrock (hydraulic-conductivity range: 0.9–20 ft/d), in the uppermost bedrock model layer (layer 5; table 6). Fractures density and aperture are known to be greatest in the uppermost parts of bedrock in the IRMA and lessen with depth (M. Jirsa, MGS, written commun., January 30, 2018). Hydraulic-conductivity variability of the Biwabik Iron Formation is large because of high concentrations of vugs in some areas relating to its metasedimentary origins (Morey, 1972; M. Liljegren, MN–DNR, oral commun., 2018).

Although the calibrated model hydraulic conductivities generally simulated water levels well, hydraulic conductivities in mine-related zones likely are too high because their simulated heads are too low (yellow and red dots below the 1:1 line, fig. 15). Further, the hydraulic conductivity of the Biwabik Iron Formation and the overlying coarse-grained glacial deposits, which supply water to mine pits, are particularly heterogeneous. The density of pilot points in these hydraulic-conductivity zones likely are inadequate to capture their real-world variability.

In glacial sediments, the most heterogeneous hydraulic-conductivity zone was Rainy Lobe sands and gravels (range: 0.8–28 ft/d, table 6). Generally, water levels from wells screened in glacial sediments had small head errors (lie close to the 1:1 line, fig.15). Most large residuals in glacial sediments are negative (modeled less than measured head), suggesting that calibration was unable to produce hydraulic conductivities high enough in some areas (fig. 15). Each of the 29 hydraulic-conductivity zones had calibrated vertical hydraulic-conductivity anisotropy. The largest anisotropies are in tailings basins (167 times less than horizontal hydraulic conductivity) and St. Louis Sublobe sand and gravel (225 times less than horizontal hydraulic conductivity). There is no trend in vertical hydraulic-conductivity anisotropy with depth in bedrock, However, because horizontal hydraulic conductivity in bedrock decreases with depth by about an order of magnitude per model layer, and because vertical hydraulic-conductivity anisotropy varies less than an order of magnitude in each bedrock unit, calibrated vertical hydraulic conductivity also decreases with depth.

The potential areal recharge array was divided into upland areas and wetland areas, each with their own multiplier (fig. 10). The wetland-area potential-recharge array multiplier was 1.1 (75 percent of the upland-area multiplier), likely reflecting the large amount of rejected recharge from the high water table in wetland areas. In a sense, the wetland-area recharge array multiplier has the effect of controlling the amount of rejected recharge. Because 20 percent of this rejected recharge is added to streamflow, the wetland-area recharge multiplier has the effect of reducing model error for stream base flow measurements. The upland-area potential-recharge multiplier was 1.5, indicating that the SWB recharge estimates were too small to produce the amount of groundwater flow through the model needed for base flow in streams in that area. Recharge multipliers equal to 1.0 would represent no change from the initial SWB potential areal recharge values. These calibrated multipliers are not equal to 1 because the groundwater model contains more complete hydrologic processes than does the SWB model. Calibrated streambed and drainbed conductance were nearly the same as the initial conductance, indicating that initial values for these parameters were appropriate for the system modeled. The model is most sensitive to the potential areal recharge multipliers.

Model Scale, Discretization, and Parameterization

This model broadly describes conditions across 3,205 mi² of the IRMA, including the upper St. Louis River Basin along the Iron Range. Hydrogeologic materials were assumed to be horizontally isotropic, but with vertical anisotropy. Model cells are 400 ft×400 ft square and range in thickness from less than 1 to 800 ft. Cells of this size may be too large to accurately answer some question at a small scale, such as the details of groundwater interactions with surface waters. Like all finite-difference models, spatial input and model results (such as the lidar land-surface elevation, glacial stratigraphy, bedrock-surface elevation, groundwater heads, and flow paths near mine pits) are averages over the model cell area. Modeled groundwater heads, for example, may be shallower than the actual heads in places where groundwater gradients change over a small distance, such as near streams, wells, and mine pits. Near these features, modeled heads likely are too high. The bedrock-surface raster, which is coarser than land-surface DEM, also was averaged across each model cell, resulting in a model bedrock surface that was above the model land surface in some areas, particularly in areas where bedrock was at or near the land surface. To correct this, the model bedrock surface was adjusted to be no higher than 1 ft below the model land surface. Pilot points used to interpolate hydraulic-conductivity fields in many zones were widely spaced (at least 32,000 ft; 80 model cells) because of the few measured heads available. The resulting calibrated hydraulic-conductivity fields produced during calibration cannot represent actual hydraulic-conductivity variability. In hydrostratigraphic units that had very few observations, the hydraulic conductivity was calibrated as a zone, where one value is applied throughout the zone. In reality, hydraulic conductivity varies throughout the zone.

Groundwater flows preferentially through coarse-grained zones in unconsolidated sediments and through fractures and vugs in bedrock. While MODFLOW’s assumption of porous-medium flow is appropriate at the regional scale of the IRMA, this assumption is less appropriate at the small scale of a mine-pit wall, for example, where flow is likely dominated by fracture flow or flow through vugs. Analyses of model results that depend on groundwater-flow paths and flow times (well capture-zone delineation or groundwater travel times) likely are inappropriate where preferential groundwater-flow paths may be present. Because all model materials are assumed to be horizontally isotropic, modeled flow paths and flow times are likely inaccurate in materials that are anisotropically fractured or deposited.

Steady State

The results of the IRMA model represent unvarying average conditions and cannot provide information about the temporal dynamics of groundwater flow and heads. The potential areal recharge array specified in the model is an average during 1995–2015 (model steady-state period) and is the average condition represented in the mining model. Water levels and streamflows used in calibrating the model were measured during different periods. Calibration heads were measured during 1900–2017, with 97 percent of these measured after 1969 and 57 percent measured during the steady-state period. Calibration flows from continuous streamgages were averaged over a shorter period: 2012–15, corresponding to the years where the largest number of streamgages were active to provide a spatially consistent dataset. Flows measured outside this period were adjusted to represent flows during 2012–5. Precipitation during the calibration flows period was close to (6 percent higher than) that of the 1995–2015 average. Model results will be less accurate during times when general climatic and weather conditions are far different from those during the model steady-state period or if hydrologic temporal stationarity cannot be assumed (for example, under climate change scenarios).

Many hydrologic stresses (for example, mine-pit dewatering stage, tailings-basin pool stage) in the mining model required an elevation, which was determined from lidar DEMs measured in 2011. To ensure consistency with these elevations, municipal well and mine-pit dewatering withdrawal rates used in the model were those reported in 2011. However, 2011 was quite dry: precipitation at Aurora, Minnesota, near the center of the IRMA, was 78 percent (5.91 in.) lower than the 1995–2015 average. The use of 2011 lidar elevations resulted in hydrologic-stress stages that are too low for the steady-state period modeled. Model results are dependent on these stages in complex ways. For example, lower mine-pit DRN stage would tend to produce higher groundwater gradients to the mine pits, increasing mine-pit discharge. This is because mine-pit water levels, at least initially, will fall faster than adjacent groundwater levels, increasing discharge to the mine pits. This discharge routed to streams would produce higher simulated base flows for calibration. A model calibrated to these higher base flows may produce lower calibrated hydraulic conductivities in aquifers discharging to streams. The precipitation mismatch between the 1995–2015 model period and the 2011 hydrologic-stress stages should be considered when interpreting model results.

Hydrologic Stresses and Boundary Fluxes

The MODFLOW packages chosen to represent mining features, surface-water bodies, and wetlands in the mining model imply assumptions and create limitations for the results of the model. Dewatering from mine pits was represented using drains (DRN package), which can only remove water from aquifers. The amount of water removed from mine pits through drain stresses is controlled by the hydraulic gradient and the drain hydraulic conductivity, which is a calibrated parameter. Reported mine-pit dewatering flows were calibration targets and were aggregated and routed to downstream base flow targets. Reported flows from mine pits were assumed to reflect actual dewatering flow rather than permitted withdrawal rates and to lose no water in transit through evapotranspiration or other mechanisms on its path to a streamgage. Flows into and out of tailings basins are poorly constrained by calibration measurements. These flows were represented with the GHB package stresses, which define a groundwater head above which water will flow out of a model cell into a tailings basin and below which water will flow into a model cell out of a tailings basin. The assumption is that the tailings basin can supply an unlimited amount of water to the model cell in response to the downward gradient. This flow is controlled by the head gradient and the GHB conductance, which was estimated on the assumptions that the bed sediments of a tailings basin are fine grained and poorly compacted. The assumed conductance value corresponds to a hydraulic conductivity of 0.1 ft/d and a bed-sediment thickness of 1 ft. The conductance is uncalibrated because the flows into or out of tailings basins to groundwater are unmeasured. At this conductance, tailings basins supply 2.7 percent of flow into groundwater in the mining model.

Wetlands are not explicitly modeled using a hydrologic stress in the mining model because wetlands cover 29 percent of the IRMA, and to do so may over-constrain the modeled heads. Wetlands are assumed to be any area where groundwater levels are within 1 ft of land surface. If the UZF package had not been used to model wetlands, the elevation of the simulated water table would have been constrained in these areas. Thus, the groundwater head field would have been distorted to force groundwater discharge that actually terminates as seepage at wetlands to discharge to downgradient stream channels instead. Such groundwater discharge to streams in the large wetland areas of the IRMA would have artificially increased modeled streamflow and lowered aquifer hydraulic conductivity to prevent cell flooding.

Fractions of the water flux simulated as rejected recharge (20 percent) or as surface seepage (80 percent) were aggregated as runoff and routed to streamgages, augmenting modeled base flow for calibration purposes. These assumed runoff percentages were based on a water balance study in a wetland-rich area of northwestern Minnesota (Cowdery and others, 2019) and were not calibration parameters. Finally, the UZF variable SURFDEP was assumed to be 1 ft based on the average land-surface undulations in wetlands in the IRMA. These UZF package assumptions, though realistic, affected the calibration of the model and should be kept in mind when interpreting or using the model results. See Feinstein and others (2020) for a more complete discussion of possible effects of SURFDEP.

The potential-recharge array used in the UZF package was made from two SWB models of different resolutions: one coarse statewide model and one higher resolution model of the St. Louis River Basin. The St. Louis River Basin SWB model used streamflows from that Basin in its calibration and as a result, the potential-recharge array is likely most accurate in this area. The potential-recharge array was multiplied by two calibrated factors: a factor for upland areas and a factor for permanent wetland areas, which were assumed to be coincident with mapped permanent wetlands from a coverage from MGS for St. Louis and Cook County and from the MN–DNR’s NWI mapped permanent wetlands for Itasca County. The unsaturated zone of permanent wetland areas is assumed to behave fundamentally differently than the unsaturated zone of uplands because unsaturated zones in the uplands are substantially thicker than those in wetland areas. This assumption justifies separate potential-recharge array multipliers. Finally, areas of open water have zero potential-recharge in SWB models, under the assumption that any water for recharge in these areas comes from the surface-water body. These assumptions likely are valid for regional groundwater flow in aquifers in the IRMA. However, in specific areas and at local scales, some of these assumptions may not be accurate.

Well withdrawals were modeled in the mining model using the MODFLOW MNW2 package. Lateral-boundary groundwater fluxes are supplied by the WEL package. Both packages will inject or extract the specified amount of water into the model at their cell locations. However, if these packages are specified to remove water from the model, they can only do so if the cells in which they are located retain some water. If a well stress removes so much water that a cell goes dry, the actual withdrawal specified in the model input file is reduced until the cell retains water. Collectively, MNW2 stresses could only produce 55 percent of the specified well withdrawals without cells drying out, although the withdrawal produced is variable among wells. This means that the modeled aquifers with calibrated hydraulic conductivities have lower hydraulic conductivities in the vicinity of wells than do the actual aquifers. Therefore, this model would be inappropriate to answer specific questions about specific wells. Lateral-boundary WEL stresses could only produce (remove) 68 percent of specified water without going dry. This means that the modeled aquifers could not transmit as much water in some areas compared to the real aquifers, implying that hydraulic conductivities are too low in modeled aquifers in those areas. For example, along the northwestern boundary of the model, streams become losing (red circles) at the model edge (fig. 14). In these cells, the aquifer cannot supply the specified WEL withdrawals, so water is induced to flow out of streams to supply some of the water. Model results in the vicinity of these unsatisfied WEL stresses are likely not reliable. However, because MNW2 and WEL stresses account for less than 5 percent of model flows, these problems likely have a small effect on the regional results of the mining model as a whole. The boundary flows were derived from the results of a 1-layer analytic-element model`. Flow differences can be expected between this simple 1-layer model and the complex, multilayer IRMA model.

Representativeness of Heads Used for Calibration

Not only are model calibration heads limited temporally, they are also limited by their spatial extent and ability to represent all geologic material in the IRMA. Wells and the data that they provide (hydraulic conductivities from aquifer/specific capacity tests, measured groundwater heads) are not evenly distributed spatially or among aquifers. Wells are drilled where water is needed, available, or where water levels or samples are needed. Water is needed near residential, commercial, and industrial development (mainly mines and surrounding infrastructure in the IRMA). Water is available in more productive aquifers or the more productive parts of them. More transmissive parts of hydrologically heterogeneous aquifers, like the Biwabik Iron Formation, are more likely to have water wells drilled in them. This bias likely results in hydraulic-conductivity estimates being higher than the average part of the aquifer actually is. Water levels from wells are then skewed to the more productive parts of aquifers where water and water data are needed. Additionally, water levels can be skewed by imprecise measurement. Nearly all of the water levels used in model calibration were measured by well drillers when the well was drilled. These measurements are often poorly located, rounded to the foot, and often do not account for the measuring point height above the ground, thus skewing water-level depths (R. Barnes, University of Minnesota and others, written commun., 2019).

Even with the limitations of the mining model, regional analyses of model results are valid and appropriate. An example of appropriate analysis is regional comparisons of the spatial patterns of groundwater flows and heads in the mining-model with those patterns in the pre-mining model (described below). However, comparisons of these results in specific areas, especially near mining features, are inappropriate because of the model limitations described previously. Another example of appropriate use of model results is using mining-model groundwater flows as specified boundary flows in another model inset in the IRMA. Care must be exercised when constructing inset models to ensure that the effects of hydrologic stresses within the inset model do not propagate to the specified flow boundaries. Other uses of the mining model are inappropriate given the model limitations. In particular, using the model to evaluate site-specific hydrologic stresses (like a new mine pit) may not produce valid results because the scale of the added stress is too small in a regional model. Such scenarios are more appropriately explored with an inset model bounded by flow from the mining model.