Groundwater/surface-water interactions in the Partridge River Basin and evaluation of hypothetical future mine pits, Minnesota

<p>The Partridge River Basin (PRB) covers 156 square miles in northeastern Minnesota with headwaters in the Mesabi Iron Range. The basin is characterized by…

Public-domain full text preserved in the Mountain Man Mining Library. Original source: pubs.usgs.gov.

U.S. Department of the Interior U.S. Geological Survey Scientific Investigations Report 2021–5038 Prepared in cooperation with the Prepared in cooperation with the Great Lakes Indian Fish & Wildlife Commission Great Lakes Indian Fish & Wildlife Commission Groundwater/Surface-Water Interactions Groundwater/Surface-Water Interactions in the Partridge River Basin and Evaluation in the Partridge River Basin and Evaluation of Hypothetical Future Mine Pits, Minnesota of Hypothetical Future Mine Pits, Minnesota

Cover—Photograph of Partridg River near Hoyt Lakes, Minnesota. Photo by Megan Haserodt, U.S. Geological Survey, 2018

Groundwater/Surface-Water Interactions in the Partridge River Basin and Evaluation of Hypothetical Future Mine Pits, Minnesota By Megan J. Haserodt, Randall J. Hunt, Michael N. Fienen, and Daniel T. Feinstein Prepared in cooperation with the Great Lakes Indian Fish & Wildlife Commission Scientific Investigations Report 2021–5038 U.S. Department of the Interior U.S. Geological Survey

U.S. Geological Survey, Reston, Virginia: 2021 For more information on the USGS—the Federal source for science about the Earth, its natural and living resources, natural hazards, and the environment—visit ://www.usgs.gov or call 1–888–ASK–USGS. For an overview of USGS information products, including maps, imagery, and publications, visit ://store.usgs.gov/. Any use of trade, firm, or product names is for descriptive purposes only and does not imply endorsement by the U.S. Government. Although this information product, for the most part, is in the public domain, it also may contain copyrighted materials as noted in the text. Permission to reproduce copyrighted items must be secured from the copyright owner. Suggested citation: Haserodt, M.J., Hunt, R.J., Fienen, M.N., and Feinstein, D.T., 2021, Groundwater/surface-water interactions in the Partridge River Basin and evaluation of hypothetical future mine pits, Minnesota: U.S. Geological Survey Scientific Investigations Report 2021–5038, 94 p., ://doi.org/ 10.3133/ sir20215038. Associated data for this publication: Haserodt, M.J., Hunt, R.J., Fienen, M.N., and Feinstein, D.J., 2021, MODFLOW-NWT and MODPATH models, capture zones and uncertainty data analysis for the Partridge River Basin, Minnesota, U.S. Geological Survey data release, ://doi.org/10.5066/P9VODOU8. U.S. Geological Survey, 2019, USGS water data for the Nation: U.S. Geological Survey National Water Information System database, accessed March 2019, at ://doi.org/10.5066/F7P55KJN. ISSN 2328-0328 (online)

Acknowledgments The authors would like to gratefully acknowledge John Coleman of the Great Lakes Indian Fish & Wildlife Commission for providing extensive technical support with his vast knowledge of the study area. The authors would also like to acknowledge Anna Baker and Tim Cowdery of the U.S. Geological Survey in Minnesota for providing useful input on the regional geology, and this study built upon their groundwater-flow modeling work of the Mesabi Range. Finally, the authors acknowledge Todd Anderson and William Eldridge of the U.S. Geological Survey Dakota Water Science Center and Alden Provost of the U.S. Geological Survey Integrated Modeling and Prediction Division Earth Systems Modeling Branch for their thoughtful reviews.

Contents Boundary Conditions 16 Streams 25 Appendix 1. Additional Data Processing Steps to Build the MODFLOW-NWT

Figures

1. Map showing the location of the Partridge River Basin, northeastern

2.  Map showing the stream network and mapped wetlands in the Partridge

3.  Map showing known mineral deposits in the Duluth Complex and existing mine

4.  Figure showing the thickness of unconsolidated material over bedrock in the

5.  Figure showing lithology and glacial geology of unconsolidated deposits in the

6.  Figure showing peat deposits in the Partridge River Basin, northeastern

9.  Figure showing cross section of the Quaternary unconsolidated material and bedrock formations from north of the Mesabi Range to the St. Louis River in

10.  Figure showing A, Layer 1 thickness and mapped wetlands, B, Layer 2

11.  Figure showing cross section model layers including layers 1–3 in the glacial

12.  Figure showing boundary conditions presenting streams, lakes, groundwater fluxes along the model edges, pumped pits, pit transfers, and groundwater

13.  Figure showing SWB-estimated net infiltration for 2011–2013 mining conditions

14.  Figure showing layer 1 hydraulic conductivity zones and represented

15.  Figure showing layer 4 hydraulic conductivity zones and represented

16.  Figure showing location of horizontal and vertical hydraulic conductivity pilot

17.  Graph showing modeled and measured heads for the 1995–2015 average

18.  Graph showing modeled and measured heads for the 2011–2013 mining

19.  Graph showing comparison of modeled and measured streamflow for the 2011–2013 mining conditions groundwater flow model of the Partridge River

20.  Map showing spatial distribution of hydraulic head target residuals for the 1995–2015 average conditions groundwater-flow model of the Partridge

21.  Map showing spatial distribution of streamflow flux target residuals and stream accumulation, for the 2011–2013 mining conditions groundwater-flow model of

22.  Figure showing horizontal hydraulic conductivity in each of the 11 model layers

23.  Figure showing vertical anisotropy in each of the 11 model layers for the

24.  Graph showing calibrated hydraulic conductivity values and expected

25.  Graph showing calibrated hydraulic conductivity values and expected literature ranges for each model layer representing the bedrock in the Duluth Complex of Keweenan Supergroup, Virginia Formation, Biwabik Iron-Formation

26.  Graph showing mass balance for the 2011–2013 mining conditions

27.  Figure showing simulated water-table elevation for the 2011–2013 mining

28.  Map showing simulated losing and gaining reaches of the Partridge River for

29.  Figure showing, A, The spatial distribution of net infiltration estimated by soil-water-balance (SWB) package , B, the SWB calibrated net infiltration after model calibration, C, the UZF recharge (SWB calibrated net infiltration minus rejected recharge), and, D, simulated rejected recharge for the 2011–2013

30.  Figure showing the spatial distribution of simulated recharge, as it is applied to

31.  Graph showing A, how much available recharge is rejected by unsaturated zone flow , B, percentage of the rejected recharge is routed to streams, C, percentage of the unsaturated zone flow seepage is routed to streams, and, D, percentage that routed groundwater runoff contributes to total streamflow for the 2011–2013 mining conditions MODFLOW-NWT model of the Partrigde

32.  Figure showing groundwater seepage from unsaturated zone flow for the 2011–2013 mining conditions MODFLOW-NWT model of the Partridge River

33.  Figure showing hypothetical mine pits and model forecasts locations for 12

34.  Graph showing model forecast values and associated uncertainty for the pit inflow rates to the Peter Mitchell pit at final buildout without new mining and

35.  Graph showing model forecast values and associated uncertainty for the pit

36.  Graph showing model forecast values and associated uncertainty for the six streamflow forecasts in the final iron mining scenario, four new mine

38.  Figure showing probabilistic capture zones for each of the four hypothetical

39.  Figure showing data worth for predicting pit inflow to the scenario 3 pit for

Tables

1.  Summary of regional hydraulic conductivity and vertical anisotropy estimates

3.  Calibrated horizontal hydraulic conductivity and vertical anisotropy values for the unconsolidated and bedrock materials in the Partridge River Basin,

4. Summary of scenario names and features for the hypothetical mine pit scenarios62

5.  Model forecast values for the wetland forecasts and associated uncertainty in Conversion Factors U.S. customary units to International System of Units Multiply By To obtain Length foot (ft) meter (m) mile (mi) kilometer (km) Area acre hectare (ha) square mile (mi2) square kilometer (km2) Flow rate cubic foot per second (ft3/s) cubic meter per second (m3/s) inch per year (in/yr) millimeter per year (mm/yr) cubic foot per second (ft3/s) gallons per minute (gal/min) Hydraulic conductivity foot per day (ft/d) meter per day (m/d) Datum Vertical coordinate information is referenced to the North American Vertical Datum of 1988 (NAVD 88). Horizontal coordinate information is referenced to the North American Datum of 1983 (NAD 83), Zone 15N. Elevation, as used in this report, refers to distance above the vertical datum.

Abbreviations Ga billion years lidar light detection and ranging DRN MODFLOW drain package GHB MODFLOW General Head Boundary MDNR Minnesota Department of Natural Resources MNW2 MODFLOW multi-node well package PRB Partridge River Basin RIV MODFLOW river package SFR2 MODFLOW streamflow routing package, version 2 SLRB Saint Louis River Basin SWB Soil Water Balance UPW MODFLOW upstream weighting package USGS U.S. Geological Survey UZF MODFLOW Unsaturated zone flow package WEL MODFLOW well package

Groundwater/Surface-Water Interactions in the Partridge River Basin and Evaluation of Hypothetical Future Mine Pits, Minnesota By Megan J. Haserodt, Randall J. Hunt, Michael N. Fienen, and Daniel T. Feinstein Abstract The Partridge River Basin (PRB) covers 156 square miles in northeastern Minnesota with headwaters in the Mesabi Iron Range. The basin is characterized by extensive wetlands, lakes, and streams in poorly drained and often thin glacial material overlying Proterozoic bedrock. To better understand the interaction between these extensive surface water features and the groundwater system, a three-dimensional, steadystate, groundwater-flow model of the PRB was developed by the U.S. Geological Survey in cooperation with the Great Lakes Indian Fish & Wildlife Commission using the finitedifference computer code MODFLOW-NWT. The model simulates steady-state base flow in streams and groundwater interactions using the streamflow routing (SFR2) package. Existing mining features including tailings basins, stockpiles, pumped mine pits, and flooded mine pits were simulated using either high hydraulic conductivity zones or the drain (DRN) package. The unsaturated zone flow (UZF) package was used to better represent the groundwater system in areas with a high water table and for wetlands often associated with such areas. UZF typically is used to represent unsaturated zone processes but also can simulate the rejection of recharge and groundwater discharge to the land surface when the water table is near land surface. The steady-state model used data from the 2011 to 2013 period when 2011 high-resolution land surface (light detecting and ranging [lidar]) data were avail­ able that reflected land-surface and water elevations from mining activity in the basin. The parameter-estimation software suite PEST_HP was used to obtain a best fit of the modeled to measured groundwater levels, streamflow, pit inflow rates, and mapped peat deposits. The PEST calibra­ tion used the target residuals from two models with the same model parameters and targets from two separate periods: (1) a 1995–2015 calibration model, which provided a larger number of calibration targets, and (2) a 2011–2013 mining conditions model, which included calibration targets that reflected condi­ tions consistent with the modeled mine-workings topography. Calibration of the PRB model resulted in ranges of glacial horizontal hydraulic conductivity parameters that generally agreed with literature values and other models of the region. Horizontal hydraulic conductivity of the bedrock was higher in the upper bedrock layers where numerous and con­ tinuous fractures have been observed and lower in the deeper bedrock layers. Average basin-wide calibrated infiltration was 5.3 inches per year. An average of 4.6 inches per year of infiltration crosses the water table and becomes recharge and 0.7 inch per year is rejected by UZF due to saturated condi­ tions at the land surface. Simulated groundwater runoff (the sum of rejected recharge and groundwater seepage to the land surface) can either be routed to streams or removed from the model as evapotranspiration. The calibrated model indicates relatively shallow groundwater-flow paths dominating and approximately 50 percent of the stream base flow coming from groundwater runoff. The 2011–2013 mining conditions model was then used to develop five model scenarios simulating the response of the groundwater and surface-water system to potential hydrologic stress. The purpose of these mine pit scenarios is to present a possible workflow to quantify a model’s uncertainty for a given model forecast and serve as a possible guide for initial data collection that may improve a future model’s ability to make such a forecast. The scenarios included one scenario with the currently existing Peter Mitchell pit at final buildout and flooded to an elevation of 1,500 feet, and four scenarios with a hypothetical, new mine pit plus the flooded Peter Mitchell at final buildout. The five model scenarios were used to forecast streamflow at six locations in the PRB, pit inflow rates for the new mine pits and the flooded Peter Mitchell pit, and the average depth to water in 12 wetlands. A linear uncer­ tainty analysis was performed using information from the PEST calibration and tools in the PyEMU python package to assess model uncertainty propagation to the model forecasts. Streamflows generally were reduced with future mining and the greatest streamflow reductions occurred from the flooded Peter Mitchell Pit, probably due to its large size. Average depth to groundwater in wetlands was most affected the closer the wetland was to a new mine pit. Linear uncertainty methods were also used to evaluate data worth, which is the ability for potential new groundwater elevation observations to reduce the uncertainty in sce­ nario forecasts. Data worth was performed for a grid of new hydraulic head observations. Overall, areas with nonnegligible

2    Groundwater/Surface-Water Interactions in the Partridge River Basin and Hypothetical Future Mine Pits data worth generally corresponded to wetland areas with no groundwater seepage to land surface from UZF. These model behaviors indicated that the land-surface boundary condition simulated by the UZF package was pinning the groundwater elevations to the land surface in areas with groundwater seep­ age (33 percent of the 2011–2013 base conditions model) such that the sensitivity to new observations in these areas was minimal. Therefore, representing wetlands as boundary conditions minimized the usefulness of data worth calculations because wetland areas were present over a large part of the model domain. Probabilistic capture zones were estimated for each of the mines in the model scenarios. A capture zone represents the area contributing recharge to a model feature, like a well or a mine pit, and can be calculated by forward tracking particles from the water table. By using Monte Carlo techniques, it is possible to generate estimated capture zones that include the probability of recharge capture given the uncertainty present in the model. Monte Carlo techniques use randomly gener­ ated model parameter sets sampled from a plausible parameter range to create many possible realizations. The resulting cap­ ture zone arrays were calculated by tallying the total number of realizations in which a particle from a model cell was captured by the feature. Probabilities from the Monte Carlo runs ranged from 1 (captured in 100 percent of the runs) near the pits to 0 (captured in 0 percent of the runs) at the edges of the capture zone. Capture zones were not always spatially continuous; for example, the capture zone for the proposed mine pits south of the flooded Peter Mitchell pit was discon­ tinuous with capture surrounding the proposed mine pit and north of the flooded Peter Mitchell pit. This northern section represents deeper groundwater flow paths that originate in the topographic high, move under the flooded pit, and discharge into the proposed pit. This pattern of capture indicates the pos­ sibility of some deeper flow through the upper fractured bed­ rock when the shallow groundwater flow system is modified. These results underscore that future site-specific applications of the base condition model require the input of site-specific data and recalibration to focus on the site of interest. Introduction The Partridge River Basin (PRB; fig. 1) drains 156 square miles on the eastern end of the Mesabi Iron Range in north­ eastern Minnesota. The basin is defined to the north by a ridge formed from the Giants Range Granite batholith. Immediately south of the ridge is the Mesabi Iron Range, which contains the Proterozoic-age Biwabik Iron-Formation of Animikie Group (hereafter Biwabik Iron-Formation) where iron mining has occurred for more than a century. Thin glacial features and bedrock outcrops form the hydrologic divide between the PRB and the Dunka River to the east. Thicker glacial deposits form the hydrologic divide with the Saint (St.) Louis River to the south. The Partridge River flows into the St. Louis River through two outlets from Colby Lake (fig. 2). One outlet is through the adjoining Whitewater Lake and then into the Whitewater Lake outlet channel and the other is through the main stem of the Partridge River directly into the St. Louis River. Lake levels in Whitewater Lake have been regulated since the lake was impounded in the 1950s to create a water storage reservoir for mining operations (Minnesota Department of Natural Resources [MDNR], 2007). Anthropogenic changes in the PRB are mostly related to mining activities along the Mesabi Range. Urban development within the PRB includes the city of Hoyt Lakes and the city of Aurora, which is on the western basin boundary (fig. 1). Residential development outside these cities is sparse. Land in the PRB generally is owned by either local govern­ ments, the State of Minnesota, the Federal government, or private companies (St. Louis County, Minnesota, 2017). In undeveloped areas, forests are extensive with numerous wetlands and lakes. Mapped wetlands, many of which contain peat deposits (Dengler and Wagner, 2017; Dengler and others, 2017), cover a third of the basin (fig. 2). Iron ore has been mined from the Biwabik Iron-Formation of the Mesabi Range in the northern PRB since the late 1800s (Tieberg, 2020) and mining operations continued to be active during the periods simulated by this study. Existing mining features include mine pits—flooded pits where mining has occurred in the past and pumped pits where mining is currently occurring—stockpiles of waste rock and overburden, and infrastructure such as buildings and roads (fig. 3). A tailings basin used for ore processing is in the study area north of the basin divide. To the south of the Mesabi Range is the Duluth Complex of Keweenawan Supergroup (hereafter Duluth Complex), which contains undeveloped mineral deposits with copper, nickel, and platinum group metal resources (Minnesota Mining, 2019). Within the PRB, there are several potential mining projects (MNDR, 2016) targeting known mineral resources in the Duluth Complex (fig. 3). At the time of this publication (2021), none of the Duluth Complex mining developments were started, but at least one had begun the permitting process. The groundwater system in this region has been described in Siegel and Ericson (1980) for the PRB and by Barr Engineering Company (Barr; 2006a, 2006b, 2011, 2014a, 2014b) for the area immediately surrounding the proposed Poly Met mine pit (within the PRB). Hydrogeology of the Mesabi Range to the west of the PRB has been investigated by Tetra Tech (2014), for the whole St. Louis River Basin by Haserodt and others (2019), and for the Mesabi Range by T.K. Cowdery (USGS, written commun., 2020) and Alderdice and others (1992). The work described here complements this previous work by characterizing the groundwater flow system across the PRB at a resolution needed to quantify groundwa­ ter/surface-water interactions. Groundwater-flow models are a way to assess groundwater-flow systems and the interaction between groundwater and surface-water features by using existing information (maps of glacial and bedrock geology, aquifer

Introduction     3 92°6' 91°48' 91°30' 47°24' 48°0' 47°42' Base map data from , 2015 Cities from ArcGIS, 2020 States from ArcGIS, 2012 Partridge River Basin from U.S. Geological Survey, 2012 Lake Superior from Minnesota Department of Natural Resources, 2019a WISCONSIN MINNESOTA LAKE SUPERIOR CANADA Hoyt Lakes Babbitt Soudan Tower Ely Winton Aurora Biwabik Study area 8 MILES 8 KILOMETERS EXPLANATION Cities and towns Model domain Partridge River Basin Mesabi Iron Ridge Figure 1.  Location of the Partridge River Basin, northeastern Minnesota.

4    Groundwater/Surface-Water Interactions in the Partridge River Basin and Hypothetical Future Mine Pits 47°42' 92°12' 92°6' 92°0' 91°54' 91°48' 47°36' 47°30' 8 MILES 8 KILOMETERS Hillshade base from Minnesota Department of Natural Resources, 2011 Partridge River Basin from U.S. Geological Survey, 2012 Hydrologic features modified from Minnesota Department of Natural Resources, 2019a Wetlands modified from U.S. Fish and Wildlife Service, 2018 EXPLANATION Lakes Streams Wetlands Model domain Partridge River Basin

P a rt ri d g e R e r E m b a r r a ss Ri e r D u n ka R er B ir h R e r St. o ui s Ri ve r Bear Island Lake Birch Lake Big Lake Seven Beaver Lake Stone Lake Butterball Lake Round Lake Pine Lake Colby Lake Whitewater Lake Figure 2.  The stream network and mapped wetlands in the Partridge River Basin, northeastern Minnesota.

Introduction     5 47°30' 47°36' 47°42' Hillshade base from Minnesota Department of Natural Resources 2011 Mineral deposits from Minnesota Department of Natural Resources, 2016 Partridge River basin outline from U.S. Geological Survey, 2012 92°12' 92°6' 92°0' 91°54' 91°48' 8 MILES 8 KILOMETERS EXPLANATION Model domain Partridge River Basin Tailings basin Mine stockpiles Existing mine pits Known mineral deposits in the Duluth Complex of Keweenawan Supergroup Figure 3.  Known mineral deposits in the Duluth Complex of Keweenawan Supergroup and existing mine pits, stockpiles, and tailings basins in the Partridge River Basin, northeastern Minnesota.

6    Groundwater/Surface-Water Interactions in the Partridge River Basin and Hypothetical Future Mine Pits test data, land-surface elevations, streamflow measurements, groundwater elevations, surface water network geometry, and groundwater recharge data) with the well-developed mathematical equations for groundwater flow. The resulting groundwater-flow model is only as good as the data used to build it and the framework and assumptions used to incorpo­ rate these data into the model. Collection of additional data was beyond the scope of this study. The work here focused on using the existing understanding to construct science-based model forecasts of interest such as streamflow under base-flow conditions and depth to water in wetlands. In addition, the utility of future data collection efforts for forecasts of interest was evaluated. A steady-state, MODFLOW-NWT (Niswonger and others, 2011) groundwater-flow model was constructed for the PRB and surrounding areas. The model domain (shown in figure 1) was selected to focus on areas in the PRB with current (existing mine pits) and potential future development (known mineral deposits; fig. 3). The perimeter boundaries of the model were chosen to include competing hydrologic sinks (major rivers) outside the PRB so that PRB divides are simu­ lated and not specified. Perimeter boundaries were also located such that they were distant from mining areas and model forecasts to limit the effects of the model boundary conditions on the model forecasts. The purpose of this modeling effort is to quantify existing exchanges between the groundwater sys­ tem and surface-water features like streams and wetlands and to explore potential changes to these exchanges due to future stresses from additional mining development in the basin. The purpose of the hypothetical mine pit scenarios is to pres­ ent a possible workflow to quantify a model’s uncertainty for a given model forecast and serve as a possible guide for initial data collection that may improve a future model’s ability to make such a forecast. This report presents the results of the PRB modeling effort and includes: (1) an overview of the geologic history, (2) an overview of the hydrologic setting and site conceptual model, (3) discussion of modeling methods, (4) description of the construction and calibration of the groundwater flow model and mining scenarios, (5) exploration of model uncer­ tainty, probabilistic capture zones, and data worth, and (6) discussion of model assumptions and limitations. The main body of this report includes an overview of all methods used to build and calibrate the models, and five appendixes provide additional details on (1) data processing for the MODFLOW packages, (2) the method used to estimate the extent of dip­ ping bedrock units at depth, (3) the method used to process streamflow targets, (4) additional information on particle tracking with MODPATH, and (5) the setup and limitation for the data worth analysis. Geologic Setting The geology of the PRB consists of Proterozoic bedrock covered in thin, sometimes absent, Quaternary unconsolidated material. The thickness of the unconsolidated material ranges from absent to almost 200 feet (ft) thick based on a statewide depth-to-bedrock dataset (Minnesota Geological Survey, 2013; fig. 4) with much of the model domain covered by deposits that are less than 5-ft thick. Thick unconsolidated deposits are present in bedrock valleys in the north and south of the PRB model domain. Figure 5 shows the extent of the many glacial lobes across the model domain and the mapped sediment lithology (Jirsa, 2016a, b). Most of the PRB is covered by till and gla­ ciofluvial material from the Rainy Lobe. The southeast corner of the model domain has material from the Brainerd Sublobe of the Rainy Lobe and a few small areas to the east and south have material from the Superior Lobe. Much of the glacial deposits are classified as a diamicton, a poorly sorted till material with a wide range of particle sizes. Scattered sand and gravel deposits are present and areassociated with glaciofluvial depositional environments.Small areas scattered across the model domain are silty material. Overlying the glacial material are extensive peat deposits (fig. 6) that formed since the last glacial period (Jirsa, 2016a, b) and are important features for local groundwater hydrology (Bay, 1968). Near mining areas, the unconsolidated material is mapped as “fill/disturbed” (fig. 6) and likely includes some added materials in addition to remaining glacial material that was present before mining began. The bedrock in the region reflects a long and complex geologic history that began approximately 2.7 billion years ago (Ga). Figure 7 shows mapped bedrock units in the model domain. The oldest bedrock in the model domain are the Archean-age (approximately 2.7-2.6 Ga) Ely Greenstone, the Lake Vermilion Formation volcaniclastic rock, and the intrusive Giants Range Granite batholith (Card, 1990; Miller and others, 2001; Miller and others, 2002; Jirsa and others, 2011). Paleoproterozoic sedimentary strata of the Animikie Group were then deposited on the Archean greenstone and Giants Range Granite from ~2.1–1.8 Ga. Bedrock units in the model domain that are part of the Animikie Group include (from oldest to youngest) the Pokegama Quartzite (quartzite, siltstone, shale, and basal conglomerate), the Biwabik IronFormation (interbedded chert, iron silicates, magnetite, and hematite), and the Virginia Formation (interbedded shale, mudstone, and siltstone that is regionally equivalent to the Thomson Formation; Jirsa and others, 2005, Miller and oth­ ers, 2002). Deformation to the Animikie Group occurred ~1.9–1.7 Ga during the Penokean orogeny, a mountain building event that resulted in the 5˚–15˚ bedrock dips to the southeast and the uplifted ridge seen today in the Mesabi Range (Miller and others, 2002). Additional deformation and metamorphism occurred during the emplacement of the Duluth Complex and is responsible for dips of up 20˚–35˚ at the far eastern end of the Mesabi Range (Miller and others,

Geologic Setting    7 Hillshade base from Minnesota Department of Natural Resources 2011 Partridge River Basin from U.S. Geological Survey 2012 Thickness of unconsolidated material from Minnesota Geological Survey, 2013. 92°18' 92°12' 92°6' 92°0' 91°54' 91°48' 91°42' 47°24' 47°30' 47°42' 47°36' 47°48' 8 KILOMETERS 8 MILES EXPLANATION Thickness of unconsolidated material, in feet Model domain Partridge River Basin 165 to 180 180 to 195 195 to 210 210 to 225 75 to 90 90 to 05 105 to 20 120 to 135 135 to 150 1 50 to 165 5 10 to 25 25 to 60 60 to 75 Figure 4.  Thickness of unconsolidated material over bedrock in the Partridge River Basin, northeastern Minnesota.

8    Groundwater/Surface-Water Interactions in the Partridge River Basin and Hypothetical Future Mine Pits Quaternary geology modified from Jirsa, 2016b Hillshade base from Minnesota Department of Natural Resources, 2011 Partridge River Basin from U.S. Geological Survey, 2012 47°30' 47°36' 47°42' 92°12' 92°6' 92°0' 91°54' 91°48' 6 KILOMETERS 6 MILES Lithology Clay Diamicton Fill/Disturbed Gravel Sand Silt Glacial lobes Brainerd Rainy Rainy and superior Superior Undifferentiated EXPLANATION Model domain Figure 5.  Lithology and glacial geology of unconsolidated deposits in the Partridge River Basin, northeastern Minnesota.

Geologic Setting    9 47°30' 47°36' 47°42' 92°12' 92°6' 92°0' 91°54' 91˚48' Hillshade base from Minnesota Department of Natural Resources, 2011 Partridge River Basin from U.S. Geological Survey, 2012 Peat from Jirsa, 2016b 6 KILOMETERS 6 MILES Partridge River Basin Model domain EXPLANATION Peat deposits Figure 6.  Peat deposits in the Partridge River Basin, northeastern Minnesota.

10    Groundwater/Surface-Water Interactions in the Partridge River Basin and Hypothetical Future Mine Pits 92˚12' 92˚6' 92˚0' 91˚54' 91˚48' 47˚30' 47˚36' 47˚42' Bedrock units modified from Jirsa and others, 2011 Hillshade base from Minnesota Department of Natural Resources, 2011 6 MILES 6 KILOMETERS Lake Vermilion Formation Giants Range Granite Lower Ely Greenstone Upper Ely Greenstone Volcanic hornfels None provided Diabase to Gabbro intrusion Duluth Complex felsic rocks Duluth Complex anorthositic rocks Bald Eagle intrusion Partridge River intrusion Virginia Formation Thompson Formation Model domain Pokegama Quartzite Biwabik Iron-Formation Western Margin intrusion Greenword Lake intrusion South Kawishiwi intrusion EXPLANATION Animikie Group bedrock units Duluth Complex of Keweenawan Supergroup bedrock units Other bedrock units Figure 7.  Detailed bedrock geology in the Partridge River Basin, northeastern Minnesota.

Geologic Setting    11 A' Bedrock units simplified from Jirsa and others (2011) Hillshade base from Minnesota Department of Natural Resources, 2011 47°30' 47°36' 47°42' 92°12' 92°6' 92°0' 91°54' 91°48' A 6 MILES 6 KILOMETERS Bedrock units representated in the groundwater flow model EXPLANATION Virginia Formation Duluth Complex of Keweenawan Supergroup Biwabik Iron-Formation Archean units Model domain Cross-section trace Figure 8.  Simplified bedrock geology in the Partridge River Basin, northeastern Minnesota with cross-section trace shown in figure 9.

12    Groundwater/Surface-Water Interactions in the Partridge River Basin and Hypothetical Future Mine Pits DATUM IS NORTH AMERICAN VERTICAL DATUM OF 1988 Elevation, in feet above North American Vertical Datum of 1988 2,400 2,200 2,000 1,800 1,600 1,400 1,200 1,000 VERTICAL EXAGGERATION: 20x NO VERTICAL EXAGGERATION A - NORTHWEST A' - SOUTHEAST Mesabi Range Partridge River St. Louis River FEET Distance in miles EXPLANATION Biwabik Iron Formation of Animikie Group Virginia Formation Duluth Complex of Keweenawan Supergroup Archean units Glacial material Bedrock fractures Figure 9.  Cross section of the Quaternary unconsolidated material and bedrock formations from north of the Mesabi Range to the St. Louis River in the south. Upper fractured bedrock aquifer is shown as black hatching in the upper 200 feet of bedrock. The lower fractured bedrock aquifer is represented with less dense black hatching. A cross section with no vertical exaggeration is also shown to demonstrate that the thickness of the unconsolidated aquifer is negligible compared to the horizontal extent of the aquifer. The location of the cross-section trace is shown on figure 8.

Hydrogeologic Setting and Conceptual Model of the Flow System    13 2002). The intrusive bedrock units to the south and east in the model domain formed during the Midcontinent Rift around ~1.1 Ga (Jirsa and others, 2005) and are part of the Duluth Complex (Miller and others, 2001; Miller and others, 2002). The Duluth Complex is the target bedrock unit for proposed nonferrous mining in the region. For the purpose of modeling groundwater flow, the bedrock was simplified into the following four major bed­ rock units: (1) the Archean units to the north, (2) the Biwabik Iron-Formation, (3) the Virginia Formation, and (4) the Duluth Complex (plan view in fig. 8; cross section in fig. 9). Hydrogeologic Setting and Conceptual Model of the Flow System The groundwater system in the PRB interacts with streams, lakes, flooded mine pits, and wetlands. Water enters the basin as rainfall or melted snow, which either routes through surface water features, evaporates, is utilized by plants, or infiltrates through the land surface to the water table. Exchanges between the groundwater and surface-water systems occur in streams, wetlands, lakes, and mine pits. Pumping wells and mine dewatering in active mine pits rep­ resent groundwater withdrawals. Pumping wells are generally screened in the glacial aquifer or the Biwabik Iron-Formation (Lindholm and others, 1979). Regionally, groundwater flow is predominantly focused in the unconsolidated deposits and in fractures in the bed­ rock (Siegel and Ericson, 1980; Haserodt and others, 2019; T.K. Cowdery (USGS, written commun., 2020)). The concep­ tual model of the study area came from existing geologic maps and studies about the unconsolidated deposits (Minnesota Geological Survey, 2013), bedrock lithology (Jirsa and others, 2011), and other hydrologic studies in the region (Cotter and others, 1965; Oakes, 1970; Stark, 1977; Olcott and others, 1978; Lindholm and others, 1979; Siegel and Ericson, 1980; Jones, 2002; Adams and others, 2004; Myers, 2014; Tetra Tech, 2014; Barr, 2014a, 2014b; Haserodt and others, 2019; T.K. Cowdery (USGS, written commun., 2020)). The conceptual model of the PBR included the following hydrostratigraphic units: a thin, sometimes absent, unconsolidated aquifer in sand and gravel outwash or till that can have peat deposits occupying low-lying areas; an upper fractured bedrock aquifer in bedrock forma­ tions where fractures are assumed to be densest near land surface; and a lower bedrock aquifer that has lower permeability and fewer fractures than the upper fractured bedrock aquifer and where permeability decreases with depth as fracture density and fracture connectivity decrease. Hydrogeologic properties in the unconsolidated aquifer are assumed to vary by glacial unit and lithology (fig. 5) and with the presence of peat deposits (fig. 6). Barr (2014a) col­ lected data on unconsolidated materials at a proposed mine site in the PRB and determined the Rainy till, the major glacial unit in the PRB, to be a highly heterogenous unit with materi­ als that ranged in particle size from dense clay to sand lenses. In addition to varying with depth, the hydrogeologic proper­ ties of the bedrock were characterized as varying spatially by major bedrock unit (fig. 8). Local heterogeneities from faults and fractures are likely within these simplified bedrock units but are not well characterized at the basin scale. Conceptually, the bedrock was assumed to behave as an equivalent porous medium consisting of zones, following the approach used in other groundwater-flow models in the region (Haserodt and others, 2019; T.K. Cowdery (USGS, written commun., 2020); Barr, 2014b). Table 1 provides a summary of hydraulic conductiv­ ity and vertical anisotropy values for the glacial and bedrock units in the PRB from the following sources: Barr (2014a); the equivalent hydraulic conductivity values in the Haserodt and others (2019) GFLOW model; the calibrated Iron Range MODFLOW-NWT model (T.K. Cowdery (USGS, writ­ ten commun., 2020); 2021b); the Myers (2014) calibrated MODFLOW model; and the Barr (2014b) calibrated site-scale MODFLOW model in the PRB. In general, reported horizon­ tal hydraulic conductivity of unconsolidated materials ranged from 10-5 to 102 feet per day (ft/day), depending on the unit and source. The ratio of the horizonal to vertical hydraulic conductivity (anisotropy) of the unconsolidated materials ranged from 1.1 to 80,000. The highest reported bedrock horizontal hydraulic conductivity values were in the Biwabik Iron-Formation and the lowest were in the Virginia Formation and Archean units (table 1). Bedrock vertical hydraulic con­ ductivity was set at one-tenth the horizontal hydraulic con­ ductivity for the bedrock in the Barr (2014b) model, and the bedrock anisotropy ranged from 0.008 to 1,400 for the other models. Generally, vertical hydraulic conductivity is assumed to be less than horizontal due to depositional layering, but the geometry of fracture networks in crystalline bedrock may not adhere to this. The extent of connectivity in the PRB between the groundwater in the bedrock and the groundwater in the uncon­ solidated aquifer depends on the extent of the bedrock fracture network and if it is connected to the shallow flow system. Adams and others (2004) used observations of mine pit wall seeps to infer that most mine pit seepage occurs at the inter­ face between the unconsolidated material and the top of bed­ rock, with only minor seepage from fractures in the Biwabik Iron-Formation. Furthermore, observations of mine pit waterlevel elevations over time showed a step-wise filling rate in pits with minimal unconsolidated material along the rim and smoother filling rates when unconsolidated materials were thicker; this was interpreted as precipitation dominated filling for the step-wise pattern and greater groundwater inflow from unconsolidated materials for the smoother pattern (Adam and

14    Groundwater/Surface-Water Interactions in the Partridge River Basin and Hypothetical Future Mine Pits Table 1.  Summary of regional hydraulic conductivity and vertical anisotropy estimates for glacial and bedrock units. Headnote: -, no data : greater than; less than Estimated horizontal hydraulic conductivity, in feet per day Horizontal:vertical hydraulic conductivity, unitless Unit Source 2.1E-5–6.7 Rainy Lobe till Stark (1977); Siegel and Ericson (1980). 0.01–31 Rainy Lobe till Barr (2014b); aquifer testing at mine site. 0.8–28 Rainy Lobe till T.K. Cowdery (USGS, written commun., 2020) calibrated model values. 0.004–362 Sand and gravel near Iron Range Jones (2002), Siegel and Ericson (1980), Stark 1977. 0.07–2.4 2.6–225 Sand and gravel T.K. Cowdery (USGS, written commun., 2020) calibrated model values. Till Haserodt and others (2019) average calibrat­ ed GFLOW value for the till in the PRB. 0.06–167 20–59,643 Upland deposits - all nonwetland material Barr (2014b) calibrated model range. 3 to 24 20–52 Surficial aquifer Myers (2014). 0.001–0.1 Peat Siegel and Ericson (1980). Peat T.K. Cowdery (USGS, written commun., 2020) calibrated model value. 0.003–224 1.1–80,000 Wetland deposits Barr (2014b). Bedrock Units Estimated horizontal hydraulic conductivity, in feet per day Horizontal:vertical hydraulic conductivity, unitless Depth, in feet Source Duluth Complex of Keweenawan Supergroup 3E-4– 0.4 Literature summarized in Barr (2014a). 0.3–4.1 Upper bedrock layer T.K. Cowdery (USGS, written commun., 2020) calibrated model value. 4.0E-04 Lowest bedrock layer T.K. Cowdery (USGS, written commun., 2020) calibrated model value. 5E-4 to 9.5 0.008–130 Myers (2014). 3E-4–4E-2 Barr (2014b) aquifer test values. 4.4E-04 Barr (2014b) calibrated model value. Virginia Formation 2E-3–0.7 Literature summarized in Barr (2014a). 0.06–1.4 Upper bedrock layer T.K. Cowdery (USGS, written commun., 2020) calibrated model values. 3.0E-04 Lowest bedrock layer T.K. Cowdery (USGS, written commun., 2020) calibrated model values. 2E-3–6.8E-1 Barr (2014b) aquifer test values. 3.1E-01 Barr (2014b) calibrated model value. Barr (2014b) calibrated model value. Archean units 4E-3–2 Literature summarized in Barr (2014a). 0.9–20 Upper bedrock layer T.K. Cowdery (USGS, written commun., 2020) calibrated model values. 4.E-03 Lowest bedrock layer T.K. Cowdery (USGS, written commun., 2020) calibrated model values. 0.005–2.6 1–1.6 Myers (2014).

Water Use    15 others, 2004). Siegel and Ericson (1980) also noted that bed­ rock in the region was likely relatively impermeable with per­ meability only from secondary features like joints, fractures, and leached zones. Fractures and joints in bedrock are thought to be most extensive in the upper 200–300 ft of rock (Siegel and Ericson, 1980). Barr (2014b) reported that estimated aquifer properties from pumping tests in the upper 300 ft of the Virginia Formation were 3–5 times as conductive as those from depths of 300–600 ft. Regional well pumping yields in the upper bedrock also provide evidence of fracture connectiv­ ity. Water supply wells in the upper crystalline bedrock can produce yields of 5–15 gallons per minute in the Giants Range Granite and Duluth Complex, as much as 30 gallons per min­ ute in the Virginia Formation, and as much as 1,000 gallons per minute in the Biwabik Iron-Formation (Siegel and Ericson, 1980). This pumping information indicates that bedrock groundwater flow can occur in the upper, more fractured rock. Although the location, extent, and magnitude of connectivity between the bedrock and the unconsolidated aquifer systems is not currently well characterized, it is important for understand­ ing the local groundwater-flow systems in the PRB and would benefit from future data collection. Water Use Water use in the PRB is largely from surface-water fea­ tures. The drinking water supply for Hoyt Lakes is Colby Lake (Hoyt Lakes, 2016; Adams and others, 2004) and the drinking water supply for Aurora is the St. James Pit, a flooded, openpit, iron ore mine that closed in the 1960s (MDNR, 2019b; Adams and others, 2004). In 2011, there were no permitted groundwater withdrawals within the basin (permits apply to wells pumping more than 10,000 gallons per day or 1 million gallons per year; MDNR, 2019c), but there were six permitted groundwater withdrawals (MDNR, 2019d) to the northwest of the basin, within the model domain. All permitted withdrawals were from glacial aquifers. Pumping rates for these permit­ ted groundwater withdrawals were modeled using the 2011 pumping data (MDNR, 2019d) and ranged from 3–97 gallons per minute. In addition to pumping for water supply, there are appropriated mine pit withdrawals for mine dewatering or mine pit level control that are then discharged into other pits or streams (Avery Guertin, Minnesota Department of Natural Resources, written commun., 2019). There are five permitted mine pit withdrawals within the study area: one in the PRB, three in the Dunka River Basin, and one permitted withdrawal in the western part of the PRB outside the model domain that is discharged to a permit point within the basin. Table 1.  Summary of regional hydraulic conductivity and vertical anisotropy estimates for glacial and bedrock units.—Continued Headnote: -, no data : greater than; less than Bedrock Units Estimated horizontal hydraulic conductivity, in feet per day Horizontal:Vertical Hydraulic Conductivity, unitless Depth, in feet Source 0–2 Barr (2014b) aquifer test values. 3.E-02 Barr (2014b) calibrated model value. Biwabik Iron-Formation of Animikie Group 0.2–16 Literature summarized in Barr (2014a). 2.6–25 Upper bedrock layer T.K. Cowdery (USGS, written commun., 2020) calibrated model values. 0.04–17 Lowest bedrock layer T.K. Cowdery (USGS, written commun., 2020) calibrated model values. 0.5 to 85 38–1400 Myers (2014). 9.E-01 Barr (2014b) calibrated model value. 3–4 Haserodt and others (2019) average cali­ brated GFLOW model value.

16    Groundwater/Surface-Water Interactions in the Partridge River Basin and Hypothetical Future Mine Pits Groundwater Flow Model Construction A finite-difference groundwater-flow model of the PRB was constructed using MODFLOW-NWT with the Newton Raphson solver (MODFLOW–NWT: Niswonger and others, 2011). This section describes the overall model construction and model inputs and outputs. Appendix 1 provides additional details on the construction of packages that required substan­ tial data processing to implement. Discretization The PRB model was discretized as an inset model of a MODFLOW-NWT model (T.K. Cowdery (USGS, written commun., 2020)) of the Mesabi Iron Range. The PRB model grid has 550 rows and 478 columns of 200 ft by 200 ft grid cells that cover the model domain area shown in figure 1.The grid resolution was chosen to be fine enough to assess ground­ water/surface-water exchanges in the PRB and represent mine pit geometries while balancing the need for a fast model run time to allow for uncertainty analysis and calibration. Vertically, the model domain was discretized into 11 layers, 4 in the unconsolidated material and 7 in the bedrock. The top of layer 1 was defined as the land surface and the bottom of layer 11 was set to a depth of 1,000 ft below the bedrock surface, which was selected to accommodate current and future mine pit depths. The geometry of the layering was controlled by the following surfaces: model top elevation (top of layer 1) is the land-surface elevation determined from light detecting and ranging (lidar) data (MDNR, 2011) and averaged across each model cell; and top of bedrock elevation (top of layer 5) is the bedrock topography surfaces (combined central and south­ eastern arrowhead regions from Dengler and Wagner, 2017; Dengler and others, 2017), resampled to an average bedrock elevation in each model cell. The thickness of the unconsolidated materials (layers 1–4) was the difference between the lidar-derived model top elevation (MDNR, 2011) and the elevation of the top of the bedrock surface (Dengler and Wagner, 2017; Dengler and others, 2017). In instances where the mapped bedrock eleva­ tion exceeded the lidar elevation, the bedrock was reset below the model top elevation. Layer top and bottom elevations were decided using the same logic used for the T.K. Cowdery (USGS, written commun., 2020) MODFLOW model and described further in appendix 1. The modeled thickness of the unconsolidated mate­ rial is shown in figure 10 for layers 1–4. A large part of the PRB has bedrock that crops out at land surface resulting in thin unconsolidated material even in layer 1. Layers 2 and 3 (figs. 10B and 10C) are absent across much of the PRB and only present in two bedrock valleys, one to the north and the other to the south of the PRB. Layer 4 is not present in most of the model domain. Figure 11 shows a cross section of model layers for a trace from A–A’ (trace shown on fig. 8). The glacial material is very thin across much of the cross section, except for thicker deposits in bedrock valleys to the north and south of the PRB. Most of the model thickness is from the bedrock. Boundary Conditions Boundary conditions are used to represent sources and sinks of water (fig. 12). The PRB model includes the following boundary conditions: a spatially variable net infiltration flux from a Soil Water Balance (SWB) model was applied across the top of the model; the areal net infiltration was locally modified by rejection and routing in areas where the water table is at land surface (simulated by MOD­ FLOW’s Unsaturated Zone Flow [UZF] package; see Hunt and others [2008] and Feinstein and others [2019] for discussions of the rejection and routing); specified flux boundary conditions simulated by MOD­ FLOW’s Well (WEL) package along the model perim­ eter; the flux rates for each perimeter cell were calcu­ lated using the regional Iron Range MODFLOW model of T.K. Cowdery (USGS, written commun., 2020), lake head-dependent flux boundary conditions were represented with MODFLOW’s River (RIV) package; stream head-dependent flux boundary conditions were represented with MODFLOW’s Streamflow Routing (SFR2) package; water use pumping represented with MODFLOW’s Multi-node Well (MNW2) package; fully pumped mine pits represented with MODFLOW’s Drain (DRN) package; and water transfers from pits represented with MOD­ FLOW’s well (WEL) package. Net Infiltration Recharge is water that moves through the ground and crosses the water table to replenish the groundwa­ ter system. Estimates of net infiltration, or water that goes past the root zone, came from a published SWB model for the St. Louis River Basin (SLRB; Smith, 2017; T.K. Cowdery (USGS, written commun., 2020)) and the Minnesota statewide SWB model for areas outside the SLRB (Smith and Westenbroek, 2015; T.K. Cowdery (USGS, written commun., 2020)). SWB uses

Groundwater Flow Model Construction    17 Hillshade base from Minnesota Department of Natural Resources, 2011 92°12' 92°6' 92°0' 91°54' 91°48' 47°42' 47°36' 47°30' Layer 1 thickness, in feet 1 as shown as underlaying hillshade 10 to 15 15 to 20 6 MILES 6 KILOMETERS EXPLANATION Model domain A Figure 10.  A, Layer 1 thickness and mapped wetlands, B, Layer 2 thickness, C, Layer 3 thickness, and D, Layer 4 thickness. Areas in gray show the underlying hillshade and represent pinched out unconsolidated materials.

18    Groundwater/Surface-Water Interactions in the Partridge River Basin and Hypothetical Future Mine Pits Hillshade base from Minnesota Department of Natural Resources, 2011 92°12' 92°6' 92°0' 91°54' 91°48' 47°42' 47°36' 47°30' 6 MILES 6 KILOMETERS Layer 2 thickness, in feet 1 as shown as underlaying hillshade 10 to 20 20 to 40 40 to 60 60 to 80 EXPLANATION Model domain B Figure 10.—  A, Layer 1 thickness and mapped wetlands, B, Layer 2 thickness, C, Layer 3 thickness, and D, Layer 4 thickness. Areas in gray show the underlying hillshade and represent pinched out unconsolidated materials.

Groundwater Flow Model Construction    19 Hillshade base from Minnesota Department of Natural Resources, 2011 92°12' 92°6' 92°0' 91°54' 91°48' 47°42' 47°36' 47°30' 6 MILES 6 KILOMETERS Layer 3 thickness, in feet 1 as shown as underlaying hillshade 20 to 40 40 to 60 60 to 80 80 to 100 EXPLANATION Model domain Figure 10.—  A, Layer 1 thickness and mapped wetlands, B, Layer 2 thickness, C, Layer 3 thickness, and D, Layer 4 thickness. Areas in gray show the underlying hillshade and represent pinched out unconsolidated materials.

20    Groundwater/Surface-Water Interactions in the Partridge River Basin and Hypothetical Future Mine Pits Hillshade base from Minnesota Department of Natural Resources, 2011 92°12' 92°6' 92°0' 91°54' 91°48' 47°42' 47°36' 47°30' 6 MILES 6 KILOMETERS Layer 4 thickness, in feet 1 as shown as underlaying hillshade 20 to 40 40 to 60 60 to 80 80 to 105 EXPLANATION Model domain D Figure 10.—  A, Layer 1 thickness and mapped wetlands, B, Layer 2 thickness, C, Layer 3 thickness, and D, Layer 4 thickness. Areas in gray show the underlying hillshade and represent pinched out unconsolidated materials.

Groundwater Flow Model Construction    21 daily climate data and readily available geographic informa­ tion system datasets of landscape and soil properties to cal­ culate a Thornthwaite-Mather SWB on a daily time step and produce an estimate of infiltration. SWB output files were pro­ cessed into data that represents 2011–2013 for the 2011–2013 mining conditions model and 1995–2015 for the 1995–2015 average conditions model. The data processing steps used to incorporate the SWB data into the model are further described in appendix 1. The spatial average SWB net infiltration rate for the model domain was 6.4 inches per year (in/yr) for 1995–2015 and 5.7 in/yr for 2011–2013. The average 2011–2013 SWB net infiltration values are shown in figure 13. The SWB net infiltration for the PRB is slightly lower than calibrated model recharge averages from Myers (2014; 7.2 in/year) and Tetra Tech (2014; 8.2 in/year), but these models simulated a differ­ ent period than the study models. Changes to Infiltration Applied by the UZF Package For area with a water table not near land surface, all net infiltration is assumed to fully reach the water table, though the timing of recharge can be affected by unsaturated zone processes (Hunt and others, 2008). However, unsaturated zone processes can also be important for steady-state models because overland flow can be appreciable when the water table is at or near land surface (see Hunt and others, 2008) because the subsurface lacks storage to accommodate infiltration applied. A large part of the PRB model domain is wetlands— areas were the water table is at or near land surface. Typically, high water tables are also present in low lying areas, near lakes, and in areas that have a stream drainage that is not explicitly represented by the streamflow routing network (Leaf and others, 2015) in the SFR2 package. With its focus on the root zone, SWB cannot account for rejection of infiltration when the water table is at the land surface. This can result in simulations of unrealistically high recharge as net infiltration estimated by SWB is added to areas with high water tables and no subsurface storage. MODFLOW groundwater models that use the MODFLOW recharge (RCH) package do not allow for groundwater seepage in areas where the groundwater elevation exceeds the land-surface elevation, a process that would reduce the extent of “flooding” in a groundwater model. Flooded cells in a groundwater flow model can produce inac­ curate results from particle tracking; result in undersimulation of large stresses (such as pit pumping) by offsetting declines with inaccurately high groundwater elevations; and could make topographic lows that should be groundwater discharge zones into local groundwater recharge areas (Leaf and others, 2015). The UZF package can be used to alleviate groundwater flooding and the problems associated with flooded model cells. The UZF package (Niswonger and others, 2006) was used to represent the processes occurring in areas with high water tables. UZF allows for rejection of recharge and groundwater seepage (collectively considered groundwater runoff) to the land surface when the water table is near the land surface; that is, UZF partitions net infiltration estimated by SWB into components of groundwater runoff and ground­ water recharge. The amount of simulated groundwater runoff is based on the hydraulic head relative to land-surface eleva­ tion (as specified in the model top elevation array). Instead of a sharp cutoff from infiltration to rejection and ground­ water seepage, UZF gradually applies the transition using a SURFDEP variable that represents natural undulations in the land surface. Groundwater runoff is zero when the simulated water table elevation is below the model top elevation minus one-half SURFDEP; groundwater runoff increases to 100 per­ cent of the SWB net infiltration when the water table elevation is at or above the model top elevation plus one-half SURFDEP. The relation of SURFDEP, hydraulic head, land-surface eleva­ tion, and groundwater runoff is discussed in detail by Feinstein and others (2019). For the PRB model, a SURFDEP of 2.5 ft was used across the model and is the median difference between the minimum and maximum lidar elevations in model cells with mapped peat deposits. Peat deposits should represent flatter, low-lying part of the model where the difference between the minimum and maximum land-surface elevations describe small surface undulations. Areas with mapped peat are also likely to have groundwater runoff because of a high water table and landscape position. UZF is designated as active or inactive using the IUZFBND. The rejected groundwater recharge and the groundwater discharge from active areas were routed to nearby stream cells using routing information supplied to UZF in the IRUNBND. Groundwater runoff that is not routed is assumed to evapotranspire to the atmosphere. Appendix 1 describes how the IRUNBND and IUZFBND were constructed for the PRB model. Regional Groundwater Flow Regional groundwater flow entering or leaving the PRB model domain perimeter was represented using the well pack­ age (WEL) with a specified flux boundary condition along the model edges (fig. 12). Layers 1–7 were assigned a flux based on groundwater-flow rates from the Iron Range model (T.K. Cowdery (USGS, written commun., 2020)). Layers 9–11 extended below the Iron Range model so the layer 8 flux from the Iron Range model was divided among layers 8–11 in the PRB model using the relative layer transmissivi­ ties. Given the relatively lower hydraulic conductivity of the deeper bedrock, it was assumed that little groundwater moved regionally through the deeper system and the assigned fluxes at this depth likely have minimal effect on the model. To reduce overall boundary effects on model forecasts, the model domain includes important river boundary conditions outside the PRB to account for the effects of competing hydrologic sinks (St. Louis River to the south, Dunka River to the east, and Embarrass River to the north; fig. 2).

22    Groundwater/Surface-Water Interactions in the Partridge River Basin and Hypothetical Future Mine Pits 2,250 2,000 1,750 1,500 1,250 1,000 VERTICAL EXAGGERATION: 20X A - NORTHWEST A' - SOUTHEAST Distance, in miles FEET Datum is North American Vertical Datum of 1988 Mesabi Range Partridge River St. Louis River EXPLANATION Layer 11 Layer 10 Layer 9 Layer 8 Layer 7 Layer 6 Layer 5 Layer 3 Layer 2 Layer 1 Biwabik Iron-Formation of Animikie Group Virginia Formation Duluth Complex of Keweenawan Supergroup Archean units Glacial material [Layer 4 is absent in this cross section] Figure 11.  Cross section showing model layers including layers 1–3 in the glacial material and 5-11 in the bedrock. Layer 4 is absent in this cross section. The cross-section trace is shown in figure 8.

Groundwater Flow Model Construction    23 Hillshade base from Minnesota Department of Natural Resources, 2011 Partridge River Basin Outline from U.S. Geological Survey, 2012 NHDPlus stream lines from McKay and others, 2012 47°30' 47°36' 47°42' 92°12' 92°6' 92°0' 91°54' 91°48' 6 KILOMETERS 6 MILES EXPLANATION Pit transfers represented with the WEL package Pumped pits represented with the DRN package Lakes represented with the RIV package Flooded mine pits represented with high hydraulic conductivity zones Partridge River Basin Model domain Edge boundary fluxes represented with the WEL Package Streams in National Hydrological Data Plus Streams represented using the Streamflow Routing (SFR2) package Pumping wells with permitted groundwater withdrawls Pit overflow location Figure 12.  Boundary conditions presenting streams, lakes, groundwater fluxes along the model edges, pumped pits, pit transfers, and groundwater pumping wells with water use permits.

24    Groundwater/Surface-Water Interactions in the Partridge River Basin and Hypothetical Future Mine Pits 47°30' 47°36' 47°42' 92°12' 92°6' 92°0' 91°54' 91°48' Hillshade base from Minnesota Department of Natural Resources, 2011 6 KILOMETERS 6 MILES EXPLANATION Average 2011 to 2013 net infiltration, in inches per year 10 to 11 11 to 12 12 to 13 13 to 14 14 Model domain Figure 13.  SWB-estimated net infiltration for 2011–2013 mining conditions model of the Partridge River Basin, northeastern Minnesota. The 1995–2015 average conditions model has the same pattern of infiltration but with slightly different values.

Groundwater Flow Model Construction    25 Streams Streams were modeled as a head-dependent flux bound­ ary condition with the Streamflow Routing (SFR2) package (Niswonger and Prudic, 2005). SFR2 allows for the interac­ tion between groundwater and streams and routes streamflow along a stream system. In this work, streamflow represents base-flow conditions, where flows in the streams reflect exchanges with the groundwater system and do not include stormflow components. The interaction between the stream and the aquifer is determined by the stage in the stream relative to the groundwater elevation (hydraulic head) in the aquifer and a conductance term based on the streambed dimensions and the vertical hydraulic conductivity of the streambed sediment. Losses from the stream to the aquifer are limited to the available streamflow in that reach. Stream stage is estimated using Manning’s Equation (Niswonger and Prudic, 2005). Information in NHDPlus version 2 (McKay and others, 2012) was used to create the SFR2 model files with the python package SFRmaker (Leaf, 2018). A few headwater stream segments were removed or edited from the NHDPlus network if they crossed mining features (like mine pits that are currently pumped dry) or if they were disconnected from the stream network. NHDPlus streamlines outside the model domain and streams represented with the SFR2 package can be seen in figure 12. Mining discharge to the stream system was simulated as specified inflows to SFR2 using reported permit data, as discussed in appendix 1. Lakes Lakes were modeled as a head-dependent flux bound­ ary condition using the River (RIV) package (Harbaugh and others, 2000) with the river stage set to the lake-surface elevation. The input datasets and formulation of the RIV package is further discussed in appendix 1. Wells In 2011, there were 6 permitted groundwater withdraw­ als (MDNR, 2019d) within the model domain (permits apply to wells pumping more than 10,000 gallons per day or 1 mil­ lion gallons per year [MDNR, 2019c; fig. 12]). Groundwater withdrawals were modeled using the MODFLOW MNW2 package to allow for wells penetrating multiple model layers (Konikow and others, 2009). Well characteristics came from the parent model (T.K. Cowdery (USGS, written commun., 2020); Cowdery and others, 2021b) and included a well radius of 0.5 ft, a skin radius of 1 ft, and a skin hydraulic con­ ductivity of 10 ft/day. Wells were all screened in the glacial aquifer in layers 2–4. Mine Pits Pumped mine pits were modeled using drain (DRN) cells (Harbaugh and others, 2000) and flooded pits were modeled as high hydraulic conductivity zones (Anderson and others, 2002; Hunt and others, 2003). DRN cells act strictly as a ground­ water sink and remove water from the groundwater model. Simulating flooded pits with high hydraulic conductivity zones allows for pit water levels to be simulated rather than speci­ fied, which is a more robust test the model. The distribution of pumped and flooded pits is shown in figure 12. Appendix 1 provides additional discussion on how pits were represented in the PRB model. The water levels in flooded pits and pumped water from dry pits were assumed to include groundwater contributions to the pits and surface runoff from precipitation events. Runoff to a pit can be approximated using annual precipita­ tion, a drainage basin area, and an assumed runoff coefficient. This volume was added to the pit using an injection WEL package well in the high hydraulic conductivity cells of the pit. Appendix 1 provides information on how the surface catchment of the pit was estimated. An initial runoff coef­ ficient of 0.4 (40 percent) was used based on the Barr (2008) runoff estimates of 0.3 (30 percent) for undisturbed areas and 0.6 (60 percent) in mining features for drainage basins near the Iron Range. The runoff coefficient was adjusted for each mine pit during calibration. Reported information on the pumped and flooded pits was used as calibration data to help constrain the bedrock properties. For pumped pits with reported pumping rates, the 2011 pumping rates (Avery Guertin, Minnesota Department of Natural Resources, written commun., 2019) were used as calibration targets. Because flooded pit elevations are not specified explicitly as a boundary condition (like the lakes or pumped pits), the flooded pit elevations from the 2011 lidar data were used as a calibration target. Hydraulic Conductivity Vertical and horizontal hydraulic conductivity were rep­ resented using zones in the upper unconsolidated layers and lower bedrock layers (layers 1–4, and 7–11, respectively) and zones with pilot points in the two upper bedrock layers (layers 5 and 6). Hydraulic conductivity values were specified using MODFLOW’s Upstream Weighting (UPW; Niswonger and others, 2011) package. Unconsolidated and Surficial Mining Features The unconsolidated material was divided into hydrau­ lic conductivity zones based on glacial units and mapped lithology (fig. 5) and on the location of waste rock, over­ burden stockpiles, tailing basins, and flooded and pumped pits. Hydraulic conductivity zones for layer 1 are shown in figure 14. The major units in layer 1 include the Rainy

26    Groundwater/Surface-Water Interactions in the Partridge River Basin and Hypothetical Future Mine Pits 47°30' 47°36' 47°42' 92°12' 92°6' 92°0' 91°54' 91°48' Hillshade base from Minnesota Department of Natural Resources, 2011 Partridge River Outline from U.S. Geological Survey, 2012 6 KILOMETERS 6 MILES EXPLANATION Layer 1 hydraulic conductivity zones Zone 1: Rainy Lobes diamicton Zone 2: Rainy Lobe or non-classified sand Zone 3: Rainy Lobe gravel Zone 4: High hydraulic conductivity over pumped pits Zone 5: High hydraulic conductivity over flooded pits Zone 6: Stockpiles Zone 7: Disturbed by mining Zone 8: Tailing basin Zone 9: Rainy lobe silt and clay Zone 11: Peat Zone 40: Brainerd Sublobe sand Zone 41: Brainerd Sublobe gravel Zone 42: Brainerd Sublobe Diamicton Zone 43: Brainerd Sublobe silt and clay Zone 44: Rainy and Superior lobe sand/ gravel Zone 45: Rainy and Superior Lobe diamicton Zone 46: Superior Lobe sand Zone 47: Superior Lobe gravel Partridge River Basin Model Domain Figure 14.  Layer 1 hydraulic conductivity zones and represented geologic units.

Groundwater Flow Model Construction    27 Hillshade base from Minnesota Department of Natural Resources, 2011 Partridge River Outline from U.S. Geological Survey, 2012 47°30' 47°36' 47°42' 92°12' 92°6' 92°0' 91°54' 91°48' 6 KILOMETERS 6 MILES EXPLANATION Layer 4 hydraulic conductivity zones Zone 1: Rainy Lobe diamicton Zone 2: Rainy Lobe or nonclassified sand Zone 4: High hydraulic conductivity over pumped pits Zone 5: High hydraulic conductivity over flooded pits Zone 7: Disturbed by mining Zone 40: Brainerd Subobe sand Zone 42: Brainerd Subobe diamicton Zone 44: Rainy and Superior Lobe silt and clay Zone 45: Rainy and Superior Lobe diamicton Zone 46: Superior Lobe gravel Partridge River Basin Model domain Figure 15.  Layer 4 hydraulic conductivity zones and represented geologic units.

28    Groundwater/Surface-Water Interactions in the Partridge River Basin and Hypothetical Future Mine Pits 47°30' 47°36' 47°42' 92°12' 92°6' 92°0' 91°54' 91˚48' Hillshade base from Minnesota Department of Natural Resources, 2011 Bedrock units simplified from Jirsa and others (2011). 6 MILES 6 KILOMETERS Bedrock units represented in the groundwater flow model Duluth Complex of Keweenawan Supergroup Virginia Formation Biwabik Iron-Formation of Animikie Group Archean units Model domain Layer 5 pilot points EXPLANATION A Figure 16.  Location of horizontal and vertical hydraulic conductivity pilot points for the upper bedrock layers.

Groundwater Flow Model Construction    29 47°30' 47°36' 47°42' 92°12' 92°6' 92°0' 91°54' 91˚48' Hillshade base from Minnesota Department of Natural Resources, 2011 Bedrock units simplified from Jirsa and others (2011). 6 MILES 6 KILOMETERS Bedrock units represented in the groundwater flow model Duluth Complex of Keweenawan Supergroup Virginia Formation Biwabik Iron-Formation of Animikie Group Archean units Model domain Layer 6 pilot points EXPLANATION B Figure 16.  Location of horizontal and vertical hydraulic conductivity pilot points for the upper bedrock layers.

30    Groundwater/Surface-Water Interactions in the Partridge River Basin and Hypothetical Future Mine Pits Lobe diamicton, sand, peat, and material mapped as dis­ turbed by mining activities. In addition to pits, other mining features represented by hydraulic conductivity zones are stockpiles and the fine tailings in a tailings basin on the PRB northern boundary. Layers 2–4 were assigned nearly identical hydraulic conductivity zones, with the exception of the tailings basin fines present in layers 1–3 but not layer 4. The zones in layers 1–4 were from the Iron Range model (T.K. Cowdery (USGS, written commun., 2020); 2021b). Layers 2–4 have many of the glacial units in layer 1 in addition to buried sands lenses. Peat was assumed to exist in layer 1 where the maximum layer thickness was 20 ft. Figure 15 shows the hydraulic conductiv­ ity zones for model layer 4. Bedrock Bedrock in the model domain consists of crystalline, dipping Proterozoic bedrock units. Jirsa and others (2011) mapped bedrock units at the top of the bedrock surface and these were simplified into four units for the model (fig. 8). Each of the four units was represented by its own hydraulic conductivity zone in every bedrock layer. These hydraulic conductivity zones from the Jirsa and others (2011) bedrock surface were translated down to lower layers using published bedrock dip angles and the Minnesota County Well Index database (County Well Index, 2018), as described in appen­ dix 2. A cross section with the dipping model layers through the PRB is shown in figure 9 with the section trace shown on figure 8. A network of vertical and horizontal hydraulic conductiv­ ity pilot points (fig. 16) within the bedrock zones in layers 5 and 6 were used to provide parameter flexibility during cali­ bration for the upper bedrock units where fractures are most prevalent. Lower bedrock units used hydraulic conductivity zones because of the lack of calibration data in these units. Model Solver The groundwater flow model used the MODFLOW-NWT solver, version 1.1.4 (Niswonger and others, 2011). MODFLOW-NWT is a Newton-Raphson formulation of MODFLOW-2005 that can robustly handle problems involving drying and rewetting of model cells during solver iterations (see Hunt and Feinstein, 2012). The PRB model had several features that could slow model run times includ­ ing steep water-table boundaries near pumped pits, thin layers of unconsolidated material over crystalline bedrock, streamflow routing, and UZF to simulate rejected recharge and groundwater seepage. MODFLOW-NWT solver settings were adjusted to optimize model run times and convergence. Solver settings are listed in the MODFLOW-NWT solver package files (*.nwt) in the accompanying model data release (Haserodt and others, 2021). Model Calibration The parameter-estimation software suite PEST_HP (Doherty, 2020) was used to obtain a best fit of the modeled to measured groundwater levels, streamflow, pit inflow rates, and mapped wetland deposits. The model calibration included (1) an average conditions model with 1995–2015 calibration data, a period where more data were available; and (2) a current mining conditions model with 2011–2013 calibration data, representing a period when the most consistent set of mining and lidar data were available. In addition to separate calibra­ tion datasets, the models have different SWB infiltration dis­ tributions, representing potential infiltration for UZF, averaged for either 2011–2013 or 1995–2015. Targets As shown in table 2 and discussed in the following sec­ tions, the 1995–2015 average conditions model was cali­ brated using 288 hydraulic head targets, 9 streamflow flux targets, and 2 composite targets representing the distribution of wetlands based on mapped peat (fig. 6) relative to areas with simulated UZF seepage. The 2011–2013 mining condi­ tions model was calibrated using 67 hydraulic head targets, 9 streamflow flux targets, and 2 pit inflow targets. Targets are from a range of sources, and required processing, as discussed in the following sections. Hydraulic Heads Hydraulic head targets are observations of measured groundwater elevations from wells and surface water features where the water level is assumed to reflect the water table. Most hydraulic head targets were derived from calibration targets used by T.K. Cowdery (USGS, written commun., 2020). The hydraulic head targets were grouped based on quality designations assigned by T.K. Cowdery (USGS, writ­ ten commun., 2020). Additional groups were made for targets that were within 0.62 miles (1,000 meters) of a mining feature where groundwater elevations likely depend on the local buildout and pit water levels. Other hydraulic head target data sources include the U.S. Geological Survey (USGS) National Water Information Systems database (U.S. Geological Survey, 2019) and mine pit water levels, as measured in the 2011 lidar data (MDNR, 2011). Group designations and data sources are listed in table 2. Fluxes Flux targets included streamflow and pit inflow targets. Eight of the nine streamflow targets are from a 2018 USGS synoptic flow survey and were pulled from the National Water Information Systems database (U.S. Geological Survey, 2019). Streamflow target site numbers can be found in the ancillary

Model Calibration    31 Table 2.  Weighting used for the target observation groups. [MDNR, Minnesota Department of Natural Resources; MN, Minnesota; %; percent; NWIS, National Water Information Sysytem (U.S. Geological Survey, 2019); USGS, U.S. Geological Survey; day/ft3, day per foot cubed; ft-1, per foot; lidar, light detection and ranging; UZF, unsaturated zone flow] Group name Model Description Weight Data source Number of targets Flux targets flux1 1995–2015 Streamflow targets including eight measure­ ments from a 2018 synoptic flow survey and one MDNR streamgage record. One of the synoptic measurements on the Dunka River was zero-weighted because it is not fully within the model domain. Flows were compensated to average 1995–2015 condi­ tions using the streamflow record at USGS St. Louis near Skibo, MN streamgage. Varies; started with a weight based on the 95% confidence interval and then adjusted to balance the contribution to the objective function from hydraulic heads and fluxes and added weight to more downstream flux targets that aggregate a larger drainage basin area. MDNR (2018), NWIS data­ base flux2 2011–2013 Streamflow targets including eight measure­ ments from a 2018 synoptic flow survey and one MDNR streamgage record. One of the synoptic measurements on the Dunka River was zero-weighted because it is not fully within the model domain. Flows were compensated to average 2011–2013 condi­ tions using the streamflow record at USGS St. Louis near Skibo, MN streamgage. Varies; started with a weight based on the 95% confidence interval and then adjusted to balance the contribution to the objective function from hydraulic heads and fluxes and added weight to more downstream flux targets that aggregate a larger drainage basin area. MDNR (2018), NWIS data­ base flux 2011–2013 Measured pit inflow rates reported in the MDNR (2011) dataset. 9 day/ft3 A. Guertin, MDNR, written com­ mun., Head targets head15 2011–2013 Hydraulic head targets representing data collected between 2011 and 2013 used in T.K. Cowdery (USGS, written commun., 2020) with a best quality designation that are near the Mesabi Iron Range. 0.04 ft-1 T.K. Cowdery (USGS, written commun., 2020) head16 2011–2013 Hydraulic head targets representing data collected between 2011 and 2013 used in T.K. Cowdery (USGS, written commun., 2020) with a best quality designation that are not near the Mesabi Iron Range. 0.2 ft-1 T.K. Cowdery (USGS, written commun., 2020) head17 2011–2013 Hydraulic head targets representing data collected between 2011 and 2013 used in T.K. Cowdery (USGS, written commun., 2020) with a good quality designation that are not near the Mesabi Iron Range. 0.14 ft-1 T.K. Cowdery (USGS, written commun., 2020) head41 2011–2013 Barr (2014b) head targets from wells installed in the bedrock for 2011–2013 model. 0.2 ft-1 Barr (2014b) head201 2011–2013 Flooded mine pit elevations from the 2011 lidar, used in the 2011–2013 model. 0.1 ft-1 From lidar (MDNR, 2011)

32    Groundwater/Surface-Water Interactions in the Partridge River Basin and Hypothetical Future Mine Pits Table 2.  Weighting used for the target observation groups.—Continued [MDNR, Minnesota Department of Natural Resources; MN, Minnesota; %; percent; NWIS, National Water Information Sysytem (U.S. Geological Survey, 2019); USGS, U.S. Geological Survey; day/ft3, day per foot cubed; ft-1, per foot; lidar, light detection and ranging; UZF, unsaturated zone flow] Head targets head21 1995–2015 Hydraulic head targets representing data collected between 1995 and 2015 used in T.K. Cowdery (USGS, written commun., 2020) with a good quality designation that are not near the Mesabi Iron Range. 0.14 ft-1 T.K. Cowdery (USGS, written commun., 2020) head22 1995–2015 Hydraulic head targets representing data collected between 1995 and 2015 used in T.K. Cowdery (USGS, written commun., 2020) with a best quality designation that are not near the Mesabi Iron Range. 0.2 ft-1 T.K. Cowdery (USGS, written commun., 2020) head24 1995–2015 Hydraulic head targets representing data collected between 1995 and 2015 used in T.K. Cowdery (USGS, written commun., 2020) with a good quality designation that are near the Mesabi Iron Range. 0.04 ft-1 T.K. Cowdery (USGS, written commun., 2020) head25 1995–2015 Hydraulic head targets representing data collected between 1995 and 2015 used in T.K. Cowdery (USGS, written commun., 2020) with a best quality designation that are near the Mesabi Iron Range. 0.05 ft-1 T.K. Cowdery (USGS, written commun., 2020) head30 1995–2015 Barr (2014b) hydraulic head targets from wells installed in the unconsolidated material for 1995–2015 model. 0.1 ft-1 Barr (2014b) head31 1995–2015 Barr (2014b) hydraulic head targets from wells installed in the bedrock for 1995–2015 model. 0.1 ft-1 Barr (2014b) head40 1995–2015 Barr (2014b) head targets from wells in­ stalled in the unconsolidated material for 2011–2013 model. 0.2 ft-1 Barr (2014b) head300 1995–2015 Miscellaneous hydraulic head targets from the NWIS database that fall within 1995–2015 period. Head target site numbers can be found in the ancillary directory calibration files provided in the Haserodt and others (2021) data release associated with this publication. 0.1 ft-1 NWIS database head301 1995–2015 Continuous record hydraulic head targets from the NWIS database that fall within 1995–2015 period. Head target site numbers can be found in the ancillary directory calibration files provided in the Haserodt and others (2021) data release associated with this publication. 0.02 ft-1 NWIS database head200 1995–2015 Flooded mine pit elevations from the 2011 lidar data, used in the 1995–2015 model. 0.04 ft-1 From lidar (MDNR, 2011) wetlands 1995–2015 Composite targets looking at the percentage of areas that have groundwater seepage, as indicated by the UZF package, compared ar­ eas with the mapped wetlands. A target was created summing both areas with mapped wetlands but no seepage and areas with seepage but no mapped wetlands. 0.1 ft-1 No source; calculated with a utility that is run during the model calibration.

Model Calibration    33 directory calibration files provided in the Haserodt and others (2021) data release associated with this publication. There is one streamgage record from the MDNR Cooperative Stream Gaging database (://​www.dnr.st​ate.mn.us/​waters/​csg/​ index.; MDNR, 2018): streamgage 3149002 that is within the PRB, near the outlet of Colby Lake. Streamflow targets were processed to represent average base-flow conditions for the model period and these processing steps are described in appendix 3. Average pit inflow rates from pit pumping records for 2011 (Avery Guertin, MDNR, written commun., 2019) were available for two of the pumped mine pits (fig. 12; WEL cells in flooded pits) that were modeled using drain cells. The measured pit inflow rates were used as a target to com­ pare with the modeled pit inflow rates from the group of drain cells representing each of these pits in the 2011–2013 mining conditions model. Composite Wetland Area Two composite wetland area targets were used for calibration with the absence or presence of wetlands based on groundwater seepage distribution output from the UZF pack­ age relative to mapped peat deposits (fig. 6; Jirsa, 2016a, b). Mapped peat deposits were assumed to be indicative of peren­ nial wetland areas. A postprocessing Fortran utility (included in the Haserodt and others [2021] data release associated with this report) was developed to read in the UZF output from the 1995–2015 average conditions model during the calibration. The utility summed model cells with no groundwater seepage that were in mapped peat deposits (“false negative wetlands”) and cells with groundwater seepage that were in areas without mapped peat deposits (“false positive wetlands”). The areas from the utility program were used to calculate a percent­ age of false-negative wetland and false-positive wetland cells. Both summary percentages were used as targets with a measured value of 0 percent, which would represent a per­ fect match between mapped wetlands and the distribution of groundwater seepage. Weighting A total of 377 weighted targets were used in the calibra­ tion including the following: 9 streamflow targets adjusted to 2011–2013 base-flow conditions, 9 streamflow targets adjusted to 1995–2015 base-flow conditions, 2 pit inflow targets to compare with pit inflows in the 2011–2013 mining conditions model, 2 composite wetland area targets to compare with UZF results from the 1995–2015 average conditions model, and 355 hydraulic head targets (67 in the 2011–2013 mining conditions model and 288 in the 1995–2015 average conditions model). The targets were grouped based on data type, source, and quality and were given a weight that determines their importance for the calibration process (table 2). The goal of the weighting process was to transfer the information in the available measurements to the model parameters as discussed in Doherty and Hunt (2010). The initial weights were assigned to each of the observations based on uncertainty in the mea­ surements according to value, location, and data source. These weights were further adjusted with multipliers on the head and flux groups to balance the relative influence of the two data types during the calibration process. Additional minor adjustments were made to flux data based on quality and (or) importance in individual measurements. To calibrate the recharge multipliers, additional weight was given to the synoptic flow survey measurement at site 473059092064901 (U.S. Geological Survey, 2019), the streamflow target farthest downstream on the Partridge River but still upstream from Whitewater Lake/Colby Lake reservoir where reservoir affects and seepage at the base of the reservoir dam may complicate a streamflow target downstream from these lakes. The MDNR streamgage 3149002, located at the outlet of Colby Lake, was lower weighted due to concerns that some flow may be seep­ ing from the reservoir dam, which would not be accounted for in the model, and that reservoir affects may complicate the base flow separation for this record. Calibration Parameters The model calibration was performed with the parameter estimation code PEST_HP, version 16.1 (Doherty, 2018a, b, 2020). The goal of the calibration is to reduce predictive uncertainty by adjusting the model parameters until the mod­ eled values acceptably match the observations, without either overfitting to noise in the data or using geologically unreason­ able values for model parameters (honoring prior information about the system). Model parameters include inputs to the model that may be constrained to a range of values based on prior knowledge but are recognized as having some degree of uncertainty due to sparse measurement, measurement error, and structural errors inherent to the model (such as parameter simplification artifacts). The following model parameters were included in the calibration of the 2011–2013 mining condi­ tions and 1996–2015 average conditions models. Two multipliers on the gridded infiltration output from SWB were applied using the PEST utility TWOAR­ RAY. One multiplier was applied to upland areas and one for wetland areas, as determined by mapped peat deposits (fig. 6). Using two multipliers for calibration allowed the groundwater-flow model to separately compensate for possible bias in the SWB results in wetland compared to upland areas.

34    Groundwater/Surface-Water Interactions in the Partridge River Basin and Hypothetical Future Mine Pits The horizontal and vertical hydraulic conductivity in the zones in the unconsolidated (model layers 1–4) and bedrock unit zones (model layers 7–11). The horizontal and vertical hydraulic conductivity values at pilot points (allowing variability of hydrau­ lic conductivities within a zone; fig. 16A and B) in the upper bedrock layers (model layers 5 and 6). The flexibility provided by pilot points was included in model layers 5 and 6 because they have model param­ eters that are potentially informed by a greater concen­ tration of hydraulic head targets in the upper bedrock than in lower layers. The runoff coefficient for each of the flooded and pumped mine pits. This coefficient accounts for uncer­ tainty in variability in the drainage area surface and any errors in delineating drainage area divides in the uneven, mined terrain. The vertical hydraulic conductivity of the streambed material (part of the stream conductance term). A single value was used for all streams in the model. The drain conductance for each pumped pit. The drain conductance reflects the hydraulic conductivity of the pit bottom material and may be different for different pits. 1,300 1,400 1,500 1,600 1,700 1,800 Measured groundwater elevation, in feet above North American Vertical Datum 1988 1,300 1,400 1,500 1,600 1,700 1,800 Modeled groundwater elevation, in feet above North American Vertical Datum 1988 EXPLANATION [Mean error: – 0.33 Mean absolute error: 6.55 NWIS is National Water Information System] 1:1 line Good quality hydraulic head observations near mining features Mine pit levels from 2011 lidar NWIS miscellaneous hydraulic head observations Good quality hydraulic head observations away from mining features Best quality hydraulic head observations near mining features NWIS daily values hydraulic head observations Barr monitoring wells (Barr, 2014b) Best quality hydraulic head observations away from mining features Figure 17.  Modeled and measured heads for the 1995–2015 average conditions groundwater flow model of the Partridge River Basin, northeastern Minnesota. The diagonal dashed 1:1 line indicates a perfect fit. Targets in gray are assumed to have greater uncertainty due to proximity to mining features and larger target symbols within a color group represent measurements that are assumed to be higher quality. Included on the plot is the mean error and the mean absolute error (in feet) for the observations plotted.

Calibration Results and Discussion    35 A total of 384 adjustable parameters were included in the calibration. Two techniques were applied to assist the calibration process: singular value decomposition (SVD) and Tikhonov regularization (see Hunt and others [2007], Doherty and Hunt [2010], and Anderson and others [2015] for more details on the methods and benefits). Singular value decom­ position ensures that the parameter estimation problem can be solved; Tikhonov regularization allows the modeler to use prior knowledge of the parameters to restrain extreme param­ eter values that result from fitting noise in the observations (often called “overfitting”). The PEST utility, ADDREG1, was used to construct the Tikhonov regularization with preferred values equal to the initial values obtained from the Iron Range model from T.K. Cowdery (USGS, written commun., 2020). The target objective function used to control the degree of fitting (PEST variable PHIMLM) was set slightly above the lowest value of the overfit objective function without Tikhonov regularization. Calibration Results and Discussion The model was considered calibrated when the modeled and observed hydraulic heads and fluxes closely matched and the model parameters had reasonable values relative to litera­ ture ranges. Because the 2011–2013 mining conditions model was used to build the hypothetical mine scenarios, this is the focus of calibration results section. Calibration residuals for hydraulic head and flux targets are shown for the 2011–2013 mining conditions models and for hydraulic head targets in the 1995–2015 average conditions model, which had most of the Measured groundwater elevation, in feet above North American Vertical Datum 1988 Modeled groundwater elevation, in feet above North American Vertical Datum 1988 1,400 1,450 1,500 1,550 1,600 1,650 1,700 1,750 1,400 1,450 1,500 1,550 1,600 1,650 1,700 1,750 1:1 line EXPLANATION [Mean error: – 0.91 Mean absolute error: 5.66] Barr monitoring wells (Barr, 2014b) Mine pit levels from 2011 lidar Best quality hydraulic head observations away from mining features Best quality hydraulic head observations near mining features Good quality hydraulic head observations away from mining features Figure 18.  Modeled and measured heads for the 2011–2013 mining conditions groundwater flow model of the Partridge River Basin, northeastern Minnesota. The diagonal dashed 1:1 line indicates a perfect fit. Targets in grey are assumed to have greater uncertainty due to proximity to mining features and larger target symbols within a color group represent measurements that are assumed to be higher quality. Included on the plot is the mean error and the mean absolute error (in feet) for the observations plotted.

36    Groundwater/Surface-Water Interactions in the Partridge River Basin and Hypothetical Future Mine Pits hydraulic head targets. Calibrated model parameters were the same in both models, except for the calibrated net infiltration, which is discussed for just the 2011–2013 mining conditions model. The fit between measured (observed) and modeled hydraulic heads and fluxes is shown on 1:1 graphs where the measured values are plotted against the modeled. Figures 17–19 show 1:1 graphs for the 1995–2015 hydraulic head targets, 2011–2013 hydraulic head targets, and 2011–2013 streamflow targets, respectively. The mean error indicates if there is bias high or low in the modeled values relative to the observed data. The mean absolute error provides information on the overall magnitude of these differences. The 1995–2015 (fig. 17) and 2011–2013 hydraulic head residuals (fig. 18) showed a small, negative mean error, meaning there is a slight bias of modeled values exceeding the observed values. High and low hydraulic head and flux residuals were distributed across the domain without areas of substantial bias (figs. 20 and 21). Hydraulic head residuals were gener­ ally largest near abrupt changes in topography related to mining features and along the northern boundary of the model domain. Base flows were slightly undersimulated near the farthest upstream headwaters of the Partridge River and oversimulated in some of the downstream tributaries. The downstream target at site 473059092064901 slightly undersimulates base flow whereas the MDNR streamgage (3149002) is appreciably oversimulated. In the parameter estimation process, the MDNR downstream streamgage (3149002) was assigned a lower weight because its location downstream from a dam confounds accurate estimation of the groundwater-derived base-flow component of measured streamflow. The low weight for parameter estimation, in turn, results in other targets obtaining better fits at the expense of this more uncertain target. Calibrated horizontal hydraulic conductivity and vertical anisotropy (horizontal ÷ vertical hydraulic conductivity) for each of the 11 model layers are shown in figures 22 and 23 and presented in table 3. The unconsolidated horizontal hydraulic conductivities ranged from 0.1 to 100 ft/day for the glacial materials and peat and from 0.08 to 219 ft/day for mining fea­ tures. The unconsolidated horizontal hydraulic conductivity is most variable in layer 1 due to the presence of mining, glacial, and peat deposits (fig. 22A). Layers 2–4 (fig. 22B-22D) are below the peat and include the glacial deposits and various mining features like pits, stockpiles, and the tailings basin. Bedrock (model layers 5–11; fig 22E-22K) showed decreasing horizontal hydraulic conductivity with depth for each bedrock unit with hydraulic conductivities that ranged from 3x10-6 to 50 ft/day. Vertical anisotropy for each of the model layers (fig. 23 and table 3) should generally be greater than 1.0 because potential layering of lower conductivity materials favors 3149002 473059092064901pr01p 4015447partridgernrb Measured streamflow, in cubic feet per second Modeled streamflow, in cubic feet per second EXPLANATION [Largest targets named] 1:1 line Zero-weighted streamflow target 2011 to 2013 streamflow target Figure 19.  Comparison of modeled and measured streamflow for the 2011–2013 mining conditions groundwater flow model of the Partridge River Basin, northeastern Minnesota. The three largest streamflow targets show their station names and can be compared to the map of streamflow target locations in figure 21. The diagonal, dashed red 1:1 line indicates a perfect fit. The gray cross symbol is for a zero-weighted flux target; this synoptic measurement flow measurement on the Dunka River was zero-weighted because it is not fully within the model domain.

Calibration Results and Discussion    37 Hillshade base from Minnesota Department of Natural Resources, 2011 Hydrologic features modified from Minnesota Department of Natural Resources, 2019 Partridge River Basin outline from U.S. Geological Survey, 2012 NHDPlus stream lines from McKay and others, 2012 Mine features modified from John Coleman, GLIFWC, written common., January 29,2018 47°30' 47°36' 47°42' 92°12' 92°6' 92°0' 91°54' 91°48' 8 MILES 6 KILOMETERS

EXPLANATION Lakes represented with the RIV package Existing mine pits Partridge River Basin Model domain Streams in NHDPlus Permitted groundwater withdrawals (2011) Hydraulic head targets simulated too high (measured modeled) – 5.0 to 0.0 – 0.0 to 5.0 – 20.0 to <–10.0 – 30.0 to 20.0 Hydraulic head targets simulated too low (measured modeled) 0.0 to 5.0 5.0 to 10.0 10.0 to 20.0 20.0 to 30.0 Figure 20.  Spatial distribution of hydraulic head target residuals for the 1995–2015 average conditions groundwater-flow model of the Partridge River Basin, northeastern Minnesota. Residuals in red indicate areas the model is undersimulating groundwater elevations and blue residuals indicate areas the model is oversimulating groundwater elevations. The symbol size is scaled to reflect the magnitude of the residual with large symbols representing larger residuals.

38    Groundwater/Surface-Water Interactions in the Partridge River Basin and Hypothetical Future Mine Pits

EXPLANATION Streamflow, in cubic feet per second 0 to 0.5 0.5 to 1 1 to 2.5 2.5 to 5 10 to 15 15 to 20 20 to 27.1 Partridge River Basin Model domain Streamflow residuals, in cubic feet per second –20 to –5 –5 to –2 –2 to –0.5 –0.5 to 0 0 to 0.5 0.5 to 2 Basemap from , 2015 Partridge River Basin outline from U.S. Geological Survey, 2012 S t

o ui s

Riv e r P a rt r d g e R e r Bi rc h R e r E m ba rrass R e r 8 KILOMETERS 8 MILES Dun k a

Ri e r 47°42' 47°36' 47°30' 92°12' 92°6' 92°0' 91°54' 91°48' 473059092064901 473203092065701 473335092020501 473147091595201 473032091563601 –1.9 – 0.7 –1.6 –19.2 Figure 21.  Spatial distribution of streamflow flux target residuals and stream accumulation, for the 2011–2013 mining conditions groundwater-flow model of the Partridge River Basin, northeastern Minnesota. Greater streamflow is indicated by larger symbols of darker color. Streamflow flux target residuals show the difference between measured and modeled values, in cubic feet per second. Red targets are undersimulated, where the modeled value is less than the measured value. Blue targets are oversimulated, where the modeled value is greater than the measured value. The flux residual, in cubic feet per second, is shown above each target and the site number is shown below each target.

Calibration Results and Discussion    39 lateral movement of water more than vertical movement. Vertical anisotropy was generally greater than 1.0 but less than 20 in the glacial materials, with higher values in the peat, the Brainerd Sublobe silt and clay, and the tailing basin deposits. The vertical anisotropy in the bedrock varied from less than 1 at a few pilot points for the Duluth Complex and Archean units in layer 5 to 440 at a pilot point for the Biwabik Iron-Formation in layer 5. Results from the parameter estima­ tion indicated that the Biwabik Iron-Formation had the highest bedrock vertical anisotropy in all bedrock layers, which is consistent with it being the only metasedimentary unit. The calibrated horizontal hydraulic conductivity values for each model zone representing the same type of uncon­ solidated material (for example, sand and gravel) are shown with ranges of literature values in figure 24. In general, the calibrated horizontal hydraulic conductivities fall within or close to the literature ranges for each material, although no literature ranges were available for the fine tailings, stockpiles, or the area mapped as disturbed by mining by Jirsa, (2016 a, b). The “disturbed by mining” zone (fig. 5) surrounds the iron mining features and represents unconsolidated materi­ als that have been moved or disturbed by mining activities. This disturbed zone hydraulic conductivity is important for local groundwater flow from the unconsolidated aquifer to the mine pits. The calibrated value of 27 ft/day for this disturbed zone fell within the literature range for sand (0.06–1,700 ft/ day; Schwartz and Zhang, 2003). The hydraulic conductivity for the bedrock layers with depth is shown in figure 25 for each of the four bedrock units–Duluth Complex, Virginia Formation, Biwabik IronFormation, and Archean units. Layers 5 and 6 in the bedrock were calibrated using pilot points, so average values are presented in figure 25. The field-measured and modeled values from Barr (2014a) and Myers (2014) represent the range in values for bedrock data collected or modeled near the domain. For the Archean units, Duluth Complex, and Virginia Formation, the upper three bedrock layers, which repre­ sent 200 ft of bedrock, have calibrated horizontal hydraulic conductivity values that plot within the literature range for fractured igneous and metamorphic rocks and the lower four bedrock layers plot at or just below the unfractured bedrock range. The Biwabik Iron-Formation has slighter higher hori­ zontal hydraulic conductivity values at all depths. Calibrated potential recharge used two calibration multi­ pliers of the SWB net infiltration array, one for upland and one for wetland areas, delineated by the presence of peat deposits (fig. 6). The 2011–2013 and the 1995–2015 net infiltration arrays used the same calibration multipliers. Calibrated SWB net infiltration for the 2011–2013 mining conditions model is Table 3.  Calibrated horizontal hydraulic conductivity and vertical anisotropy (defined as the horizontal hydraulic conductivity/vertical hydraulic conductivity) values for the unconsolidated and bedrock materials in the Partridge River Basin, northeastern Minnesota. A. Unconsolidated deposits

Estimated horizontal hydraulic conductivity, in feet per day Horizontal:vertical hydraulic conductivity, unitless Unit Rainy Lobe diamicton. Brainerd Sublobe diamicton. Rainy and Superior Lobe diamicton. Rainy Lobe silt and clay. Brainerd Sublobe silt and clay. Rainy Lobe or nonclassified sand. Brainerd Sublobe sand. Superior Lobe sand. Rainy Lobe gravel. Brainerd Sublobe gravel. Rainy and Superior Lobe sand and gravel. Superior Lobe gravel. Peat. Stockpiles. Disturbed by mining. Tailings basin fines.

40    Groundwater/Surface-Water Interactions in the Partridge River Basin and Hypothetical Future Mine Pits discussed in detail with the UZF results in the Model Results and Discussion section. The average calibrated SWB infil­ tration across the model domain for the 2011–2013 mining conditions model was 5.3 in/yr. Calibrated runoff coefficients for the mine pits ranged from 0.38–0.80, with an average of 0.43. The calibrated vertical hydraulic conductivity of the streambed material was 1.00 ft/day. The average drain conductance for the pumped pits was 1.97x105 square feet per day (ft2/day). Using a width and length of the model cell (200 ft) and an assumed pit bed thickness of 1 ft, this drain conductance represents a vertical hydraulic conductivity of 4.9 ft/day. Table 3.  Calibrated horizontal hydraulic conductivity and vertical anisotropy (defined as the horizontal hydraulic conductivity/vertical hydraulic conductivity) values for the unconsolidated and bedrock materials in the Partridge River Basin, northeastern Minnesota.— Continued B. Bedrock units

Estimated horizontal hydraulic conductivity, in feet per day Horizontal:vertical hydraulic conductivity, unitless Model layer Duluth Complex of Keweenan Supergroup 0.02–10 0.08–39 0.027–0.034 36–49 3.E-03 4.E-04 4.E-05 4.E-05 3.E-06 Virginia Formation 0.6–14 1–26 0.03–0.04 76–88 3.E-03 3.E-04 3.E-05 3.E-05 3.E-06 Biwabik Iron-Formation of Animikie Group 4–50 12–440 0.2–0.4 53–100 8.E-03 4.E-03 3.E-03 4.E-04 Archean units 1–13 0.4–7 0.1–0.2 6–14 4.E-03 4.E-04 4.E-04 4.E-05

Calibration Results and Discussion    41 Model row Model row A. Layer 1 horizontal hydraulic conductivity Model row Model row Log horizontal hydraulic conductivity, in feet per day EXPLANATION B. Layer 2 horizontal hydraulic conductivity C. Layer 3 horizontal hydraulic conductivity D. Layer 4 horizontal hydraulic conductivity Figure 22.  Horizontal hydraulic conductivity in each of the 11 model layers for the groundwater-flow model of the Partridge River Basin, northeastern Minnesota. Note, the figure scale is in log scale but the values presented on the figure are actual values that are not log transformed.

42    Groundwater/Surface-Water Interactions in the Partridge River Basin and Hypothetical Future Mine Pits Model row Model row E. Layer 5 horizontal hydraulic conductivity Model row Model row 1e –6 1e –5 1e –4 1e –3 1e –2 F. Layer 6 horizontal hydraulic conductivity G. Layer 7 horizontal hydraulic conductivity H. Layer 8 horizontal hydraulic conductivity Horizontal hydraulic conductivity, in feet per day EXPLANATION Figure 22.—Horizontal hydraulic conductivity in each of the 11 model layers for the groundwater-flow model of the Partridge River Basin, northeastern Minnesota. Note, the figure scale is in log scale but the values presented on the figure are actual values that are not log transformed.

Calibration Results and Discussion    43 Model row Model row I. Layer 9 horizontal hydraulic conductivity Model row Model row 1e –6 1e –5 1e –4 1e –3 1e –2 J. Layer 10 horizontal hydraulic conductivity K. Layer 11 horizontal hydraulic conductivity Horizontal hydraulic conductivity, in feet per day EXPLANATION Figure 22.—Horizontal hydraulic conductivity in each of the 11 model layers for the groundwater-flow model of the Partridge River Basin, northeastern Minnesota. Note, the figure scale is in log scale but the values presented on the figure are actual values that are not log transformed.

44    Groundwater/Surface-Water Interactions in the Partridge River Basin and Hypothetical Future Mine Pits Model row Model row A. Layer 1 vertical anisotrophy Model row Model row B. Layer 2 vertical anisotrophy C. Layer 3 vertical anisotrophy D. Layer 4 vertical anisotrophy Vertical anisotrphy, unitless EXPLANATION Figure 23.  Vertical anisotropy (defined as the horizontal hydraulic conductivity/vertical hydraulic conductivity) in each of the 11 model layers for the groundwater flow model of Partridge River Basin, northeastern Minnesota.

Calibration Results and Discussion    45 Model row Model row E. Layer 5 vertical anisotrophy Model row Model row F. Layer 6 vertical anisotrophy G. Layer 7 vertical anisotrophy H. Layer 8 vertical anisotrophy Vertical anisotrphy, unitless EXPLANATION Figure 23.—Vertical anisotropy (defined as the horizontal hydraulic conductivity/vertical hydraulic conductivity) in each of the 11 model layers for the groundwater flow model of Partridge River Basin, northeastern Minnesota.

46    Groundwater/Surface-Water Interactions in the Partridge River Basin and Hypothetical Future Mine Pits I. Layer 9 vertical anisotrophy Model row Model row Model row Model row J. Layer 10 vertical anisotrophy K. Layer 11 vertical anisotrophy Vertical anisotrphy, unitless EXPLANATION Figure 23.—Vertical anisotropy (defined as the horizontal hydraulic conductivity/vertical hydraulic conductivity) in each of the 11 model layers for the groundwater flow model of Partridge River Basin, northeastern Minnesota.

Model Results and Discussion    47 Model Results and Discussion This section discusses the model results from the 2011–2013 mining conditions model. This model and the 1995–2015 average conditions model (Haserodt and others, 2021) were used for parameter estimation; however, only the 2011–2013 model results are discussed in this report because they represented mining conditions that are of primary impor­ tance for simulations of future mining scenarios. Groundwater inflows to and outflows from the model (model groundwater budget) are shown in figure 26. Sources of groundwater inflow included recharge, regional groundwater coming through the model boundary, wells injecting water into mine pits to represent runoff into pits,and streams and lakes that were losing water to the aquifer. Major outflows from the groundwater system were groundwater seepage to streams, wetlands, and ponds at land surface; dis­ charge to streams and lakes; seepage into pumped mine pits; and water leaving through the model boundaries. The simulated water table across the PRB is shown in figure 27. Groundwater flows from areas with high hydraulic heads to lower hydraulic heads approximately perpendicular to the equal-elevation contour lines. The topographic high along the Iron Range is a groundwater divide in the northern bound­ ary of the PRB. Groundwater flow in the PRB was generally southeast from the Iron Range and west from the topographic high in the eastern PRB. The groundwater flow reflects a sys­ tem in which water moves from topographic highs and other recharge areas to streams and low-lying wetland areas (fig. 27). Streams, represented by SFR2 cells (fig. 28), can either be gaining water from the groundwater system (gaining stream) or losing water to the groundwater system (losing stream). As expected in a humid temperate climate without appreciable pumping, streams are mostly gaining in the PRB model domain, with the exception of reaches near the pumped pits and a few reaches near the model perimeter where bound­ ary conditions may be influencing local groundwater/surfacewater interaction. Fine tailings Silt and clay Peat Diamicton Sand Gravel Disturbed by mining Stockpiles 10–2 10–1 Horizontal hydraulic conductivity, in feet per day Geologic unit Barr (2014b) upland deposits calibrated values Literature range from Schwartz and Zhang (2003) Barr (2014b) wetland deposits calibrated values Calibrated values EXPLANATION Figure 24.  Calbrated hydraulic conductivity values (black diamond) and expected literature ranges (blue, red, and gray rectangles) for the unconsolidated materials within the model domain. Literature values came from range of upland and wetland deposits in the Barr (2014b) model and from values in Schwartz and Zhang (2003). Calibrated values are also presented in table 3.

48    Groundwater/Surface-Water Interactions in the Partridge River Basin and Hypothetical Future Mine Pits As described in the “Boundary Conditions” section, MODFLOW’s UZF package partitions the infiltration cal­ culated by the SWB model into groundwater recharge and groundwater runoff, where groundwater runoff is composed of infiltration rejected because of the lack of subsurface stor­ age and seepage at land surface resulting from groundwater heads above the land surface. The groundwater runoff is either routed or removed from the model. Thus, stream base flow can be derived from both groundwater runoff and direct discharge to the stream channel. Large parts of the PRB were characterized by having a shallow water table, peat, and a thin mantle of unconsolidated material overlying bedrock. As a result, a substantial part of the simulated streamflow within the PRB was derived from groundwater runoff calculated by the UZF package. Figure 29 shows the distribution of net infiltration estimated by SWB, SWB net infiltration available to UZF as potential recharge after the calibration multipliers were applied, the amount of recharge applied by UZF to the water table, and rejected recharge. The average, calibrated SWB infiltration across the model domain for the 2011–2013 min­ ing conditions model was 5.3 in/yr and the average applied UZF recharge (water that crosses the water table after rejected recharge is removed) was 4.6 in/yr. Figure 30 breaks the applied UZF recharge into recharge on upland and wetland areas. On average, uplands received 6.0 in/yr of recharge and wetlands received 2.6 in/yr of recharge. The UZF package accounting of recharge, rejected recharge, and groundwater seepage allowed the relative impor­ tance of these components to be evaluated. For the 2011–2013 simulation, UZF rejected 12.9 percent of the calibrated SWB infiltration (fig, 31A), and of this rejected water, 42.0 percent was then routed to a stream and 58.0 percent was removed from the model (fig. 31B). For the groundwater seepage com­ ponent, 30.7 percent was routed and 69.3 percent was removed from the model via assumed evapotranspiration (fig. 31C). These components were evaluated in terms of importance for simulated stream base flow. For 2011–2013 simulated base flow, 50.4 percent comes from direct groundwater discharge to the channel, 37.5 percent is from routed groundwater seepage, and 12.1 percent is from routed rejected recharge (fig. 31D). Approximately one-half the simulated base flow was derived from groundwater discharge directly to streams and one-half was contributed by UZF-routed groundwater runoff.Spatially, groundwater seepage was simulated across 33 percent of the model domain (fig. 32) and typically occurred in topographic lows, where peat deposits were common (Jirsa, 2016b). Hypothetical Mine Pit Scenarios and Model Forecasts The calibrated 2011–2013 mining conditions model was used to create five hypothetical mine pit model scenarios (fig. 33) for known mineral deposits (MDNR, 2016) within the basin. The spatial extent of these deposits was digitized from the MDNR (2016) published map, “Exploration for metallic mineral resources in Minnesota Duluth Complex.” Most of these deposits have not entered the mine permit phase and some may never be developed. These model scenarios were not designed to predict affects from any specific future mine within the basin; to do so would require a groundwater-flow model with detailed information about the proposed mine and calibration data near the mine. The purpose of these mining scenarios was to present a possible workflow to quantify model uncertainty for a given model forecast and serve as a potential guide for initial data collec­ tion that may improve a future model’s ability to make such a forecast. Looking at possible mining scenarios across the basin highlighted how the variability in the groundwater flow system, the distribution of current calibration data, and the location of model forecasts affected the reliability with which the model can make the predictions of interest. The location of the hypothetical mine pits and the flooded extent of the Peter Mitchell pit are shown in figure 33. The names and features in the hypothetical mine pit scenarios are summarized in table 4. The hypothetical mine pits were modeled using drain cells with a conductance of 2x105 ft2/day, based on the cali­ brated drain conductance values of 1.9x105–2.0x105 ft2/day. The drain elevations for all the hypothetical mine pits were set to 660 ft below the model top elevation (land surface). The deposit depths are unpublished for most of the deposits the hypothetical mine pits were targeting; therefore, a depth of 660 ft was used based on the proposed depths for the two deeper NorthMet pits in the PRB (MDNR, 2013). Because the Peter Mitchell pit is an important hydrologic feature in the PRB headwaters for the 2011–2013 mining con­ ditions model, it was modified to represent future conditions in the mine pit scenario simulations. If developed, future mines in the basin would likely be operational as the Peter Mitchell pit is phasing out or after it has been reclaimed and flooded. Therefore, the Peter Mitchell pit was modeled in the scenarios as a flooded mine pit with a water elevation of 1,500 ft and a bottom topography digitized from Barr (2010). This repre­ sented the expected final buildout for the Peter Mitchell pit with the flooded elevation set to the longterm expected water level for the mine pit (Barr, 2010). With a water elevation of 1,500 ft, much of the Peter Mitchell pit will be lower than groundwater elevations in some nearby areas and the pit will act as a hydrologic sink. Near the eastern end where local groundwater elevations are lower than 1,500 ft, however, it may act as a source to the groundwater system. Therefore, the flooded Peter Mitchell pit was modeled using MODFLOW’s General Head Boundary (GHB) condition with a specified

Hypothetical Mine Pit Scenarios and Model Forecasts    49 Horizontal hydraulic conductivity, in feet per day 10–1 10–3 10–5 10–7 10–9 Model layer A. Archean units 10–1 10–3 10–5 10–7 10–9 Horizontal hydraulic conductivity, in feet per day Model layer B. Biwabik Iron-Formation of Animikie Group Range of horizontal hydraulic conductivit Barr (2014b) and Myers (2014) literatur field values lowest est. at 1.4e-6 ft/da Divide between literature fractured and unfractured bedrock ranges Duffield (2019) fractured and unfractured igneous and metamorphic rocks EXPLANATION [est. means estimated] Range of horizontal hydraulic conductiv from Barr (2014b) and Myers (2014) literature and field values Divide between literature fractured and unfractured bedrock ranges Duffield (2019) fractured and unfractured igneous and metamorphic rocks EXPLANATION Figure 25.  Calibrated hydraulic conductivity values (blue bars) and expected literature ranges for each model layer representing the bedrock (model layers 5–11) in the Duluth Complex of Keweenan Supergroup, Virginia Formation, Biwabik Iron-Formation of Animikie Group, and Archean units. Included on the plots are a red rectangle representing the range in modeled and field data from Myers (2014) and Barr (2014a) and a gray rectangle representing the literature values presented by Duffield (2019) for fractured and unfractured igneous and metamorphic rocks. A dashed line shows a divide between the lowest literature (Duffield, 2019) values for fractured bedrock and the upper value for the unfractured bedrock. Calibrated bedrock hydraulic conductivity values are also presented in table 3.

50    Groundwater/Surface-Water Interactions in the Partridge River Basin and Hypothetical Future Mine Pits Horizontal hydraulic conductivity, in feet per day 10–1 10–3 10–5 10–7 10–9 Model layer C. Virginia Formation 10–1 10–3 10–5 10–7 10–9 Horizontal hydraulic conductivity, in feet per day Model layer D. Duluth Complex of Keweenan Supergroup Range of horizontal hydraulic conductivity from Barr (2014b) and Myers (2014) literature and field values Divide between literature fractured and unfractured bedrock ranges Duffield (2019) fractured and unfractured igneous and metamorphic rocks EXPLANATION Range of horizontal hydraulic conductivity from Barr (2014b) literature and field values Divide between literature fractured and unfractured bedrock ranges Duffield (2019) fractured and unfractured igneous and metamorphic rocks EXPLANATION Figure 25.—Figure 25.  Calibrated hydraulic conductivity values (blue bars) and expected literature ranges for each model layer representing the bedrock (model layers 5–11) in the Duluth Complex of Keweenan Supergroup, Virginia Formation, Biwabik Iron-Formation of Animikie Group, and Archean units. Included on the plots are a red rectangle representing the range in modeled and field data from Myers (2014) and Barr (2014a) and a gray rectangle representing the literature values presented by Duffield (2019) for fractured and unfractured igneous and metamorphic rocks. A dashed line shows a divide between the lowest literature (Duffield, 2019) values for fractured bedrock and the upper value for the unfractured bedrock. Calibrated bedrock hydraulic conductivity values are also presented in table 3.

Model Forecast Results and Associated Uncertainty    51 elevation of 1,500 ft (set in the layer where the excavated pit bottom would occur) and with a conductance of 2x105 ft2/day. GHB cells were selected instead of DRN cells to allow for both groundwater inflow and outflow to the pit. For both the hypothetical mine pits and the flooded Peter Mitchell Pit, a high hydraulic conductivity (1,000 ft/day in both horizontal and vertical) was used in cells that repre­ sented air (excavated rock) or water in the pits. For cells in the flooded Peter Mitchell pit where the land surface in the 2011 lidar data were below the flooded elevation of 1,500 ft, the model top elevation was adjusted to just above 1,500 feet to prevent unrealistic circulation between flooded-pit GHB cells and UZF. The UZF routing array was also modified to remove infiltration from the nodes of the hypothetical mine pit areas. The 2011–2013 net precipitation of 1.4 in/yr was applied as net infiltration over the hypothetical mine pits. Although this rate was consistent with the recharge applied over flooded pits in the base model, no attempt was made to estimate surface drainage areas for the new pits because the geometry of the potential mine pit areas was unknown. Streams intersecting the hypothetical mine pits remained in the model in their original configuration and given a very low conductance, effectively modeling them as an imperme­ able conduit that transferred water from the upstream side of the new pit to the downstream side. This conduit method assumes that a stream corridor affected by a mine pit would be rerouted around the open pit. The premining stream con­ figuration was retained, however, because the location of the rerouted stream was unknown. Additionally, the current mine discharge of 2.3 cubic feet per second of pumped pit water into the headwaters of the Partridge River was removed to reflect the assumption that dewatering would not be required after the Peter Mitchell pit floods. Model Forecast Results and Associated Uncertainty Model forecasts were selected to assess the effects of hypothetical mine pits on the groundwater and surface-water systems. Forecasts were selected in conjunction with Great Lakes Indian Fish & Wildlife Commission and included the following: the average depth to groundwater across 12 mapped wetlands in the model domain (fig. 33); streamflow at 6 locations in the PRB (fig. 33); the pit inflow to the Peter Mitchell pit in the final iron mining scenarios and each of the 4 new mine pit sce­ narios (table 4); and the pit inflow to new pits in the 4 mining scenarios (fig. 33). – 9.9 –10.3 –18.0 – 33.3 – 0.4 – 80.6 Wells (boundary and pit runoff) Pumped mine pits Lakes Streams Recharge Groundwater seepage Water supply wells Budget item –100 – 50 Flux in or out, in cubic feet per second Inflows to groundwater system Outflows from groundwater system EXPLANATION Figure 26.  Mass balance for the 2011–2013 mining conditions groundwater-flow model of the Partridge River Basin, northeastern Minnesota. Positive volumetric fluxes are shown in blue and represent inflows to the groundwater system and negative flows are shown in red and represent outflows from the groundwater system. The water supply wells had such a low flux it is not visible in this graph.

52    Groundwater/Surface-Water Interactions in the Partridge River Basin and Hypothetical Future Mine Pits , , , , , , ,50 , , 1,5 , , , ,7 , , , , , , , 47°42' 47°36' 47°30' 92°12' 92°6' 92°0' 91°54' 91°48' Hillshade base from Minnesota Department of Natural Resources, 2011 NHDPlus stream lines from McKay and others, 2012 , , , , ,, 6 KILOMETERS 6 MILES EXPLANATION Pumped pits represented with the DRN package Lakes represented with the RIV package Simulated water-table elevation — Shows altitude of water table. Contour interval is 1,00 feet. Datum is NAVD88 Simulated water-table elevation — Shows altitude of water table. Contour interval is 25 feet. Datum is NAVD88 Model domain Streams represented using the Streamflow Routing (SFR2) package 1,700 1,625 Figure 27.  Simulated water-table elevation for the 2011–2013 mining conditions groundwater flow simulation.

Model Forecast Results and Associated Uncertainty    53 47°30' 47°36' 47°42' 92°12' 92°6' 92°0' 91°54' 91°48' Hillshade base from Minnesota Department of Natural Resources, 2011 Partridge River Basin Outline from U.S. Geological Survey, 2012 NHDPlus stream lines from McKay and others, 2012 6 MILES 6 KILOMETERS EXPLANATION Pumped pits represented with the DRN package Lakes represented with the RIV package Partridge River Basin Model domain Gaining and losing stream reaches Stream cells losing water to the aquifer Stream cells gaining water from the aquifer Pumping wells with permitted groundwater withdrawals Figure 28.  Simulated losing and gaining reaches of the Partridge River for the 2011–2013 mining conditions model.

54    Groundwater/Surface-Water Interactions in the Partridge River Basin and Hypothetical Future Mine Pits Model column (Mean: 0.7 inches per year, Max: 12.8 inches per year, Min: 0.0 inches per year) D. Rejected recharge Model column Model row (Mean: 4.6 inches per year, Max: 15.2 inches per year, Min: 0.0 inches per year) C. UZF recharge (Mean: 5.3 inches per year, Max: 15.2 inches per year, Min: 0.0 inches per year) B. Calibrated SWB recharge Model row (Mean: 5.7 inches per year; Max: 13.6 inches per year, Min: 0.0 inches per year) A. Original SWB recharge Recharge rate, in inches per year EXPLANATION [SWB soil water balance UFZ unsaturated zone flow max. maximum min. minimum] Figure 29.  A, The spatial distribution of net infiltration estimated by soil-water-balance (SWB) package , B, the SWB calibrated net infiltration after model calibration, C, the UZF recharge (SWB calibrated net infiltration minus rejected recharge), and, D, simulated rejected recharge for the 2011–2013 simulation.

Model Forecast Results and Associated Uncertainty    55 Model column Model column Model row Mean: 6.0 inches per year (lakes and mine pits removed) C. UZF recharge in uplands Mean: 2.6 inches per year B. UZF recharge in wetlands Model row Mean: 4.6 inches per year A. UZF recharge Recharge rate, in inches per year EXPLANATION [UFZ unsaturated zone flow] Figure 30.  The spatial distribution of simulated recharge (soil-water-balance calibrated net infiltration minus rejected recharge), as it is applied to wetland areas and upland areas, for the 2011–2013 simulation.

56    Groundwater/Surface-Water Interactions in the Partridge River Basin and Hypothetical Future Mine Pits Percentage of total groundwater seepage from unsaturated zone flow package that is routed to streams Unsaturated zone flow recharge 87.1% Rejected recharge 12.9% A. Percentage (%) of calibrated recharge from SWB that is rejected recharge and percent that is applied as unsaturated zone flow recharge D. Percentage of total streamflow inflows from groundwater discharge to streams and from routed seepage and rejected recharge from the unsaturated zone flow package Total streamflow from routed unsaturated zone flow seepage 37.5% Total streamflow from groundwater discharge 50.4% Total streamflow from routed rejected recharge 12.1% Total non-routed unsaturated zone flow seepage 69.3% Total routed unsaturated zone flow seepage 30.7% B. Percentage of total rejected recharge that is routed to streams Total non-routed unsaturated zone flow rejected recharge 58.0% Total routed unsaturated zone flow rejected recharge 42.0% Figure 31.  A, how much available recharge is rejected by unsaturated zone flow , B, percentage of the rejected recharge is routed to streams, C, percentage of the unsaturated zone flow seepage is routed to streams, and, D, percentage that routed groundwater runoff contributes to total streamflow for the 2011–2013 mining conditions MODFLOW-NWT model of the Partrigde River Basin, northeastern Minnesota.

Model Forecast Results and Associated Uncertainty    57 47°30' 47°36' 47°42' 92°12' 92°6' 92°0' 91°54' 91°48' Peat deposits from Jirsa, 2016b Hillshade base from Minnesota Department of Natural Resources, 2011 8 MILES 8 KILOMETERS Peat deposits Model domain EXPLANATION Groundwater seepage, in cubic feet per day Figure 32.  Groundwater seepage from unsaturated zone flow for the 2011–2013 mining conditions MODFLOW-NWT model of the Partridge River Basin, northeastern Minnesota.

58    Groundwater/Surface-Water Interactions in the Partridge River Basin and Hypothetical Future Mine Pits 47°30' 47°36' 47°42' 92°12' 92°6' 92°0' 91°54' 91°48' Base map data from , 2015 and Minnesota Department of Natural Resources, 2011 Wetlands modified from U.S. Fish and Wildlife Service, 2018 Hydrologic features modified from Minnesota Department of Natural Resources, 2019 8 MILES 8 KILOMETERS EXPLANATION Lakes Wetland forecast locations and forecast number and forecast number Wetlands Hypothetical mine pit and scenario number Peter Mitchell pit at final buildout 1,500 foot flooded footprint Streams Model domain Streamflow flux forecast locations Figure 33.  Hypothetical mine pits and model forecasts locations for 12 wetlands and 6 stream reaches in the Partridge River Basin, northeastern Minnesota.

Model Forecast Results and Associated Uncertainty    59 Scenario 1 Scenario 2 Scenario 3 Scenario name Scenario 4 Final iron mining Groundwater flux into the pit, in cubic feet per second Predicted value Estimated uncertainty (1 standard deviation) EXPLANATION Figure 34.  Model forecast values and associated uncertainty for the pit inflow rates to the Peter Mitchell pit at final buildout without new mining and then in the scenarios with each of the four hypothetical mine pits. The solid bars in the plots represent the modeled value of the forecast and the error bars represent one standard deviation above and below that value. Scenario 1 Scenario 2 Scenario 3 Scenario name Scenario 4 Groundwater flux into the pit, in cubic feet per second Predicted value Estimated uncertainty (1 standard deviation) EXPLANATION Figure 35.  Model forecast values and associated uncertainty for the pit inflow rates to each of the hypothetical mine pits in scenarios 1–4. The solid bars in the plots represent the modeled value of the forecast and the error bars represent one standard deviation above and below that value.

60    Groundwater/Surface-Water Interactions in the Partridge River Basin and Hypothetical Future Mine Pits Steamflow at forecast location 1 Steamflow at forecast location 2 Steamflow, in cubic feet per second Steamflow at forecast location 3 Steamflow at forecast location 4 Steamflow at forecast location 5 Steamflow at forecast location 6 Scenario name Final iron mining Base model Final iron mining Base model Predicted value: groundwater discharge to steam + routed groundwater runoff Predicted value: pumped mine pit water Estimated uncertainty (1 standard deviation) EXPLANATION Figure 36.  Model forecast values and associated uncertainty for the six streamflow forecasts in the final iron mining scenario, four new mine scenarios, and 2011–2013 base model. Streamflow represents base flow conditions and graphs are ordered from up- to downstream along the Partridge River. The 2011–2013 base model included a 2.3 cubic feet per second stream input from water that is pumped out of a mine pit. This pit is then part of the flooded Peter Mitchell pit in the final iron mining and mine scenarios 1–4; the pit pumping is assumed to cease for the scenarios. The solid bars in the plots represent the modeled value of the forecast and the error bars represent one standard deviation above and below that value.

Model Forecast Results and Associated Uncertainty    61 47°36' 47°42' 92°0' 91°54' 91°48' Base map data from , 2015 Minnesota Department of Natural Resource, 2011 4 MILES 4 KILOMETERS EXPLANATION Peter Mitchell pit at final buildout Lakes represented with the RIV package Probabilistic capture zone 0.125 to 0.25 0.25 to 0.375 0.375 to 0.5 0.5 to 0.625 0.625 to 0.75 0.75 to 0.875 0.875 to 1 Streams represented using the Streamflow Routing (SFR2) package Figure 37.  Probabilistic capture zones for the final iron mining scenario.

62    Groundwater/Surface-Water Interactions in the Partridge River Basin and Hypothetical Future Mine Pits The model prediction and estimated uncertainty for the 12 average depths to water in wetlands, 6 base-flow, and 5 pit inflow forecasts are presented in figures 34–36 and table 5. Each forecast shown includes uncertainty estimated using “linear” or First-Order, Second-Moment analysis (Doherty, 2015; Welter and others, 2015; Anderson and others, Chapter 10, 2015; White and others, 2016). The methods have the advantage of relatively low computational burden but assume that uncertainty extends equally around the forecasted value, which may not reflect realistic system responses. Forecast uncertainty was estimated using the python pack­ age, PyEMU (://github.com/​jtwhite79/​pyemu; White and others, 2016) that builds upon the linear uncertainty methods from Doherty (2015). The solid bars in the plots represent the modeled value of the forecast and the error bars represent one standard devia­ tion above and below that value. Error bars extending below zero were terminated at zero in the plots to censor unrealis­ tic results that stem from the First-Order, Second-Moment assumptions. The scenario number corresponds to the mine pit numbers shown in figure 33. The “final iron mining” scenario was the scenario with just the flooded Peter Mitchell pit and represented flooded pit conditions at the final mine buildout for this pit. The “base model” was the calibrated 2011–2013 mining conditions model. The estimated groundwater flux into the Peter Mitchell pit (fig. 34) was similar for all scenarios. The hypothetical mine pits in scenarios 1–4 (fig. 35) have similar predicted groundwater inflow rates with a lower rate in the scenario 2 pit; however, the difference was within the error. The streamflow forecasts (fig. 36) reflect direct ground­ water flow discharging into the stream cells and routed UZF seepage and rejected recharge under base-flow conditions. For the base model, the streamflow also included the 2.3 cubic feet per second pumped mine pit water discharged to the head­ waters of the Partridge River but was zero in the mine sce­ narios. To distinguish between the mine pit water and the other streamflow contributions (groundwater discharge and routed UZF water), the pumped pit water was plotted separately. The streamflow forecast for scenario 3 at location 4 was zero, rep­ resenting the case where the stream was forecasted to go dry. Except for forecast location 2, all mine pit scenarios caused reduced stream base flow relative to the base model. The big­ gest change to base flow occured from the Peter Mitchell pit at final buildout. This change is seen by comparing the base model and the final iron mining scenario; adding the hypothet­ ical pits in scenarios 1–4 resulted in more small reductions in streamflow relative to the final iron mining scenario. However, there is a noticeable decline in streamflow at forecast location 2 under scenario 1, which represented upstream potential min­ ing near this tributary. The forecasts for the average depth to water at wetlands (table 5) were calculated as the change in average depth to water from the model top elevation (relative to the base model) for all model cells intersected by the wetland areas shown in figure 33. Table 5 also presents the standard devia­ tion in water level changes across all model cells within each wetland forecast area and the estimated uncertainty in the base model averages. In general, the forecast uncertainty for the base model was lowest for wetlands where the average depth to water was less than one-half of the UZF SURFDEP term of 2.5 ft. At depths less than half SURFDEP, UZF is controlling the groundwater elevation by generating groundwater runoff and limiting the presence of pressurized areas of the model with hydraulic heads simulated well above the land surface. The sensitivity of changing model parameters (calculated in the Jacobian) on the modeled values at observations is likely small if these changes do not move the water table deep enough to fall below one-half of SURFDEP. Because there are no UZF parameters in the calibration, the model sensitivity to how UZF is affecting groundwater elevations is probably poorly captured in these wetland uncertainties. In addition to linear uncertainty in the base model forecasts, table 5 presents the average and standard deviation of the change in water depths across the forecast wetlands, relative to the base model. Table 4.  Summary of scenario names and features for the hypothetical mine pit scenarios. [lidar, light detecting and ranging; ft, feet] Scenario name Configuration of the Peter Mitchell pit Hypothetical mine pits Base model Pumped and excavated to level shown in the 2011 lidar. None. Final iron mining At final buildout topography and flooded to 1,500 ft, as dis­ cussed in Barr (2010). None. Scenario 1 At final buildout topography and flooded to 1,500 ft, as dis­ cussed in Barr (2010). Pit 1; fig. 33. Scenario 2 At final buildout topography and flooded to 1,500 ft, as dis­ cussed in Barr (2010). Pit 2; fig. 33. Scenario 3 At final buildout topography and flooded to 1,500 ft, as dis­ cussed in Barr (2010). Pit 3; fig. 33. Scenario 4 At final buildout topography and flooded to 1,500 ft, as dis­ cussed in Barr (2010). Pit 4; fig. 33.

Probabilistic Capture Zones    63 Wetland depths were similar across most of the model sce­ narios, with the biggest changes for forecast location 9 under scenario 4 and forecast 3 under all mining scenarios because it is close to the flooded Peter Mitchell pit. In general, as expected, wetlands closer to mine pits showed greater changes in average water depths. Probabilistic Capture Zones A capture zone represents the area contributing recharge to any hydrologic sink such as a well, spring, or a mine pit (Franke and others, 1999) and can be calculated by forward tracking particles from the water table (for example, Hunt and Steuer, 2000; Juckem and others, 2014). By using Monte Carlo techniques, it is possible to generate estimated capture zones that include the probability that recharge originating in each model cell is captured by a given model feature. Monte Carlo techniques use randomly generated model parameter sets that are sampled from a plausible parameter range to create numer­ ous model runs. The resulting capture zone generated by each model run can be tabulated, and a probability of capture can be calculated by dividing the total number of times a particle from a model cell was captured by the total number of realiza­ tions tried. This results in a capture probability ranging from 0 to 1 for each model cell. MODPATH, version 7 (Pollock, 2016), was used to do forward particle tracking with a particle placed in every model cell and assessing where each particle terminated. Appendix 4 provides additional discussion on the MODPATH and Monte Carlo setups. Probabilistic capture zones (figs. 37 and 38) were esti­ mated for each of the mine pits in scenarios 1–4 and for the final iron mining scenario of the flooded Peter Mitchell pit. The probabilistic capture zone had a value of one (captured in 100 percent of the Monte Carlo scenarios) for the immediate area around the Peter Mitchell pit (fig. 37). The probabilistic capture zones for the new mine pits in scenarios 1–4 (fig. 38) were generally confined to areas closer to the pit than the capture zones for the flooded Peter Mitchell pit, despite the new mine pits having lower water levels then the flooded Peter Mitchell pit. Scenario 1–4 pits were smaller and in deep crystalline bedrock with a low hydraulic conductivity relative to the upper layers of Biwabik Iron-Formation that surrounds the Peter Mitchell pit and likely have more prevalent frac­ turing. Pits in scenarios 3 and 4, immediately to the south of the flooded Peter Mitchell pit, had a line of particles that originated in the ridge formed by the Giants Range Granite to the north (fig. 38). These particle paths represented flow that originated as recharge in a topographic high, traveled into the deeper bedrock below the Peter Mitchell pit, and discharged into mine scenarios 3 and 4 pits. The probabilistic capture zone analysis in figures 37 and 38 also showed some low probability locations of capture relatively close to the mine pit Table 5.  Model forecast values for the wetland forecasts and associated uncertainty in the base model. A column is then included for each mine scenario with an average change in depth to water in the wetland, relative to the base model, and standard deviation of the changes for each model cell representing the wetland forecast. Negative values represent a lower water table for a wetland, relative to the base model, and a positive value represents a higher water table in the mine scenario than the base model. Wetland forecast number Average depth to water below model top in base model, in feet Estimated uncertainty for average depth to water below model top in base model, in feet Average change in final iron mining scenario wetland depth to water, relative to base model, with standard deviation, in feet Average change in scenario 1 wetland depth to water, relative to base model, with standard deviation, in feet Average change in scenario 2 wetland depth to water, relative to base model, with standard deviation, in feet Average change in scenario 3 wetland depth to water, relative to base model, with standard deviation, in feet Average change in scenario 4 wetland depth to water, relative to base model, with standard deviation, in feet 0.00 ± 0.00 0.00 ± 0.00 0.00 ± 0.00 0.00 ± 0.00 0.00 ± 0.00 0.00 ± 0.00 0.00 ± 0.00 0.00 ± 0.00 0.00 ± 0.00 0.00 ± 0.00 -40.87 ± 8.29 -40.87 ± 8.29 -40.88 ± 8.29 -49.68 ± 6.08 -42.04 ± 8.2 -7.47 ± 6.26 -7.47 ± 6.26 -7.84 ± 6.4 -8.11 ± 6.21 -7.47 ± 6.26 -0.05 ± 0.03 -0.05 ± 0.03 -0.05 ± 0.03 -0.06 ± 0.03 -0.05 ± 0.03 -0.03 ± 0.01 -0.03 ± 0.01 -0.03 ± 0.01 -0.03 ± 0.01 -0.03 ± 0.01 0.14 ± 0.04 0.14 ± 0.04 0.14 ± 0.04 0.14 ± 0.04 0.13 ± 0.04 -0.03 ± 0.04 -0.03 ± 0.04 -0.03 ± 0.04 -0.03 ± 0.04 -0.03 ± 0.04 -0.01 ± 0.01 -0.01 ± 0.01 -0.01 ± 0.01 -0.01 ± 0.01 -19.03 ± 10.83 0.00 ± 0.00 0.00 ± 0.00 0.00 ± 0.00 0.00 ± 0.00 -0.03 ± 0.06 0.00 ± 0.00 0.00 ± 0.00 0.00 ± 0.00 0.00 ± 0.00 0.00 ± 0.00 0.00 ± 0.01 0.00 ± 0.01 0.00 ± 0.01 0.00 ± 0.01 0.00 ± 0.01

64    Groundwater/Surface-Water Interactions in the Partridge River Basin and Hypothetical Future Mine Pits scenario stress, which may be a spurious result arising from issues with MODPATH’s handling of dry cells reported by the MODFLOW-NWT groundwater code (USGS, 2014). Data Worth Data worth analysis estimates how potential new obser­ vations can reduce uncertainty in a forecast; here the linear uncertainty analysis framework was extended to perform data worth analysis following the methods of Doherty (2015) and Fienen and others (2010). In this approach, data worth is the normalized difference between the forecast uncertainty with and without the new hypothetical data and represents how new data may improve the calibration of the base model such that the forecast values in the model scenario may become more certain. Appendix 4 provides additional information on how data worth was set up for the PRB model scenarios. Data worth results for most of the forecasts generally indicated a scattered pattern with relatively low data worth values across the model domain (fig. 39). This contrasts with the more classic data worth patterns in Leaf and others (2015), where data worth was generally highest in large, continu­ ous areas near the forecast location. Comparing the results between this study where data worth showed a more diffuse, scattered result against the more classic continuous pattern shown in Leaf and others (2015) indicates that either numeri­ cal instabilities or how the wetland areas were simulated may have complicated the application of data worth to this model. Most areas that showed cohesive patterns of data worth were best explained as wetland areas (peat deposits) with no UZF seepage. Areas with UZF seepage (fig. 32; 33 percent percent of the model domain for the 2011–2013 mining conditions model) showed very little data worth for any forecast in any scenario. The UZF package in MODFLOW limits the possibil­ ity for groundwater elevations above the model top. This limit is realistic in that the land surface does in fact act as a bound­ ary, but applying UZF this way also limits sensitivity esti­ mates that are calculated as the connection between changing parameters and a detectable numerical response in the ground­ water elevation. Unless the sensitivities from small parameter adjustments were calculated from a model solution where the modeled hydraulic head was below the top of the model, UZF would have just continued to pin the groundwater elevation at the land surface and any observations at or near UZF seepage areas would likely show very little change in modeled ground­ water elevation. This means that any new data in areas with or near strong UZF seepage would show little data worth because new information in such areas cannot be expected to inform parameter values or forecast estimates made with the model. UZF provided a tool to model groundwater flow in this poorly drained environment where wetlands in very thin unconsolidated materials overlie fairly low conductivity crystalline bedrock. An unforeseen challenge in this model­ ing approach was the complication of applying a data worth analysis to provide meaningful insights on future data collec­ tions because the model effectively fixed groundwater eleva­ tions in many of the wetland areas. In reality, with wetlands so prevalent across the landscape, the water table is likely known (within a few feet of land surface) in a large part of the PRB and the collection of new data in these wetland areas may not provide the model much new information that would reduce uncertainty. A better approach to planning future data collec­ tion may be to use estimated capture zones, as discussed in this study, or by simply applying hydrologic first principles to select monitoring locations relative to the locations of existing and new stresses on the system. Assumptions and Limitations Several simplifying assumptions were made that lim­ ited the accuracy of model results. The construction and subsequent utility of any model depends on the purpose of the model. For the PRB model described here, the purpose focused on simulations of regional groundwater flow across the PRB. This purpose, and the lack of site-specific data, limits model utility to basin-wide evaluations. As such, it is best suited to understanding general groundwater flow and ground­ water/surface-water interaction patterns across the basin at base-flow conditions. Similarly, all groundwater models are a necessarily simplified representation of a vastly compli­ cated natural flow system. For the model presented here, the simplifications used and incomplete knowledge of the sys­ tem resulted in several assumptions and limitations that are important for future use of the model. These model assump­ tions were made using available data and sound application of hydrologic principles based on modeling experience where data on the system were unavailable. Some key assumptions of the current PRB model are described in this section. A groundwater-flow system characterized by fractured, crystalline bedrock aquifer, thin or absent overlying unconsoli­ dated materials, numerous wetlands, and mining features with dramatic elevation changes is difficult to simulate. The mag­ nitude and direction of groundwater flow in bedrock is likely affected by the prevalence of fractures and their connectivity. Fracture location and connectivity is poorly understood in the model area, although previous studies and anecdotal evidence from shallow bedrock wells indicate that fracture flow could be appreciable. The fracture system was assumed to have decreasing hydraulic conductivity with depth. In addition, piecewise-constant hydraulic conductivity in the deeper bed­ rock units is a simplification of a more complex flow system. The assumption of piecewise-constant hydraulic conductivity in the bedrock and that the bedrock behaves as an equivalent porous medium may be reasonable at a basin-scale but this unsimulated heterogeneity in the bedrock may be important for understanding local groundwater flow at a local scale (for example, at a mine pit). Lastly, sloping bedrock units were created from a map of the bedrock surface units with a single

Assumptions and Limitations    65 Base map data from , 2015 Minnesota Department of Natural Resources, 2011 47°30' 47°36' 92°6' 92°0' 91°54' A 4 MILES 4 KILOMETERS EXPLANATION Lakes represented with the RIV package Probabilistic capture zone 0.125 - <0.25 0.25 - <0.375 0.375 - <0.5 0.5 - <0.625 0.625 - <0.75 0.75 - <0.875 0.875 - Scenario 1 hypothetical mine pit Peter Mitchell pit at final buildout Streams represented using the Streamflow Routing (SFR2) package Figure 38.  Probabilistic capture zones for each of the four hypothetical mine scenarios.

66    Groundwater/Surface-Water Interactions in the Partridge River Basin and Hypothetical Future Mine Pits Base map data from , 2015 47°30' 47°36' 92°6' 92°0' 91°54' B 4 MILES EXPLANATION Lakes represented with the RIV package Probabilistic capture zone 0.125 - <0.25 0.25 - <0.375 0.375 - <0.5 0.5 - <0.625 0.625 - <0.75 0.75 - <0.875 0.875 - Scenario 2 hypothetical mine pit Peter Mitchell pit at final buildout Streams represented using the Streamflow Routing (SFR2) package Figure 38.—Probabilistic capture zones for each of the four hypothetical mine scenarios.

Assumptions and Limitations    67 Base map data from , 2015 Mi t D t t f N t l R 47°42' 47°36' 92°6' 92°0' 91°54' 4 MILES EXPLANATION Lakes represented with the RIV package Probabilistic capture zone 0.125 - <0.25 0.25 - <0.375 0.375 - <0.5 0.5 - <0.625 0.625 - <0.75 0.75 - <0.875 0.875 - Scenario 3 hypothetical mine pit Peter Mitchell pit at final buildout Streams represented using the Streamflow Routing (SFR2) package Figure 38.—Probabilistic capture zones for each of the four hypothetical mine scenarios.

68    Groundwater/Surface-Water Interactions in the Partridge River Basin and Hypothetical Future Mine Pits Base map data from , 2015 Mi t D t t f N t l R 47°42' 47°36' 92°6' 92°0' 91°54' 91°48' D 4 MILES EXPLANATION Lakes represented with the RIV package Probabilistic capture zone 0.125 - <0.25 0.25 - <0.375 0.375 - <0.5 0.5 - <0.625 0.625 - <0.75 0.75 - <0.875 0.875 - Scenario 4 hypothetical mine pit Peter Mitchell pit at final buildout Streams represented using the Streamflow Routing (SFR2) package Figure 38.—Probabilistic capture zones for each of the four hypothetical mine scenarios.

Summary     69 assumed slope and manual adjustment to boundaries using borehole data. Additional three-dimensional information regarding bedrock stratigraphy would improve model repre­ sentation of these units. Surface infiltration of precipitation was the dominant source of water to the groundwater system in the PRB. A key assumption was that the pattern of net infiltration calculated by SWB was representative of actual groundwater recharge after calibration multipliers were applied to account for bias in the SWB results. In addition, the PRB conceptual model assumed that shallow flow through a poorly drained land surface was important for simulating stream base flow. The use of UZF reduces unrealistic overpressurization of the groundwater system in low-lying areas by conveying shallow flows to streams as an alternative to only allowing deeper flowpaths that discharge directly into streams. The relative importance of Dunnian overland flow, groundwater seeps, and direct dis­ charge of groundwater to stream base flow is not well under­ stood in systems with large wetland areas. With no calibration data to support or reject the existence of dominant shallow flow, choosing to model the system with UZF and where water is and is not routed are major assumptions. It was assumed that the calibrated Iron Range model simulates flow realistically enough in the PRB to extract rea­ sonable perimeter boundary conditions. Perimeter boundary conditions were set sufficiently distant to explicitly simulate competing hydrologic sinks outside the PRB when available, such as the Dunka, St. Louis, and Embarrass Rivers. However, perimeter boundary conditions may have appre­ ciable importance when the location of potential future stress does not have an explicitly simulated competing sink between the stress and the perimeter boundary conditions. It was also assumed that steady-state perimeter boundary conditions were sufficiently representative across the range of scenarios tested. Although the model was calibrated with available hydraulic head and streamflow data, the number, spatial cover­ age, and quality of the data was not uniform. Hydraulic head observations were concentrated near mining features and north of the Iron Range. Most of the PRB, particularly the central and southern parts, had little hydraulic head data coverage. Of the available hydraulic head data, a lack of vertical gradient information, such as data from piezometers, limited the abil­ ity to represent flow between the unconsolidated aquifer and bedrock aquifer, and the shallow and deep bedrock aquifers. Furthermore, much of the streamflow data came from a single synoptic flow survey. Future modeling efforts would benefit from continuous streamflow data in unmanaged parts of the basin. Simulation of existing mine features was necessarily simplistic. The PRB model assumed that 2011 lidar eleva­ tions adequately represented steady-state mining conditions. Given the dynamic nature of mining, 2011 conditions can only generally represent an average, modern mining condition. The actual depths of pits that were flooded in the 2011 lidar data were unknown and thus assumed. Although appropriate for regionally focused simulations, the tailings basin in the model domain was represented in a rudimentary way—a rela­ tively low hydraulic conductivity unit at land surface with no water being moved as part of the tailings processing. The use of estimated catchments for the mine pits to simulate runoff contributions to the pit water levels was an approxima­ tion and required assumptions of the drainage divides and a runoff coefficient. Finally, the mine conditions simulated by mining sce­ narios can only be simplistically represented because they simulate hypothetical pit configurations in identified mineral deposits. As such, the mining scenario simulations should be considered a starting point for evaluation of future mining, a demonstration of modeling methods, and a tool to assist baseline data collection before activities begin. Because the model purpose drives decisions regarding appropriateness of assumptions used, simulation of mine effects for changes to current mine operations or a new mining feature would require augmentation, refinement, and recalibration of the current PRB model with the new model purpose in mind. Summary The Partridge River Basin (PRB) is a mineral-rich area with historic iron mining and potential for new mines in copper-nickel deposits. The basin is characterized by extensive wetlands, lakes, and streams in poorly drained and often thin glacial material overlying Proterozoic bedrock. To better understand the interaction between these extensive surface water features and the groundwater system, a threedimensional, steady-state, groundwater-flow model of the PRB was developed by the U.S. Geological Survey in cooperation with the Great Lakes Indian Fish & Wildlife Commission using the finite-difference computer code MODFLOW-NWT. The purpose of this modeling effort was to quantify exist­ ing regional exchanges between the groundwater and surface water system and to evaluate the magnitude of potential changes to these exchanges from additional mining develop­ ment in the basin. New mining development included the development of four known mineral deposits in the PRB and the existing Peter Mitchell pit in a final buildout configuration and flooded. The primary outcomes from this modeling effort include the following. Two representations of the PRB system were devel­ oped: a long-term 1995–2015 average conditions model and a 2011–2013 mining conditions model. The models used similar properties based on avail­ able published information on the bedrock and glacial geology, hydrologic features, mining conditions, and wetland distribution. The models were calibrated to best available hydraulic head and streamflow measure­ ments. The 2011–2013 model better represented mine site conditions and was the primary tool used for

70    Groundwater/Surface-Water Interactions in the Partridge River Basin and Hypothetical Future Mine Pits 47°30' 47°36' 47°42' 92°12' 92°6' 92°0' 91°54' 91°48' Hillshade base from Minnesota Department of Natural Resources, 2011 8 MILES 8 KILOMETERS EXPLANATION Data worth for scenario 3 pit inflow, unitless 0.05 to 0.1 0.1 to 0.15 0.15 to 0.2 0.2 to 0.3 0.3 to 0.4 Groundwater seepage Peat deposits Lakes Partridge River Basin Streams represented using the Streamflow Routing (SFR2) package Model domain Hypothetical mine pit for scenario 3 Figure 39.  Data worth for predicting pit inflow to the scenario 3 pit for potential new hydraulic head observations in model layer 1.

Summary     71 evaluating regional groundwater flows and interactions between surface-water features and the groundwater system. Calibration of the PRB model resulted in ranges of glacial horizontal hydraulic conductivity parameters that generally agreed with literature values and other models of the region. Horizontal hydraulic conductiv­ ity of the bedrock was higher in the upper bedrock layers where numerous and continuous fractures have been observed and lower in the deeper bedrock layers. This report documents models that aimed to realisti­ cally represent the source of water to stream base flow. Use of MODFLOW’s UZF package allowed for explicit simulations of the shallow and deep groundwater flow system in the wetland-rich parts of the system where Quaternary unconsolidated deposits are thin or absent, drainage is poor, and the bedrock aquifer is fractured. The rejection of recharge and allowance of groundwater seepage resulted in increased importance of shorter, shallower flow paths that provided approximately 50 percent of the base flow in streams. The groundwater-flow system was best characterized where hydraulic head targets were concentrated near iron mining features in the northwestern half of the Partridge River Basin. The southeastern half of the model domain lacks calibration data, and simulation uncertainty in these areas can be expected to be higher. The 2011–2013 mining conditions model was then used to develop five model scenarios simulating the response of the groundwater and surface-water system to potential hydrologic stress. The purpose of these mine pit scenarios is to present a possible workflow to quantify a model’s uncertainty for a given model forecast and serve as a possible guide for initial data collection that may improve a future model’s ability to make such a forecast. The five model scenarios were used to forecast streamflow at six locations in the PRB, pit inflow rates for the new mine pits and the flooded Peter Mitchell pit, and the average depth to water in 12 wetlands. An uncertainty analysis was performed on streamflow, pit inflow, and average wetland depth-to-water fore­ casts in the five mine scenarios. Streamflows generally were reduced with future mining, and the greatest streamflow reductions occurred from the flooded Peter Mitchell Pit, probably due to its large size. Forecast streamflow changes were greatest in the headwaters and more minimal at downstream locations. Average depth to groundwater in wetlands was most affected the closer the wetland was to a new mine pit. Probabilistic capture zones were estimated for the five mining scenarios using Monte Carlo techniques. Generally, pit capture was dominantly expressed in areas adjacent to the pit. However, capture zones were not always continuous, particularly for pits south of the flooded Peter Mitchell pit where pit capture is expected to occur near the proposed pit and north of the flooded Peter Mitchell pit with deeper groundwater flowpaths that originate on a topographic high, move under the flooded pit, and discharge into the proposed pit. The PRB simulations indicated the possibility of deeper flow through the bedrock. Data-worth analysis of potential, new hydraulic head observations was complicated by using UZF to treat wetlands as a boundary condition and by pos­ sible instabilities in the model solution. These results highlight the dominance of groundwater runoff and wetland prevalence in the PRB model and the lack of sensitivity of new data in areas where groundwater seepage is occurring. In reality, collecting new data in these wetland areas would likely show a water level at or near land surface anyway and may not constrain the model solution much more than the assumptions made by allowing groundwater runoff. New data collection may be better guided by the capture zone analysis or by applying hydrologic first principles to design a monitoring network. This model was developed to focus on regional flow across the PRB, thus assumptions used were evaluated with regional flow as the primary concern. Given the regional focus of the model, site-specific applications of the model to areas of a local scale within the basin would likely require additional site-specific data col­ lection and model recalibration. The hypothetical mining scenarios were developed to inform the initial mine permitting process of poten­ tial new mining projects in the basin. The hypo­ thetical mining scenarios were purposely simplified and not at all predictive of conditions for a specific future mine. The capture zones may be a useful first assessment of where to consider baseline monitor­ ing locations. The workflow described in this report provides a template for assessing potential new mining projects and integrating additional information as it becomes available.

72    Groundwater/Surface-Water Interactions in the Partridge River Basin and Hypothetical Future Mine Pits References Cited Adams, J.L., Leibfried, R.T., and Herr, E.S., 2004, East range hydrology project: Minnesota Department of Natural Resources Division of Lands and Minerals: Minnesota Department of Natural Resources Division of Land and Minerals Report, 102 p. Alderdice, L., Adams, J.L., and Eger, P.A., 1992, Managing the hydrologic impacts of mining on Minnesota’s Mesabi Iron Range: Proceedings America Society of Mining and Reclamation, accessed April 24, 2019, at ://www.asmr.us/​Portals/​0/​Documents/​Conference-​ Proceedings/​1992/​0108-​Alderdice.pdf. Anderson, M.P., Hunt, R.J., Krohelski, J.T., and Chung, K., 2002, Using high hydraulic conductivity nodes to simulate seepage lakes: Ground Water, v. 40, no. 2, p. 119–124. Anderson, M.P., Woessner, W.W., and Hunt, R.J., 2015, Applied Groundwater Modeling—Simulation of Flow and Advective Transport (2d ed): Online only, Academic Press, Inc., 564 p. Barr Engineering Company [Barr], 2006a, Drill hole monitor­ ing and data collection, hydrogeological investigation— Phase 1, PolyMet NorthMet Mine Site: Prepared for PolyMet Mining Inc., RS–02, variously paged, accessed November 16, 2017, at ://fi​les.dnr.st​ate.mn.us/​lands_​ minerals/​northmet/​water-​approp/​references/​mine-​site-​ hydro-​phase1-​rs-​02.pdf. Barr Engineering Company [Barr], 2006b, Drill hole monitor­ ing and data collection, hydrogeological investigation— Phase II, PolyMet NorthMet Mine Site: Prepared for PolyMet Mining Inc., RS–10, variously paged, accessed December 14, 2017, at ://fi​les.dnr.st​ate.mn.us/​lands_​ minerals/​northmet/​water-​approp/​references/​mine-​site-​ hydro-​phase2-​rs-​02.pdf. Barr Engineering Company [Barr], 2008, Long-range hydrol­ ogy study: Prepared for Northshore Mining Company, accessed April 28, 2020 at ://​www.leg.st​ate.mn.us/​ docs/​2015/​other/​150681/​PFEISref_​1/​Barr%202008o.pdf. Barr Engineering Company [Barr], 2010, Peter Mitchell Pit Concept Mitigation Plan Aquatic Habitat Enhancement through In-Pit Disposal of Overburden and Rock: Prepared for Northshore Mining Company, 49 p. accessed March 6, 2019, at ://​www.leg.st​ate.mn.us/​docs/​2015/​other/​ 150681/​PFEISref_​2/​MDNR%202011p.pdf. Barr Engineering Company [Barr], 2011, Comparison of hydrogeologic setting—Canisteo Pit, Minntac Mine near Kinney, and NorthMet Mine Site: Barr technical memo­ randum, 13 p., accessed November 16, 2017, at ://​ www.leg.st​ate.mn.us/​docs/​2015/​other/​150681/​PFEISref_​1/​ Barr%202011i.pdf. Barr Engineering Company [Barr], 2014a, Hydrogeology of the fractured bedrock in the vicinity of the NorthMet Project: Prepared for PolyMet Mining Inc., variously paged, accessed July 28, 2017, at ://fi​les.dnr.st​ate.mn.us/​lands_​ minerals/​northmet/​water-​approp/​references/​hydrogeology-​ fractured-​bedrock-​v3.pdf. Barr Engineering Company [Barr], 2014b, Groundwater modeling of the NorthMet Mine Site—Supporting docu­ ment for water modeling data package volume 1—Mine site: Prepared for PolyMet Mining Inc., variously paged, accessed December 14, 2017, at ​://gis.lic​.wisc.edu/​ wwwlicgf/​glifwc/​PolyMet/​FEIS/​reference/​PolyMet%20 2015m_​Attach_​B_​GW-​modeling.pdf. Bay, R.R., 1968. The hydrology of several peat deposits in northern Minnesota, U.S.A, in Proceedings of the Third International Peat Congress: Quebec, Canada, National Research Council of Canada, p. 212–218, accessed December 12, 2020, at ​s://www.fs​.usda.gov/​treesearch/​ pubs/​12896. Card, K.D., 1990, A review of the Superior Province of the Canadian Shield, a product of Archean accre­ tion: Precambrian Research, v. 48, np. 1–2. p. 99–156 [Also available at ://doi.org/​10.1016/​0301-​ 9268(90)90059-​Y]. Cotter, R.D., Young, H.L., Petri, L.R., and Prior, C.H., 1965, Ground and surface water in the Mesabi and Iron Range area, northeastern Minnesota: U.S. Geological Survey Water-Supply Paper 1759–A, 36 p., 1 pl. [Also available at ://doi.org/​10.3133/​wsp1759A.]. County Well Index, 2018, The Minnesota County Well Index: Minnesota Geological Survey, University of Minnesota web page, accessed December 11, 2018, at ://m​gsweb2.mng​ s.umn.edu/​cwi_​doc/​cwidoc.htm. Dengler, E., and Wagner, K., 2017, Preliminary surficial geologic map of southeastern Arrowhead area, Lake and St. Louis Counties, northeastern Minnesota: Minnesota Geological Survey Open-File Report OFR2016-04, 1 p., accessed April 15, 2019, at ://​conservanc​y.umn.edu/​ handle/​11299/​183258. Dengler, E.L., Wagner, K.G., and Hamilton, J., 2017, Preliminary Quaternary geology and sand distribution models of the central Arrowhead Area, Lake and St. Louis Counties, northeastern Minnesota: Minnesota Geological Survey Open-File Report OFR2016-04, 1 p., accessed April 15, 2019, at ://​conservanc​y.umn.edu/​handle/​ 11299/​183258.

References Cited    73 Doherty, J.E., and Hunt, R.J., 2010, Approaches to highly parameterized inversion—a guide to using PEST for groundwater-model calibration: U.S. Geological Survey Scientific Investigations Report 2010–5169, 59 p. accessed April 28, 2020, at ht​tps://pubs​.usgs.gov/​sir/​2010/​5169/​pdf/​ GWPEST_​sir2010-​5169.pdf. Doherty, J., 2015, Calibration and uncertainty analysis for complex environmental models: Brisbane, Australia, Watermark Numerical Computing, 227 p. Doherty, J., 2018a, PEST, Model-independent parameter estimation—user manual part II: PEST utility support soft­ ware (7th ed.): Brisbane, Australia, Watermark Numerical Computing, 267 p., accessed April 28, 2020, at ://w​ ww.pesthom​epage.org/​Downloads.php [Also available at ​://pesthom​epage.org/​documentation.]. Doherty, J., 2018b, PEST, Model-independent parameter estimation—user manual part I: PEST, SENSAN and Global Optimisers (7th ed.): Brisbane, Australia, Watermark Numerical Computing, 393 p., accessed April 28, 2020, at ://w​ww.pesthom​epage.org/​Downloads.php [Also avail­ able at ​://pesthom​epage.org/​documentation.]. Doherty, J., 2020, PEST for highly parallelized computing environments: Brisbane, Australia, Watermark Numerical Computing, 88 p., accessed April 28, 2020, at ://w​ ww.pesthom​epage.org/​Downloads.php. Duffield, G.M., 2019, Representative values of hydraulic prop­ erties: AQTESOLV web page, accessed June 4, 2020, at ​ ://www.aqt​esolv.com/​aquifer-​tests/​aquifer_​properties.htm. Fienen, M.N., Doherty, J.E., Hunt, R.J., and Reeves, H.W., 2010, Using prediction uncertainty analysis to design hydrologic monitoring networks—Example applications from the Great Lakes Water Availability Pilot Project: U.S. Geological Survey Scientific Investigations Report 2010–5159, 44 p., ht​tps://pubs​.usgs.gov/​sir/​2010/​5159/​. Feinstein, D.T., Hart, D.J., Gatzke, S., Hunt, R.J., Niswonger, R.G., and Fienen, M.N., 2019, A simple method for simulat­ ing interactions with fens to forecast development effects: Ground Water, v. 58, no. 4, p. 524–534. [Also available at ://doi.org/​10.1111/​gwat.12931]. Franke, O.L., Reilly, T.E., Pollock, D.W., and LaBaugh, J.W., 1999, Estimating areas contributing recharge to wells, lessons from previous studies: U.S. Geological Survey Circular 1174, 20 p., ://doi.org/​10.3133/​cir1174. Harbaugh, A.W., Banta, E.R., Hill, M.C., and McDonald, M.G., 2000, MODFLOW-2000, the U.S. Geological Survey modular ground-water model—User guide to modulariza­ tion concepts and the ground-water flow process: U.S. Geological Survey Open-File Report 00-92, 121 p. Haserodt, M.J., Hunt, R.J., Cowdery, T.K., Leaf, A.T., and Baker, A.C., 2019, Simulation of the regional groundwa­ ter flow system in the St. Louis River Basin, Minnesota: U.S. Geological Survey Scientific Investigations Report 2019–5033, 41 p., ://doi.org/​10.3133/​sir20195033. Haserodt, M.J., Hunt, R.J., Fienen, M.N., and Feinstein, D.J., 2021, MODFLOW-NWT and MODPATH models, capture zones and uncertainty data analysis for the Partridge River Basin, Minnesota, U.S. Geological Survey data release, ://doi.org/​10.5066/P9VODOU8. Hoyt Lakes, 2016, Hoyt Lakes 2016 Drinking Water Report: Hoyt Lakes Consumer Confidence Report, accessed October 11, 2019, at h​ttp://hoyt​lakes.com/​2016%​ 20Consumer​%20Confide​nce%20Repo​rt%20.pdf/​. Hunt, R.J., and Steuer, J.J., 2000, Simulation of the recharge area for Frederick Springs, Dane County, Wisconsin: USGS Water-Resources Investigations Report 00-4172, 33 p. Hunt, R.J., Haitjema, H.M., Krohelski, J.T., and Feinstein, D.T., 2003, Simulating ground water-lake interactions— Approaches and insights: Ground Water, v. 41, no. 2, p. 227–237. Hunt, R.J., Doherty, J., and Tonkin, M.J., 2007, Are models too simple? Arguments for increased parameterization: Ground Water, v. 45, no. 3, p. 254–262. Hunt, R.J., Prudic, D.E., Walker, J.F., and Anderson, M.P., 2008, Importance of unsaturated zone flow for simulating recharge in a humid climate: Ground Water, v. 46, no. 4, p. 551–560. 10.1111/j.1745-6584.2007.00427.x. Hunt, R.J., and Feinstein, D.T., 2012, MODFLOW-NWT – Robust handling of dry cells using a Newton Formulation of MODFLOW-2005: Ground Water, v. 50, no. 5, p. 659–663. 10.1111/j.1745-6584.2012.00976.x. Jirsa, M.A., Chandler, V.W., and Lively, R.S., 2005, Bedrock geology of the Mesabi Iron Range, Minnesota: Minnesota Geological Survey Miscellaneous Map Series Map M-163, accessed September 18, 2020, at ://​conservanc​ y.umn.edu/​handle/​11299/​58043. Jirsa, M.A., Boerboom, T.J., Chandler, V.W., Mossler, J.H., Runkel, A.C., and Setterholm, D.R., 2011, Geologic map of Minnesota—Bedrock geology: Minnesota Geological Survey State Map Series S–21, 1 sheet plus supplemental files, scale 1:500,000, accessed July 31, 2017, at ://​ conservanc​y.umn.edu/​handle/​11299/​10146­6. Jirsa, M.A., 2016a, Preliminary geologic maps of Lake and St. Louis Counties, northeastern Minnesota: Minnesota Geological Survey OFR 16-4, 1 p., accessed October 11, 2019 at ://​conservanc​y.umn.edu/​handle/​11299/​183258.

74    Groundwater/Surface-Water Interactions in the Partridge River Basin and Hypothetical Future Mine Pits Jirsa, M.A., 2016b, Preliminary geologic maps of Lake and St. Louis Counties, northeastern Minnesota—Draft merged surficial geology: Minnesota Geological Survey OFR 16-4, 1 p., accessed March 27, 2019, at ht​tp://hdl.h​andle.net/​ 11299/​183258. Jones, P.M., 2002, Characterization of ground-water flow between the Canisteo Mine Pit and surrounding aquifers, Mesabi Iron Range, Minnesota: U.S. Geological Survey Water-Resources Investigations Report 02–4198, 34 p. [Also available at ://doi.org/​10.3133/​wri024198.]. Juckem, P.F., Fienen, M.F., and Hunt, R.J., 2014, Simulation of groundwater flow and interaction of groundwater and surface water on the Lac du Flambeau Reservation, Wisconsin: U.S. Geological Survey Scientific Investigations Report 2014–5020, 34 p., ://dx.doi.org/​10.3133/​ sir20145020. Konikow, L.F., Hornberger, G.Z., Halford, K.J., and Hanson, R.T., 2009, Revised multi-node well (MNW2) package for MODFLOW ground-water flow model: U.S. Geological Survey Techniques and Methods 6–A30, 67 p. Leaf, A.T., Fienen, M.N., Hunt, R.J., and Buchwald, C.A., 2015, Groundwater/surface-water interactions in the Bad River Watershed, Wisconsin: U.S. Geological Survey Scientific Investigations Report 2015–5162, 110 p., ://dx.doi.org/​10.3133/​sir20155162. Leaf, A.T., 2018, SFRmaker: GitHub software web page, accessed November, 27, 2018, from ://github.com/​ aleaf/​SFRmaker. Lindholm, G.F., Ericson, D.W., Broussard, W.L., and Hult, M.F., 1979, Water resources of the St. Louis River watershed, northeastern Minnesota: U.S. Geological Survey Hydrologic Investigations Atlas HA–586, 3 sheets, scale 1:250,000. [Also available at ://doi.org/​ 10.3133/​ ha586.]. McKay, L., Bondelid, T., Dewald, T., Johnston, J., Moore, R., and Rea, A., 2012, NHDPlus Version 2—User guide: NHDPlus web page, accessed February 4, 2019, at ://www.horizon-​systems.com/​NHDPlus/​ NHDPlusV2_​home.php. Miller, J.D., Green, J.C., Severson, M.J., Chandler, S.A., and Peterson, D.M., 2001, Geologic map of the Duluth Complex and related rocks, northeastern Minnesota: Minnesota Geological Survey Miscellaneous Map Series Map M-119, 1 p., accessed October 11, 2019, at ://​conservanc​ y.umn.edu/​handle/​11299/​183. Miller, J.D., Green, J.C., Severson, M.J., Chandler, V.W., Hauck, S.A., Peterson, D.M., and Wahl, T.E., 2002, Geology and mineral potential of the Duluth Complex and related rocks of Northeastern Minnesota: Minnesota Geological Survey Report of Investigations 58, 207p., accessed October 10, 2019, at ://​conservanc​y.umn.edu/​ bitstream/​handle/​11299/​58804/​ri_​58%5b1%5d.pdf?​ sequence=​3&isAllowed=​y. Minnesota Department of Natural Resources [MDNR], 2007, Lake Information Report, Name: Whitewater: Minnesota DNR Lake information report, accessed April 2019 at htt​ ps://www.l​eg.mn.gov/​docs/​2015/​other/​150681/​PFEISref_​2/​ MDNR%202007b.pdf. Minnesota Department of Natural Resources [MDNR], 2011, LiDAR elevation, Arrowhead Region, NE Minnesota, 2011: Minnesota IT Services Geospatial Information Office data­ set, accessed June 2017, at ://ww​w.mngeo.st​ate.mn.us/​ chouse/​elevation/​lidar.. Minnesota Department of Natural Resources [MDNR], United States Army Corps of Engineers [USACE], and United States Forest Service [USFS], 2013, Northmet mining project and land exchange supplemental draft environment impact statement: MDNR web page, accessed May 28, 2020, at ://​www.dnr.st​ate.mn.us/​input/​ environmentalreview/​polymet/​seis-​toc.. Minnesota Department of Natural Resources [MDNR], 2016, August 2016 – Exploration for metallic mineral resources in Minnesota Duluth Complex: MDNR web page, accessed April 2019 at ://​www.dnr.st​ate.mn.us/​lands_​minerals/​ metallic_​nf/​exp_​areas.. Minnesota Department of Natural Resources [MDNR], 2018, Partridge River nr Hoyt Lakes, 2.3mi us of CSAH110 (03149002) dataset: MDNR web page, accessed December 17, 2018, at ://​www.dnr.st​ate.mn.us/​waters/​ csg/​site_​report.?​mode=​getsitereport&site=​03149002. Minnesota Department of Natural Resources [MDNR], 2019a, DNR hydrography dataset: MDNR web page, accessed April 2019 at htt​ps://gisda​ta.mn.gov/​dataset/​water-​dnr-​ hydrography. Minnesota Department of Natural Resources [MDNR], 2019b, St. James Mine Pit: MDNR web page, accessed October 11, 2019, at ://​www.dnr.st​ate.mn.us/​lands_​minerals/​ waters_​program/​st-​james.. Minnesota Department of Natural Resources [MDNR], 2019c, Water use permits: MDNR web page, accessed October 11, 2019, at ://​www.dnr.st​ate.mn.us/​waters/​watermgmt_​ section/​appropriations/​permits..

References Cited    75 Minnesota Department of Natural Resources [MDNR], 2019d, Minnesota Department of Natural Resources Water Appropriations Permit Program data: MDNR web page, accessed April 30, 2019, at ://​www.dnr.st​ate.mn.us/​ waters/​watermgmt_​section/​appropriations/​wateruse.. Minnesota Geological Survey, 2013, Bedrock topography and depth to bedrock: Minnesota Geological Survey dataset, accessed March 30, 2017, at ​://www.mng​s.umn.edu/​ service.htm. Minnesota Mining, 2019, Minnesota’s vast mineral resources: Minnesota Mining web page, accessed April 2019 at h​ ttp://www.​miningminn​esota.com/​duluth-​complex/​. Myers, T., 2014, Technical memorandum: Twin Metals min­ ing and the Boundary Waters Canoe Area Wilderness, risk assessment for underground metals mining: Save the Boundary Waters web page, accessed October 15, 2019, at ://​www.saveth​eboundaryw​aters.org/​sites/​default/​files/​ attachments/​myers_​2014_​-​_​final_​flow_​transport_​analysis_​ tm080914_​-​_​risks_​to_​bwcaw.pdf. Niswonger, R.G., and Prudic, D.E., 2005, Documentation of the Streamflow-Routing (SFR2) Package to include unsatu­ rated flow beneath streams—A modification to SFR1: U.S. Geological Survey Techniques and Methods Book 6, chap. A13, 50 p. Niswonger, R.G., Prudic, D.E., and Regan, R.S., 2006, Documentation of the Unsaturated-Zone Flow (UZF1) Package for modeling unsaturated flow between the land surface and the water table with MODFLOW–2005: U.S. Geological Survey Techniques and Methods Book 6, chap. A19, 62 p. Niswonger, R.G., Panday, S., and Ibaraki, M., 2011, MODFLOW–NWT—A Newton formulation for MODFLOW–2005: U.S. Geological Survey Techniques and Methods Book 6, chap. A37, 44 p. Oakes, E.L., 1970, Geology and ground-water resources on the Grand Rapids area, north-central Minnesota: U.S. Geological Survey Hydrologic Investigations Atlas HA–322, 2 sheets, scale 1:48,000. Olcott, P.G., Ericson, D.W., Felsheim, P.E., and Broussard, W.L., 1978, Water resources of the Lake Superior water­ shed, northeastern Minnesota: U.S. Geological Survey Hydrologic Investigations Atlas HA–582, 2 sheets, scale 1:500,000. Pollock, D.W., 2016, User guide for MODPATH Version 7—A particle-tracking model for MODFLOW: U.S. Geological Survey Open-File Report 2016–1086, 35 p., ://dx.doi.org/​10.3133/​ofr20161086. Schwartz, F.W., and Zhang, H., 2003, Fundamentals of ground water: New York, John Wiley & Sons, 583 p. Siegel, D.I., and Ericson, D.W., 1980, Hydrology and water quality of the copper-nickel study region, northeastern Minnesota: U.S. Geological Survey Open-File Report 80–739, 87 p., 2 pl., accessed November 16, 2017, at ://​conservanc​y.umn.edu/​handle/​11299/​189248. Smith, E.A., 2017, Soil-water-balance model data sets for the St. Louis River drainage basin, northeast Minnesota, 1995–2010: U.S. Geological Survey data release, accessed June 2018 at ://doi.org/​10.5066/​F7Z60MJ0. Smith, E.A., and Westenbroek, S.M., 2015, Potential ground­ water recharge for the State of Minnesota using the Soil Water-Balance model, 1996–2010: U.S. Geological Survey Scientific Investigations Report 2015–5038, 85 p., accessed November 16, 2017, at ://doi.org/​10.3133/​sir20155038. St. Louis County, Minnesota, 2017, Parcels, Saint Louis County, MN: St. Louis County, Minnesota GIS & Mapping Portal geospatial dataset, accessed April 2019 at ://open-​data-​sl​cgis.hub.a​regis.com/.​ Stark, J.R., 1977, Surficial geology and ground-water geol­ ogy of the Babbitt- Kawishiwi area, northeastern Minnesota with planning implications: Madison, Wisconsin, University of Wisconsin-Madison, M.S. thesis, 127 p. Tetra Tech, 2014, Upper St. Louis River watershed mining area hydrology: Prepared by Tetra Tech for the Minnesota Pollution Control Agency, 88 p., accessed December 14, 2017, at ://​www.pca.st​ate.mn.us/​sites/​default/​files/​wq-​ iw10-​12p.pdf. Tieberg, A., 2020, A brief history of Minnesota’s Mesabi Iron Range: MINNPOST (Minnesota), May 26, 2020, accessed November 25, 2020, at ​://www.min​post.com/​ mnopedia/​2020/​05/​a-​brief-​history-​of-​minnesotas-​mesabi-​ iron-​range/​. U.S. Fish & Wildlife Service [USFWS], 2018, National wetlands inventory – Version 2 – Surface water and wet­ lands inventory: U.S. Fish and Wildlife geospatial dataset, accessed April 24, 2019, at ://ww​w.mngeo.st​ate.mn.us/​ chouse/​water_​wetlands.. U.S. Geological Survey, 2012, 1:24,000-scale hydrologic units of the United States: U.S. Geological Survey dataset, accessed January 25, 2012, h​://www​.usgs.gov/​core-​ science-​systems/​ngp/​national-​hydrography/​watershed-​ boundary-​dataset?​qt-​science_​support_​page_​related_​con=​ 4#qt-​science_​support_​page_​related_​con. U.S. Geological Survey, 2014, Compatibility of MODPATH with MODFLOW—NWT. U.S. Geological Survey readme, accessed August 24, 2020, at ://water.usgs.gov/​ water-resourses/​software/​MODPATH​/UsingModpath​With​ Modflow-NWT.txt​.

76    Groundwater/Surface-Water Interactions in the Partridge River Basin and Hypothetical Future Mine Pits U.S. Geological Survey, 2019, USGS water data for the Nation: U.S. Geological Survey National Water Information System database, accessed March, 2019, at ://doi.org/10.5066/F7P55KJN. Welter, D., White, J., Doherty, J., and Hunt, R., 2015. PEST++ version 3, a parameter estimation and uncertainty analysis software suite optimized for large environmental models: U.S. Geological Survey Techniques and Methods Report 7-C12, 54 p. White, J.T., Fienen, M.N., and Doherty, J.E., 2016, A python framework for environmental model uncertainty analysis: Environmental Modelling & Software, v. 85, p. 217–228.

Appendix 1.  Additional Data Processing Steps to Build the MODFLOW-NWT Packages    77 Appendix 1.  Additional Data Processing Steps to Build the MODFLOW-NWT Packages A finite-difference groundwater-flow model of the Partridge River Basin (PRB) was constructed using MODFLOW-NWT with the Newton Raphson solver Niswonger and others, 2011). Several MODFLOW packages were used to simulate the different components of the hydro­ logic system. This appendix provides supporting information on how the more complicated packages were formulated for this model including the following: logic used to develop the model layers for the discretization (DIS) package, processing of Soil-Water-Balance (SWB) net infiltration data, delineation of active and routed areas for the Unsaturated Zone Flow (UZF) package, the setup of the Stream Flow Routing (SFR2) package with mine water inputs, setup of the river (RIV) package to represent lakes, and how pumped and flooded mine pits were represented. Model Discretization (DIS Package): Layering Logic The discretization of the PRB model layers were decided using the same logic used for the T.K. Cowdery (USGS, writ­ ten commun., 2020) MODFLOW model and described in this section. The T.K. Cowdery (USGS, written commun., 2020) MODFLOW model used four bedrock layers; three additional bedrock layers (seven bedrock layers total) were added in the PRB model to accommodate pit depths for the future mining scenarios. Layer thickness was assigned using the follow­ ing logic. Unconsolidated Layers Layer 1 had a maximum thickness of 20 feet (ft) or was set to the total unconsolidated thickness if the model top elevation minus top of bedrock elevation was less than 20 ft. The layer 1 bottom was the model top elevation minus the assigned thickness. Layer 2 had a maximum thickness of 80 ft or was set to the total remaining unconsolidated thickness if the difference between the bottom elevation of layer 1 and the top of bedrock was less than 80 ft. Layer 2 was assigned a 1-ft thickness if the entire unconsolidated thickness was represented by layer 1. Layer 3 had a maximum thickness of 100 ft or was set to the total remaining unconsolidated thickness if the bottom elevation of layer 2 minus the top of bedrock elevation was less than 100 ft. Layer 3 was assigned a thickness of 1 ft for instances where the entire uncon­ solidated thickness was represented by layers 1 and 2. Layer 4 had a maximum thickness of 104 ft or was set to the total remaining unconsolidated thickness if the bottom elevation of layer 3 minus the top of bedrock elevation was less than 104 ft. If all the unconsolidated thickness was represented by layers 1, 2, and 3, then layer 4 was assigned a thickness of 1 ft. Bedrock Layers Layer 5 had a uniform thickness of 50 ft. The top elevation of layer five was the top of the bedrock elevation (Dengler and Wagner, 2017; Dengler and others, 2017). Remaining layers below layer 5 were set to the bottom elevation of the overlying layer minus the thickness assigned to that layer. Layer 6 was assigned a uniform thickness of 50 ft. Layer 7 was assigned a uniform thickness of 100 ft. Layers 8–11 were each assigned a uniform thickness of 200 ft. Processing of the SWB Net Infiltration Estimates of net infiltration, or water that goes past the root zone, came from a published SWB model for the St. Louis River Basin (SLRB; Smith, 2017; Cowdery and others, 2021b) and the Minnesota statewide SWB model for areas outside the SLRB (Smith and Westenbroek, 2015 and Cowdery and others, 2021b). The statewide SWB data were available at a 0.62-mile (1-kilometer) resolution and SLRB SWB data were available at a 0.062-mile (0.1-kilometer) resolution. The combined SWB array was resampled to an average net infiltration value at each model grid cell using the zonal_stats function from the rasterstats python package (Perry, 2015). Net infiltration values for the 1995–2015 condi­ tions model were the average of the 1995–2015 SWB yearly data; net infiltration for the mining conditions model are the average the 2011–2013 SWB yearly data. Net infiltration over

78    Groundwater/Surface-Water Interactions in the Partridge River Basin and Hypothetical Future Mine Pits lakes, which are modeled as river cells, was set to 0 inch per year. Net infiltration over flooded pits, pumped pits, and the tailings basin pond was set to the average net precipitation rate, defined as the average annual precipitation minus average annual potential evapotranspiration, as calculated in the SLRB SWB model input files. The net precipitation in the SLRB for 1995–2015 was 1.2 inches per year and for 2011–2013 was 1.4 inches per year. Delineation of Active and Routed Areas for the UZF Package The UZF package (Niswonger and others, 2006) was used to represent the processes occurring in areas with high water tables. UZF simulated rejection of recharge and ground­ water seepage (collectively considered groundwater runoff) to the land surface when the water table is near the land surface. Model areas where UZF was active and used to simulate groundwater runoff and recharge are “specified as active” in the IUZFBND. The entire PRB model domain used “specified as active” except for lakes, which were simulated as RIV cells, and flooded mine pits, where net precipitation was applied with surface-water runoff using an injection well. The rejected groundwater recharge and the groundwater dis­ charge from active areas were routed to nearby SFR2 stream cells. SFR2 segment numbers were assigned to the UZF IRUNBND array using the catchment areas from NHDPlus v2 (McKay and others, 2012). IRUNBND cells without an SFR segment number are not routed and represent closed depres­ sions that were identified in the land surface using geographic information system (GIS) utilities; groundwater runoff from these closed depressions was assumed to be evapotranspired and therefore was removed from the model. Additional areas not routed included areas near mining features, which were assumed to have altered drainage patterns, and areas where the landscape is hummocky and not likely to contribute to surface drainage. The smoothness of the land surface was determined using the roughness tool in QGIS (ht​tps://docs​.qgis.org/​2.8/​ en/​docs/​user_​manual/​processing_​algs/​gdalogr/​gdal_​analysis/​ roughness.) with the lidar digital elevation model (DEM) resampled to 131-ft (40-meter) cells. Areas with a rough­ ness value greater than 0.5 were not routed. All areas within 492 ft (150 meters) of an SFR2 segment, regardless of either proximity to mining features or roughness, were assumed to be routed. Catchments containing no SFR2 segments were not routed. Figure 1.1 shows the IRUNBND array of where groundwater runoff is being routed to surface water or, alterna­ tively, assumed to evapotranspire to the atmosphere. Updates to the SFR2 Package to Include Mine Inputs Streams were modeled as a head-dependent flux bound­ ary condition with the SFR2 package (Niswonger and Prudic, 2005). Information in NHDPlus version 2 (McKay and others, 2012) was used to create the SFR2 model files with the python package SFRmaker (Leaf, 2018; an updated and more acces­ sible version of SFRmaker can be found at Leaf and others, 2021). There are two mining features in the model domain that discharge to the stream system. These reported discharges were included in the SFR2 package as stream inflows to an SFR2 cell. The inflow rates were based on reported values. Pit overflow from a flooded mine pit to Wyman Creek was modeled as a stream input in the stream reach at the pit overflow location (white circle with black center shown on fig. 12). The reported overflow rates for 2011 ranged from less than 0.02 cubic foot per second (ft3/s) to greater than 2.2 ft3/s with an average flow of 1 ft3/s (Richard Clark, Minnesota Pollution Control Agency, written commun., 2019). Based on a downstream measured flow of 1 ft3/s during the U.S. Geological Survey synoptic flow survey in 2018 (NWIS, 2019) and the assumption of some base flow, pit overflow to Wyman Creek was specified in the model at 0.5 ft3/s. The pit transfer from the pumped pit at the headwaters of the Partridge River (shown as a single red WEL cell in fig. 12) was mod­ eled as inflow to the Partridge River SFR2 cells; the inflow rate was 2.3 ft3/s based on the total reported gallons pumped in 2011 (Avery Guertin, Minnesota Department of Natural Resources, written commun., 2019). Representation of Lakes with the RIV Package Lakes were modeled as a head-dependent flux bound­ ary condition using the RIV package (Harbaugh and others, 2000) with the river stage set to the lake-surface elevation. Lake extents were from the Minnesota Department of Natural Resources Hydrography Dataset (Minnesota Department of Natural Resources [MDNR], 2019); any “Lake or Pond” feature more than 20 acres was included in the model. Lake stages in the RIV cells were set as the minimum lidar elevation (MDNR, 2011) in a model cell occupied by a lake. The lake bottom was set at 1 ft below the stage with a 1-ft thick lakebed. The conductance of the lakebed depends on the cell area, the thickness of the lakebed, and the vertical hydraulic conductivity of the lakebed material. A model cell was 100 percent or 0 percent lake. Therefore, the area across the bottom face of the river cell is 40,000 square feet (ft2; 200 ft x 200 ft), the lakebed thickness was set to 1 ft, and the vertical hydraulic conductivity of the lakebed was 0.01 foot per day (ft/day), resulting in a conductance value for each cell of 400 square feet per day. Several drainage lakes in the model

Appendix 1.  Additional Data Processing Steps to Build the MODFLOW-NWT Packages    79 47°42' 47°36' 47°30' 92°12' 92°6' 92°0' 91°54' 91°48' Hillshade base from Minnesota Department of Natural Resources, 2011 Catchment boundaries from McKay and others, 2012 6 KILOMETERS 6 MILES EXPLANATION Model domain National Hydrography Dataset Plus stream catchments Areas not routed are shown as the underlying hillshade Routed areas Streams represented using the Streamflow Routing (SFR2) package Figure 1.1.  IRUNBND for the unconsolidated zone flow package showing NHDPlus v2 (McKay and others, 2012) catchments with streams represented by streamflow routing cells in the Partridge River Basin model. Purple shading represents areas within catchments where rejected recharge and groundwater seepage was routed to the stream segment. Areas without shading show the underlying hillshade digital elevation model (DEM) and represent areas where no routing was simulated.

80    Groundwater/Surface-Water Interactions in the Partridge River Basin and Hypothetical Future Mine Pits domain occur within the stream network and were modeled as RIV cells with SFR2 cells passing through the center of the lake, as shown in figure 12. This allowed streamflow routing through lakes. To prevent spurious circular exchange between adjacent SFR2 and RIV cells, the vertical hydraulic conductiv­ ity of the streambed in SFR2 cells surrounded by RIV cells was set to 1x10-6 ft/day, causing limited groundwater exchange between these features while still allowing routed streamflow across the lakes. Groundwater discharged to RIV cells was not included in the routed streamflow, and streams that flowed through lakes may underestimate the total base flow. Representation of Mine Pits The PRB contains numerous pumped and flooded mine pits that were included in the model. The geometry of the pits was determined using a MDNR mining feature shapefile (John Coleman, Great Lakes Indian Fish & Wildlife Commission, written commun., 2018). The pit extents were modified using the pit extents visible in the 2011 lidar data (MDNR, 2011). Pumped mine pits were modeled using drain (DRN) cells (Harbaugh and others, 2000) and flooded pits were modeled as high hydraulic conductivity zones (Anderson and others, 2002; Hunt and others, 2003). Drain elevations were set to the minimum lidar elevation (MDNR, 2011) in a model cell occupied by a pumped pit. Pit drain cells were assigned to a layer based on the drain elevation relative to the model layer elevations. Because model top elevation was the average lidar elevation and the drain cell elevations were the minimum lidar eleva­ tion, most of the drain cells were below layer 1. Cells in the layers above a drain cell, representing air in the pumped pit, were assigned a high vertical and horizontal hydraulic conduc­ tivity. Model cells containing drain cells were also assigned a high hydraulic conductivity. Using a hydraulic conductivity for the pit cell representing rock properties resulted in inaccu­ rate ponding of water in the layer above the drain cell. Using a high hydraulic conductivity in the cell containing the pit drain cell resulted in a pit seepage face where water discharging along the pit walls above the drain cell was conveyed down to the drain cell. Flooded pits were modeled using a high vertical and horizontal hydraulic conductivity. The depth of flooded pits was unknown because lidar cannot penetrate water. Therefore, flooded pits were assumed to be in model layers 1 through 6. The bottom of layer 6 was 100 ft into the bedrock so the simu­ lated pit depth was 100 ft into the bedrock. The flooded pit depth was selected based on the depths of nearby pumped pits, as determined from lidar elevations. Because the flooded pit elevations are not specified explicitly as a boundary condition like the lakes or pumped pits, the flooded pit elevations from the 2011 lidar data were used as a calibration target to help constrain the bedrock properties. Partially flooded pits are pits flooded with water but also with a reported pit pumping rate. Partially flooded pits were modeled as high conductivity flooded pits with a single pumping well in the pit. The well pumping rate was assigned based on the 2011 reported pit pumping rates (Avery Guertin, MDNR, written commun., 2019). There were two partially pumped pits in the model domain: the Dunka pit on the far northeast end of the Iron Range and a part of the Peter Mitchell pit. The two partially pumped pits are shown in figure 12 as flooded pits with a red WEL cell. The water levels in flooded pits and pumped water from dry pits were assumed to include groundwater contributions to the pits and surface runoff from precipitation events. Runoff to a pit can be approximated using annual precipita­ tion and an assumed runoff coefficient, or the percentage of precipitation that becomes runoff in the surface catchment of the pit. Runoff coefficients are dependent on the driange basin cover (developed, field, wooded, etc.), soil type, and slope and are generally used to estimate runoff for a par­ ticular storm event of a given intensity using an equation like the Rational Method (Washington State Department of Transportation, 2019). A runoff coefficient multiplied by the average annual precipitation was used to estimate a steadystate runoff to each mine pit. An average annual precipitation rate of 33.3 inches/year (in/yr; calculated from the average of 1995–2015 [33.2 in/yr] and 2011–2013 [33.4 in/yr] average annual precipitation from the SLRB SWB model) was used for the 2011–2013 and 1995–2015 model runs. The runoff rate was multiplied by the assumed pit catchment area to calculate a volume of runoff. This volume was added to the pit using an injection WEL package well in the high hydraulic conductivity cells of the pit. The pit catchments were digitized using land elevations contours from 2011 lidar data (MDNR, 2011). The digiti­ zation was an approximate process because of the uneven topography of the mined surfaces. An initial runoff coefficient of 0.4 (40 percent) was used based on the Barr (2008) runoff estimates of 0.3 (30 percent) for undisturbed areas and 0.6 (60 percent) in mining features for basins near the Iron Range.

Appendix 1.  Additional Data Processing Steps to Build the MODFLOW-NWT Packages    81 References Cited Anderson, M.P., Hunt, R.J., Krohelski, J.T., and Chung, K., 2002, Using high hydraulic conductivity nodes to simulate seepage lakes: Ground Water, v. 40, no. 2, p. 117–124. Barr Engineering Company [Barr], 2008, Long-range hydrol­ ogy study: Prepared for Northshore Mining Company, accessed April 28, 2020 at ://​www.leg.st​ate.mn.us/​ docs/​2015/​other/​150681/​PFEISref_​1/​Barr%202008o.pdf. Dengler, E., and Wagner, K., 2017, Preliminary surficial geologic map of southeastern Arrowhead area, Lake and St. Louis Counties, northeastern Minnesota: Minnesota Geological Survey Open-File Report OFR2016-04, 1 p., accessed April 15, 2019, at ://​conservanc​y.umn.edu/​ handle/​11299/​183258. Dengler, E.L., Wagner, K.G., and Hamilton, J., 2017, Preliminary Quaternary geology and sand distribution models of the central Arrowhead Area, Lake and St. Louis Counties, northeastern Minnesota: Minnesota Geological Survey Open-File Report OFR2016-04, 1 p., accessed April 15, 2019, at ://​conservanc​y.umn.edu/​handle/​ 11299/​183258. Harbaugh, A.W., Banta, E.R., Hill, M.C., and McDonald, M.G., 2000, MODFLOW-2000, the U.S. Geological Survey modular ground-water model—User guide to modulariza­ tion concepts and the Ground-Water Flow Process: U.S. Geological Survey Open-File Report 00–92, 121 p. Hunt, R.J., Haitjema, H.M., Krohelski, J.T., and Feinstein, D.T., 2003, Simulating ground water-lake interactions— Approaches and insights: Ground Water, v. 41, no. 2, p. 227–237. Leaf, A.T., 2018, SFRmaker: GitHub software web page, accessed November, 27, 2018, from ://github.com/​ aleaf/​SFRmaker. Leaf, A.T., Reeves, H.W., and Fienen, M.N., 2021, SFR­ maker: U.S. Geological Survey code database web page, ://doi.org/​10.5066/​P9U2T031. McKay, L., Bondelid, T., Dewald, T., Johnston, J., Moore, R., and Rea, A., 2012, NHDPlus Version 2—User guide: NHDPlus web page, accessed February 4, 2019, at ://www.horizon-​systems.com/​NHDPlus/​ NHDPlusV2_​home.php. Minnesota Department of Natural Resources [MDNR], 2011, LiDAR elevation, Arrowhead Region, NE Minnesota, 2011: Minnesota IT Services Geospatial Information Office data­ set, accessed June 2017, at ://ww​w.mngeo.st​ate.mn.us/​ chouse/​elevation/​lidar.. Minnesota Department of Natural Resources [MDNR], 2019, DNR hydrography dataset: MDNR web page, accessed April 2019 at htt​ps://gisda​ta.mn.gov/​dataset/​water-​dnr-​ hydrography. Niswonger, R.G., and Prudic, D.E., 2005, Documentation of the Streamflow-Routing (SFR2) Package to include unsatu­ rated flow beneath streams—A modification to SFR1: U.S. Geological Survey Techniques and Methods Book 6, chap. A13, 50 p. Niswonger, R.G., Prudic, D.E., and Regan, R.S., 2006, Documentation of the Unsaturated-Zone Flow (UZF1) Package for modeling unsaturated flow between the land surface and the water table with MODFLOW–2005: U.S. Geological Survey Techniques and Methods Book 6, chap. A19, 62 p. Niswonger, R.G., Panday, S., and Ibaraki, M., 2011, MODFLOW–NWT—A Newton formulation for MODFLOW–2005: U.S. Geological Survey Techniques and Methods, Book 6, chap. A37, 44 p. Perry, M.T., 2015, rasterstats: Python package web page, accessed May 21, 2019, at ​://pythonh​osted.org/​ rasterstats/.​ Smith, E.A., 2017, Soil-water-balance model data sets for the St. Louis River drainage basin, northeast Minnesota, 1995–2010: U.S. Geological Survey data release, accessed June 2018 at ://doi.org/​10.5066/​F7Z60MJ0. Smith, E.A., and Westenbroek, S.M., 2015, Potential ground­ water recharge for the State of Minnesota using the Soil Water-Balance model, 1996–2010: U.S. Geological Survey Scientific Investigations Report 2015–5038, 85 p., accessed November 16, 2017, at ://doi.org/​10.3133/​sir20155038. Washington State Department of Transportation (WSDOT), 2019, Hydraulics manual: WSDOT web page, accessed April 28, 2020, at ​://www.wsd​ot.wa.gov/​Publications/​ Manuals/​M23-​03.htm.

82    Groundwater/Surface-Water Interactions in the Partridge River Basin and Hypothetical Future Mine Pits Appendix 2.  Estimation of Dipping Bedrock Units Bedrock in the Partridge River Basin (PRB) model domain consists of crystalline, dipping Proterozoic bedrock units. Jirsa and other (2011) mapped bedrock units at the top of the bedrock surface. Maps of the bedrock beneath the bed­ rock surface were not available; therefore, the bedrock surface from Jirsa and others (2011) was translated vertically down through the seven lower bedrock model layers by assuming a 10˚ dip to the southeast between the Biwabik Iron-Formation of Animikie Group (hereafter Biwabik Iron-Formation) and the Virginia Formation and between the Archean formations and the Biwabik Iron-Formation, and assuming a 30˚ dip to the southeast between the Virginia Formation and the Duluth Complex of Keweenawan Supergroup (hereafter Duluth Complex). The dip angles were based on reported dips of 5–15˚ for the contact between the Biwabik Iron-Formation and Virginia Formation and 20–35˚ for the contact between the Virginia Formation and the Duluth Complex (Miller and others, 2002). Bedrock layer thicknesses were assigned using the layering scheme discussed in the "Discretization” section of the report and appendix 1. Bedrock borehole data from the Minnesota County Well Index database (County Well Index, 2018) was used to adjust hydraulic conductivity zones in the bedrock layers. The County Well Index contains information on borehole depth, top and bottom of bedrock units, top of borehole loca­ tion, and borehole dip and azimuth. The borehole data were used to calculate the depth and location of the top and bottom of each bedrock unit in the model domain. The depth and location in dipping boreholes were corrected for the borehole dip. A bedrock unit was assigned to a model layer if the layer midpoint was within the range of the mapped bedrock top and bottom elevations for a bedrock unit. Starting at the mapped top of the bedrock surface and moving down one layer at a time, assumed dips were applied to the units using a python script. The resulting bedrock hydraulic conductivity zones for the model layer generated by the python script were plotted in a GIS platform with the borehole points that fell within the specific model layer. Areas along unit boundaries where the borehole data disagreed with the estimated bedrock unit were adjusted in GIS to create an updated bedrock distribution. The updated bedrock layer zones were then used to estimate the underlying bedrock zones in the python script. This process was continued downward until the bottom model layer was created. About 4,500 boreholes were used to correct the bedrock boundaries in the model layers. Not every borehole point was matched, and a few borehole points seemed mislabeled. The updating process produced a 96 percent agreement between the borehole data and model layering. The model layer and borehole data are provided for layer 8 in figure 2.1. References Cited County Well Index, 2018, The Minnesota County Well Index: Minnesota Geological Survey, University of Minnesota web page, accessed December 11, 2018, at ://m​gsweb2.mng​ s.umn.edu/​cwi_​doc/​cwidoc.htm. Jirsa, M.A., Boerboom, T.J., Chandler, V.W., Mossler, J.H., Runkel, A.C., and Setterholm, D.R., 2011, Geologic map of Minnesota—Bedrock geology: Minnesota Geological Survey State Map Series S–21, 1 sheet plus supplemental files, scale 1:500,000, accessed July 31, 2017, at ://​ conservanc​y.umn.edu/​handle/​11299/​101466. Miller, J.D., Green, J.C., Severson, M.J., Chandler, V.W., Hauck, S.A., Peterson, D.M., and Wahl, T.E., 2002, Geology and mineral potential of the Duluth Complex and related rocks of Northeastern Minnesota: Minnesota Geological Survey Report of Investigations 58, 207p., accessed October 10, 2019, at ://​conservanc​y.umn.edu/​ bitstream/​handle/​11299/​58804/​ri_​58%5b1%5d.pdf?​ sequence=​3&isAllowed=​y.

Appendix 2.  Estimation of Dipping Bedrock Units    83 47°42' 47°36' 47°30' 92°12' 92°6' 92°0' 91°54' 91°48' Model domain Partridge River Basin Duluth Complex of Keweenan Supergroup Virginia Formation Biwabik Iron-Formation of Animikie Group Archean units Duluth Complex of Keweenan Supergroup Virginia Formation Biwabik Iron-Formation of Animikie Group Archean units EXPLANATION County Well Index borehole data from dipping holes County Well Index borehole data from vertical holes Partridge River Basin Outline from USGS, 2012 Hillshade from Minnesota Department of Natural Resources, 2011 Borehole data modified from County Well Index, 2018. Archean units Biwabik Iron Formation of Animikie Group Virginia Formation Duluth Complex of Keweenan Supergroup 6 KILOMETERS 6 MILES Figure 2.1.  Layer 8 model bedrock hydraulic conductivity zones and bedrock data from the County Well Index.

Appendix 3.  Streamflow Target Processing Streamflow targets are used in the model calibration pro­ cess to compare model-calculated streamflow with measured streamflow data. Modeled streamflow represents base-flow conditions and does not include storm runoff or snowmelt. Available streamflow data include eight 2018 USGS synoptic flow survey measurements available on the National Water Information System database (NWIS, 2019; station number listed in the Haserodt and others, [2021] data release) and one continuous streamgage record from the Minnesota Department of Natural Resources (MDNR) Cooperative Stream Gaging database (://​www.dnr.st​ate.mn.us/​waters/​csg/​index.; MDNR, 2018) streamgage 3149002. Base-flow separation techniques were used to estimate continuous base flow from the continuous streamflow record streamgage 3149002. Synoptic streamflow measurements only represent stream conditions on a single day and were adjusted to reflect average base-flow conditions during the model period. This appendix describes how both the continuous record and the synoptic streamflow measurements were processed to represent baseflow conditions. Data from MDNR streamgage 3149002 is the only con­ tinuous streamflow record within the Partridge River Basin; the streamgage data start in June 2009. For MDNR streamgage 3149002, average base flows were computed from continuous records by using the Institute of Hydrology Base Flow Index (BFI) method (Institute of Hydrology, 1980; Wahl and Wahl, 1988). These base flows were averaged from 2009 to 2015 for the 1995–2015 average conditions model and from 2011 to 2013 for the 2011–2013 mining conditions model. Synoptic flow measurements are measurements made on a single day throughout the basin, and these measure­ ments were converted to an equivalent base flow using the following methods. The time of the measurement likely does not represent an average base flow for the model time peri­ ods, which include several years. Therefore, single synoptic measurements were adjusted to reflect average conditions using the base-flow record for the nearest streamgage with continuous data for the full model period relative to the day the synoptic measurement was collected. No streamgage data were available in or near the Partridge River Basin for the entire 1995–2015 period. Due to potential effects from the impounded Colby Lake and the period of record, the MDNR streamgage at the Colby Lake outlet was not used to adjust the synoptic flow data. The closest streamgage not affected by a dam is the USGS streamgage (NWIS, 2019) for the St. Louis River near Skibo (04015438) with a period of record starting in 2011. The synoptic streamflow data were adjusted using the ratio of the daily estimated base flow at the Skibo streamgage on the day of the synoptic survey compared to the average base flow during 2011–2015 (for the 1995–2015 average conditions model) or 2011–2013 (for the 2011–2013 mining conditions model). This ratio was then used to adjust the syn­ optic flow values to a value that represented average base-flow conditions during the time of the study. References Cited Haserodt, M.J., Hunt, R.J., Fienen, M.N., and Feinstein, D.J., 2021, MODFLOW-NWT and MODPATH models, capture zones and uncertainty data analysis for the Partridge River Basin, Minnesota, U.S. Geological Survey data release, ://doi.org/​10.5066/​P9VODOU8. Institute of Hydrology, 1980, Low flow studies report no. 3— Research report: Wallingford, Oxon, United Kingdom, Institute of Hydrology Report no. 3, p. 12–19. Minnesota Department of Natural Resources [MDNR], 2018, Partridge River nr Hoyt Lakes, 2.3mi us of CSAH110 (03149002) dataset: MDNR web page, accessed December 17, 2018, at ://​www.dnr.st​ate.mn.us/​waters/​ csg/​site_​report.?​mode=​getsitereport&site=​03149002. U.S. Geological Survey, 2019, USGS water data for the Nation: U.S. Geological Survey National Water Information System database, accessed March, 2019, at ://doi.org/​10.5066/​F7P55KJN. Wahl, K.L., and Wahl, T.L., 1988, Effects of regional ground­ water level declines on streamflow in the Oklahoma Panhandle in Proceedings of the Symposium on WaterUse Data for Water Resources Management: Bethesda, Maryland, American Water Resources Association, p. 239–249. 84    Groundwater/Surface-Water Interactions in the Partridge River Basin and Hypothetical Future Mine Pits

Appendix 4.  MODPATH and Monte Carlo Setup for Capture Zone Analysis    85 Appendix 4.  MODPATH and Monte Carlo Setup for Capture Zone Analysis A capture zone represents the area contributing recharge to any hydrologic sink such as a well, spring, or a mine pit (Franke and others, 1999) and can be calculated by forward tracking particles from the water table (for example, Hunt and Steuer, 2000; Juckem and others, 2014). The Partridge River Basin (PRB) project used Monte Carlo techniques to gener­ ate estimated capture zones that include the probability of capture of recharge derived from each model cell for a given model feature. Monte Carlo techniques use randomly gener­ ated model parameter sets that are sampled from a plausible parameter range to create numerous model runs. This appendix describes the MODPATH settings used for the particle track­ ing and the Monte Carlo setup used to estimate probabilistic capture zones in the PRB. MODPATH, version 7 (Pollock, 2016), was used to do forward particle tracking for the PRB model scenarios with a particle placed in every model cell. The MODPATH option to “stop at weak sink” was used for “internal” sinks (for exam­ ple, wells, drains, and general head cells). The “stop at weak sink” option is used to determine if a particle will stop or con­ tinue through a sink cell that has both inflows and outflows (a “weak sink cell”). In contrast, a “strong sink cell” has inflows and no outflows. For “external” cells intersecting the water table (streams, land surface), the IFACE option in MODPATH was selected which eliminates the weak-sink problem and allows for an unambiguous determination of capture (Pollock, 2016). When a particle terminates at an external face (typically the top face in correspondence with a stream, the water table, or the land surface), it is considered captured. For the Monte Carlo realizations, random model param­ eter sets were generated using the ParameterEnsemble. from_gaussian_draw function in the PyEMU python pack­ age (://github.com/​jtwhite79/​pyemu; White and others, 2016). This function utilizes the parameter estimation files used for the based model calibration with PEST (Doherty, 2018a, 2018b, 2020). A posterior covariance calculated from the Schur Complement was used to adhere to associations between correlated parameters in the generated model param­ eter sets. The parameters were drawn from reasonable bounds for each model parameter based on literature values. The PEST control file with parameter bounds used for the Gaussian parameter draw is in the ancillary directory of model archive associated with this report (Haserodt and others, 2021). For each model realization, PEST was used to run model scenarios of the generated parameter sets with a single forward run of the two base models (1995–2015 average conditions and 2011–2013 mining conditions models) used in the model calibration and the models representing the five mining model scenarios. The reasonableness of each realization parameter set was evaluated using the misfit (or “phi”) between the modeled measured calibration targets. Realizations with suf­ ficiently large misfits were removed from the probabilistic capture zone analysis (also known as realization “condition­ ing”). Out of 148 realizations, 127 had a sufficiently good fit with a phi value below 1,400; as reference, the phi of the selected base model calibration for this project was 995 and represents a model that fit the calibration data reasonably well without extreme parameter values. References Cited Doherty, J., 2018a, PEST, Model-independent parameter estimation—user manual part II: PEST utility support soft­ ware (7th ed.): Brisbane, Australia, Watermark Numerical Computing, 267 p., accessed April 28, 2020, at ://w​ ww.pesthom​epage.org/​Downloads.php [Also available at ​://pesthom​epage.org/​documentation.]. Doherty, J., 2018b, PEST, Model-independent parameter estimation—user manual part I: PEST, SENSAN and Global Optimisers (7th ed.): Brisbane, Australia, Watermark Numerical Computing, 393 p., accessed April 28, 2020, at ://w​ww.pesthom​epage.org/​Downloads.php [Also avail­ able at ​://pesthom​epage.org/​documentation.]. Doherty, J., 2020, PEST for highly parallelized computing environments: Brisbane, Australia, Watermark Numerical Computing, 88 p., accessed April 28, 2020, at ://w​ ww.pesthom​epage.org/​Downloads.php. Franke, O.L., Reilly, T.E., Pollock, D.W., and LaBaugh, J.W., 1999, Estimating areas contributing recharge to wells, lessons from previous studies: U.S. Geological Survey Circular 1174, 20 p. Haserodt, M.J., Hunt, R.J., Fienen, M.N., and Feinstein, D.J., MODFLOW-NWT and MODPATH models, capture zones and uncertainty data analysis for the Partridge River Basin, Minnesota, US Geological Survey data release, ://doi.org/​10.5066/​P9VODOU8. Hunt, R.J., and Steuer, J.J., 2000, Simulation of the recharge area for Frederick Springs, Dane County, Wisconsin: USGS Water-Resources Investigations Report 00–4172, 33 p. Juckem, P.F., Fienen, M.F., and Hunt, R.J., 2014, Simulation of groundwater flow and interaction of groundwater and surface water on the Lac du Flambeau Reservation, Wisconsin: U.S. Geological Survey Scientific Investigations Report 2014–5020, 34 p. Pollock, D.W., 2016, User guide for MODPATH Version 7—A particle-tracking model for MODFLOW: U.S. Geological Survey Open-File Report 2016–1086, 35 p. White, J.T., Fienen, M.N., and Doherty, J.E., 2016, A python framework for environmental model uncertainty analysis: Environmental Modelling & Software, v. 85, p. 217–228.

Appendix 5.  Data Worth Setup Data worth analysis estimates how potential new obser­ vations can reduce uncertainty in a forecast; here the linear uncertainty analysis framework was extended to perform data worth analysis following the methods of Doherty (2015) and Fienen and others (2010). In this approach, data worth is the normalized difference between the forecast uncertainty with and without the new hypothetical data and represents how new data may improve the calibration of the base model such that the forecast values in the model scenario may become more certain. This appendix describes how data worth was set up. The data worth was performed using the get_added_ obs_importance function in the python package PyEMU (://github.com/​jtwhite79/​pyemu; White and others, 2016). Potential hydraulic head observations were given a weight of 0.2 feet-1 or an assumed standard deviation in the hydraulic head observation of 5 feet (​1 foot

0.2  ​); this is the same weight assigned to the highest quality hydraulic head observations in the calibration dataset. Because most of the hydraulic conduc­ tivities in the model were assigned using a zone approach— with the exception of the hydraulic conductivity of the bed­ rock in layers 5 and 6 where pilot points were used—a grid of pilot points (fig. 5.1) representing multipliers on the hydraulic conductivity arrays was added using PyEMU. Without these pilot points, the data worth for a forecast may produce spurious results overly affected by the bound­ ary between two zones; this artifact of the model zoning is discussed in detail in Fienen and others (2010). Adding pilot point multipliers across the hydraulic conductivity array provides additional parameter flexibility by considering effects of new data on different areas within a zone. The new data considered in the data worth analysis for the Partridge River Basin model were a grid of hypothetical potential new head observations in the center of every fifth model cell starting three cells in from the model upper right corner and in each of the odd layers (1, 3, 5, 7, 9, and 11). The data worth analysis was performed for the following model forecasts: (1) the average depth to water for 12 wet­ lands, (2) 6 streamflow locations, and (3) 5 pit inflow forecasts (fig. 33). The forecasts were evaluated in each of the five mining scenarios (table 4). A data worth calculation was made for each potential hydraulic head observation in all the odd model layers, for every forecast under each mining scenario. The result is a grid of data worth with a value of 0–1 assigned to the potential monitoring points. A value of 0 represents no decrease in forecast uncertainty by adding that hydraulic head observation and 1 represents a 100 percent reduction in forecast uncertainty by adding that hydraulic head observation. The data worth at each point was then interpolated to make a raster of data worth across the model domain for a given layer, forecast, and mining scenario. References Cited Doherty, J., 2015, Calibration and uncertainty analysis for complex environmental models: Brisbane, Australia, Watermark Numerical Computing, 227 p. Fienen, M.N., Doherty, J.E., Hunt, R.J., and Reeves, H.W., 2010, Using prediction uncertainty analysis to design hydrologic monitoring networks—Example applications from the Great Lakes Water Availability Pilot Project: U.S. Geological Survey Scientific Investigations Report 2010–5159, 44 p., ht​tps://pubs​.usgs.gov/​sir/​2010/​5159/​. White, J.T., Fienen, M.N., and Doherty, J.E., 2016, A python framework for environmental model uncertainty analysis: Environmental Modelling & Software, v. 85, p. 217–228. 86    Groundwater/Surface-Water Interactions in the Partridge River Basin and Hypothetical Future Mine Pits

Appendix 5.  Data Worth Setup Hillshade base from Minnesota Department of Natural Resources, 2011 47°42' 47°36' 47°30' 92°12' 92°6' 92°0' 91°54' 91°48' 8 KILOMETERS 8 MILES Lakes represented with the RIV package Streams represented using the Streamflow Routing (SFR2) package Model domain EXPLANATION Hydraulic conductivity pilot points Hypothetical head observations Figure 5.1.  Pilot point and potential head observation locations for the data worth analysis.

For additional information, contact: Director, Upper Midwest Water Science Center U.S. Geological Survey 8505 Research Way Middleton, WI 53562 ://www.usgs.gov/centers/umid-water

Haserodt and others—Groundwater/Surface-Water Interactions in the Partridge River Basin and Hypothetical Future Mine Pits —SIR 2021–5038 ISSN 2328-0328 (online) ://doi.org/ 10.3133/ sir20215038