Water Resources of Michigan

Ground-Water Flow and Contributing Areas to Public -Supply Wells in Kingsford and Iron Mountain, Michigan

US Geological Survey Water-Resources Investigation 00-4226
Lansing, Michigan 2001

By: Carol L. Luukkonen and D. B. Westjohn

Table of Contents including Figures, Maps, Graphs, Tables, Appendix, Conversion Factors and Vertical Datum, and Additional Information.
http://mi.water.usgs.gov/pubs/WRIR/WRIR00-4226/WRIR00-4226TOC.php

ABSTRACT

The cities of Kingsford and Iron Mountain are in the southwestern part of Dickinson County in the Upper Peninsula of Michigan. Residents and businesses in these cites rely primarily on ground water from aquifers in glacial deposits. Glacial deposits generally consist of an upper terrace sand-and-gravel unit and a lower outwash sand-and-gravel unit, separated by lacustrine silt and clay and eolian silt layers. These units are not regionally continuous, and are absent in some areas. Glacial deposits overlie Precambrian bedrock units that are generally impermeable. Precambrian bedrock consists of metasedimentary (Michigamme Slate, Vulcan Iron Formation, and Randville Dolomite) and metavolcanic (Badwater Greenstone and Quinnesec Formation) rocks. Where glacial deposits are too thin to compose an aquifer usable for public or residential water supply, Precambrian bedrock is relied upon for water supply. Typically a few hundred feet of bedrock must be open to a wellbore to provide adequate water for domestic users. Ground-water flow in the glacial deposits is primarily toward the Menominee River and follows the direction of the regional topographic slope and the bedrock surface.
To protect the quality of ground water, Kingsford and Iron Mountain are developing Wellhead Protection Plans to delineate areas that contribute water to public-supply wells. Because of the complexity of hydrogeology in this area and historical land-use practices, a steady-state ground-water-flow model was prepared to represent the ground-water-flow system and to delineate contributing areas to public-supply wells. Results of steady-state simulations indicate close agreement between simulated and observed water levels and between water flowing into and out of the model area. The 10-year contributing areas for Kingsford's public-supply wells encompass about 0.11 square miles and consist of elongated areas to the east of the well fields. The 10-year contributing areas for Iron Mountain's public-supply wells encompass about 0.09 square miles and consist of elongate areas to the east of the well field.

INTRODUCTION

In 1997, the U.S. Geological Survey, in cooperation with the cities of Kingsford and Iron Mountain, began a 2-year study to describe ground-water flow and determine the areas that contribute water to the cities' public-supply wells. The cities of Kingsford and Iron Mountain are in the southwestern part of Dickinson County, in the Upper Peninsula of Michigan (fig. 1). Both of these cities currently have four municipal wells in glacial aquifers that provide ground water to residents and businesses for water supply (fig. 2). In an effort to protect the quality of ground water that is withdrawn by these wells, the cities of Kingsford and Iron Mountain are developing Wellhead Protection Plans (WHPP). As part of WHPP, the cities will delineate areas that contribute water to public-supply wells. Kingsford's wells are located near Cowboy Lake and near Ford Airport (fig. 2). Iron Mountain's wells are located about 2 miles north of Kingsford's well fields. Because of the complexity of hydrogeology in these areas and historical land-use practices, a computer model is needed to describe the ground-water-flow system and to delineate areas that contribute water to public-supply wells.
This report describes the hydrology of the Kingsford-Iron Mountain area and the simulation of ground-water flow in the Kingsford-Iron Mountain area. During this study, water levels were measured in 28 domestic, irrigation, and observation wells to determine ground-water-flow directions. Thirteen aquifer tests performed previously to this study were reanalyzed and one aquifer test was conducted during this study to determine hydraulic conductivity information for the ground-water-flow model. Particle-tracking analysis was used in conjunction with flow simulations to delineate the land surface and subsurface areas that contribute water to the public-supply wells. Model simulations are based on the U.S. Geological Survey Modular Three-Dimensional Finite-Difference Ground-Water Flow Model, MODFLOW-96 (McDonald and Harbaugh, 1988, 1996), which is a generalized code for ground-water-flow modeling. Areas contributing water to supply wells were determined using a particle-tracking post-processor package for MODFLOW (MODPATH, Pollock, 1989).

Figure 1. Location of Kingsford and Iron Mountain in the Upper Peninsula of Michigan.

Figure 2. Location of public -supply wells in the cities of Kingsford and Iron Mountain, Michigan.

Previous Studies

Numerous studies have contributed to an understanding of the geology and hydrology of this area. Kingsford and Iron Mountain have evaluated potential areas for well-field development. These evaluations include test borings, aquifer tests, and geophysical surveys (Keck Consulting Services, Inc., 1974; Coleman Engineering Co., 1975, 1980, 1982, and 1987; Lynn M. Miller, P.E., 1976; Sundberg, Carlson, and Associates, Inc., 1989; and McNamee, Porter and Seeley, 1981 and 1990).
Two studies in the area (Engineering for Earth, Water, and Air Resources [EWA], 1986; 1987) focused on hydrogeologic conditions at a site in Kingsford that was used for manufacturing facilities. In this area, glacial kettles (depressions in the land surface formed by melting of ice blocks and subsidence of overlying deposits) were used as retention ponds for waste products. EWA's investigations included installation of monitoring wells and soil borings, and development of a ground-water-flow model to determine the potential effect of these waste-disposal pits on water quality, and the potential effect on Kingsford's water supply. Subsequent investigations have been conducted and indicate that although environmental contamination is present near historic industrial sites, it is unlikely that these areas pose any threat to the quality of water at either the Kingsford or Iron Mountain well fields.

Acknowledgments

The authors gratefully acknowledge the assistance of Darryl Wickman, Tony Edlebeck,

George Groeneveld, and Al Muntz, city of Kingsford; and Jim Urbany and John Ellis, city of Iron Mountain, for providing information and access to the public-supply wells and providing other sources of information about the cities; Chuck Thomas, Drinking Water and Radiological Protection Division, Michigan Department of Environmental Quality, for providing copies of reports and information about geologic and hydrologic conditions in the areas; and Shirley Businski, Land and Water Management Division, Michigan Department of Environmental Quality, for providing information on public-supply wells. Sincere thanks also are given to homeowners and golf course owners, who permitted access to their wells for water-level measurements.

DESCRIPTION OF STUDY AREA

The cities of Kingsford and Iron Mountain are bounded by the Menominee River to the south and west, and by Trader Hill to the north (fig. 2). The Kingsford-Iron Mountain study area is within the Menominee River drainage basin which includes Dickinson County, Mich. and Florence County, Wisc. Historically, iron ore and timber industries were important to the area, while currently lumber and paper industries, tourism, and recreational activities support the area's economy. The geologic and hydrologic characteristics that affect ground-water flow in the study area are described below.

Geologic Setting

Most of the Kingsford-Iron Mountain area is underlain by Precambrian metasedimentary (Michigamme Slate) and metavolcanic (Badwater Greenstone and Quinnesec Formation) rocks. The Michigamme Slate consists of fine- to medium-grained metagraywacke and slate in massive or graded beds, whereas the Badwater Greenstone consists of massive, chloritized mafic volcanic rocks (Sims, 1990). The Quinnesec Formation consists of mafic and felsic volcanic rocks (Oakes and others, 1973). In the vicinity of Pine Mountain, Millie Hill, and Trader Hill (fig. 3), bedrock consists of metasedimentary rocks of the Vulcan Iron Formation, with alternating iron-rich and quartz-rich layers, and the Randville Dolomite, consisting of massive dolomite with detrital quartz grains and thick-bedded crystalline dolomite (Sims, 1990).
Bedrock is overlain by Pleistocene glacial deposits, except where bedrock is exposed along the Menominee River to the south and to the north near Spread Eagle, Wisconsin (Clayton, 1986). There also are isolated areas where bedrock is exposed (for example, Ford Airport), but these are of small areal extent. Glacial deposits range from 0 to about 380 ft in thickness. The Marenisco Moraine (Martin, 1955) forms a distinct ridge through the Kingsford-Iron Mountain area (fig. 3). Approximately parallel to this moraine and south of the Menominee River lies another moraine informally referred to as the Menominee Moraine.
Within the study area, glacial deposits generally consist of an upper terrace sand-and-gravel unit and a lower outwash sand-and-gravel unit separated by intermediate lacustrine silt and clay and eolian-silt layers (figs. 4, 5, and 6). These units are not regionally continuous and are absent in some isolated areas within Kingsford as indicated by boring logs. Towards the west and north, information obtained from well and geophysical logs indicates that upper terrace and lower outwash sand-and-gravel units are separated primarily by silt and fine sand.
The upper sand-and-gravel unit ranges from 0 to 160 ft in thickness and consists of medium to fine-grained sand with varying amounts of silt and gravel (Engineering for Earth, Water, and Air Resources, Inc., 1987). The confining unit ranges from 0 to 150 ft in thickness and consists of well sorted and well compacted fine silt and clay. In some areas, the silt and clay are absent, possibly due to partial burial and protrusion of stranded ice blocks that prevented silt deposition (Engineering for Earth, Water, and Air Resources, Inc., 1987). The lower sand-and-gravel unit ranges from 0 to 237 ft in thickness and consists of moderately well sorted fine to coarse sand with layers and lenses of silt and gravel (Engineering for Earth, Water, and Air Resources, Inc., 1987). The lower sand-and-gravel unit is underlain by a lodgement till or the Vulcan Iron Formation over much of the area. The lodgement till and Vulcan Iron Formation units range from 0 to 250 ft in thickness. West of the Menominee River, glacial deposits consist primarily of till or fluvial sediments that change conspicuously within relatively short distances (Clayton, 1986). The sediments consist of sand and gravel, which were deposited by glacial melt-water streams. The till consists of reddish brown or dark brown gravelly sandy loam with abundant surface boulders (Clayton, 1986, plate 1).

Figure 3. Location of well fields and physical featues in the Kingsford-Iron Mountain area, Michigan

Figure 4. Location of lines of sections, test wells (TW), and supply wells (PW) in the Kingsford-Iron Mountain area, Michigan

Figure 5. North-south section view of Kingsford's well field, Michigan. (Line of section show in figure 4)

Figure 6. North-south section view of Iron Mountain's well field, Michigan (Line of section shown in figure 4)

Hydrologic Setting

Precipitation averages about 29 in/yr over the Kingsford-Iron Mountain area. Mean monthly precipitation ranges from a minimum in February of 0.9 in to a maximum in August of 3.7 in. The study area receives about 62.8 in/yr of snowfall (Midwestern Climate Center, 1999). Precipitation averages about 30 in/yr in the Pine-Popple River Basin, west of the Menominee River, and is highest from late spring to early autumn (Oakes and others, 1973).
Two glacial aquifers (upper and lower aquifers) underlie the study area, based on previous studies and available well and boring logs. The upper aquifer in the glacial deposits consists of the upper terrace sand-and-gravel unit, whereas the lower aquifer in the glacial deposits consists of the lower outwash sand-and-gravel unit. The intermediate silt layer is a fairly effective barrier to hydraulic communication between the two aquifers. This lack of hydraulic communication is strongly indicated by substantial head differences between the two aquifers; water levels in shallow wells typically are 20-60 ft higher than water levels in deep wells (D.B. Westjohn, oral commun., 2000). Where the silt layer is absent in the Kingsford area, the upper and lower sand-and-gravel units compose one aquifer (Engineering for Earth, Water, and Air Resources, Inc., 1987). Hydraulic conductivities for the upper aquifer in the Kingsford area range from 9 to 0.9 ft/d based on tests performed by Engineering for Earth, Water, and Air Resources, Inc. (1987). Hydraulic conductivities average about 9 ft/d for the lower aquifer and range from 0.001 to 0.5 ft/d for the intermediate confining unit (Engineering for Earth, Water, and Air Resources, Inc., 1987).
An understanding of the hydraulic characteristics of glacial aquifers in the public-supply well field areas has been obtained from aquifer tests performed during investigations of potential areas for well-field development (fig. 2, table 1). Aquifer tests of wells in the lower aquifer indicate ground-water leakage from overlying units. Iron Mountain's Well #4 is screened within the lower aquifer; however, in the immediate area surrounding well #4 the upper aquifer is missing.
Ground-water flow within the Kingsford-Iron Mountain area is primarily towards the surface-water bodies to the west and southwest and follows the topographic slope and the bedrock surface. Generally, the topographic slope follows the bedrock surface except within the central part of Kingsford where there is a bedrock low. Ground-water flow is probably confined primarily within the permeable sand-and-gravel units because Precambrian bedrock units appear to be largely unfractured and limit downward migration of ground water. However, bedrock consists of metasedimentary rocks in the vicinity of Pine Mountain and Millie Hill, and ground-water flow is locally controlled by these units where they are fractured.

CONCEPTUAL MODEL

The conceptual model describes ground-water flow within the Kingsford-Iron Mountain and surrounding areas and includes a definition of the position and thickness of aquifer and confining units, boundaries, and the flow system. Geologic units within the Kingsford-Iron Mountain area can be divided into layers that control ground-water flow.
Four units were used to conceptualize the ground-water-flow system although not all units are present throughout the study area. These units consist of an upper permeable sand-and-gravel unit, an underlying confining unit consisting of silt and clay, a lower permeable sand-and-gravel unit, and a lowermost unit consisting of sandstone, dolomite, or lodgement till (where they are present). The bedrock of the Vulcan Iron Formation and the Randville Dolomite and the lodgement till of the glacial unit are poorly permeable and were combined in the lowermost unit because of their similar hydraulic properties.

Although these units are poorly permeable, some domestic wells withdraw water from the bedrock; however, a few hundred feet of borehole must be open to a wellbore to provide an adequate supply. Geologic information obtained from logs of public-supply and domestic wells and geophysical logs were used to delineate the extent of the ground-water-flow system units.
Ground-water-flow systems have physical boundaries, which are formed by the presence of an impermeable body of rock or a large body of water, or hydrologic boundaries, which include ground-water divides and streamlines. Regional physical boundaries to the flow system are formed by the boundaries of the Menominee River drainage basin that extends west into Wisconsin and east into Michigan. Boundaries for the local-flow system in the Kingsford-Iron Mountain area consist primarily of bedrock and surface-water bodies in Michigan and Wisconsin that affect flow to the Menominee River. Locally, ground-water flow is affected by physical boundaries formed by bedrock highs south of the Menominee River and north of Iron Mountain's well fields, and by surface water bodies including Lake Antoine on the east, and the Pine and Little Popple Rivers on the west. Ground-water-flow divides are assumed to coincide with the surface-water divides which follow topographic highs and the impermeable bedrock in this area. Local hydrologic boundaries consist of streamlines, which are perpendicular to the ground-water-flow direction along the east and northeast. The lower boundary is formed by the upper surface of the Michigamme Slate, Badwater Greenstone, and Quinnesec Formations and is considered to be poorly permeable.
A preliminary understanding of ground-water-flow directions can be made on the basis of surface-water elevations taken from topographic maps of the area, ground-water levels measured during this study, and previous reports covering part of the study area. For this study, water levels were measured in 28 domestic, irrigation, and observation wells. Water levels recorded on well logs and average drawdown measurements obtained for each of the public-supply wells also were analyzed to understand ground-water-flow directions. Regional ground-water flow is towards the Menominee River and follows the direction of the regional topographic slope and the bedrock surface (fig. 7).

Figure 7. Ground water levels in the upper aquifer in the Kingsford - Iron Mountain area, Michigan

SIMULATION OF GROUND-WATER FLOW

Simulation of ground-water flow is made possible by first developing a conceptual model of the flow system and then developing a numerical model that is consistent with the conceptual model. The model area consists of the Kingsford-Iron Mountain area and surrounding parts of Dickinson County, Michigan and Florence County, Wisconsin. This larger area (larger than the cities of Kingsford and Iron Mountain) was modeled so as to eliminate boundary effects on the ground-water-flow solution in the interior portion of the model area by allowing natural physical and hydrologic boundaries to be chosen, as well as to investigate whether ground water moved beneath the Menominee River as possibly indicated by the high change in ground-water levels from Pine Creek to the Menominee River.

Numerical Model

The U.S. Geological Survey Modular Three-Dimensional Finite-Difference Ground-Water Flow Model (MODFLOW-96, McDonald and Harbaugh, 1996) was the computer code used to simulate ground-water flow in the Kingsford-Iron Mountain area.

This code allows the simulation of steady-state ground-water flow in two dimensions with leakage between model layers; no ground-water storage or temporal discretization terms are required. All water entering the model area through the boundaries or as recharge is assumed to leave the model area through the boundaries, rivers, or wells. The Kingsford and Iron Mountain wells have been in operation for 8 to 25 years; therefore the ground-water flow system is considered to be in steady state and water is not lost or gained from storage.
Various assumptions were made in the numerical simulation of ground-water flow. Recharge to ground water was assumed to differ over the model area depending on the presence or absence of glacial till. Recharge is assumed to be zero in areas where till or bedrock is at or very near the surface. In other areas recharge initially was assumed to be 10 in/yr. The assumption of no recharge to areas of till and bedrock highs was investigated during model calibration. Evapotranspiration of ground water also was assumed to be zero because recharge was applied to the upper surface of the model (water table), where the effects of evapotranspiration are expected to be small, and because of the regional scale of the model. Ground water was assumed to flow primarily through the glacial sand-and-gravel units. Glacial deposits are known to vary considerably with varying lithologies and thicknesses over short distances. This variability makes exact representation of the detailed hydrogeology impossible in the numerical model. Thus, hydraulic properties of the glacial deposits are generalized to represent the regional ground-water-flow system. The underlying lodgement till and bedrock units of the Vulcan Iron Formation and Randville Dolomite were simulated as a low permeability layer. The lowermost Precambrian bedrock units were simulated as impermeable boundaries. North of the Chapin Mine, ground water is withdrawn from thHamilton Mine Shaft that extends deep into the underlying bedrock units and is discharged into Lake Antoine (fig. 3). A probable hydraulic connection is present between the bedrock and Lake Antoine based on that water levels in the mine shaft only can be lowered approximately 2 feet even with extensive pumping and the high gradient from the lake to the shaft. It is possible that the withdrawal of ground water from this mine shaft affects local ground-water levels in the immediate vicinity; however, this area is over 1.5 mi from the public-supply well fields and is not thought to affect hydrologic conditions there. Therefore, these conditions were not simulated in the model.


Model Grid and Layers

The modeled area is rectangular and consists of about 64 mi2 (fig. 8). The area is 6.7 mi long (north-south) and 9.5 mi wide (east-west). This area is horizontally discretized into a variably spaced grid of cells in 213 columns and 204 rows. In the central portion of the model, each cell is 100 by 100 ft. Cell spacing increases initially by a factor of 1.25 for 5 rows/columns and then by a factor of 1.5 out to a maximum grid spacing of 1,029 ft. This maximum spacing was selected so as not to exceed a maximum row to column spacing ratio of 10:1. Each grid cell represents the average aquifer properties in the volume of aquifer represented by the cell; any variations in properties that are within a grid cell only can be represented in the model averages.
The model area is vertically discretized into three model layers. Layer 1 represents the upper aquifer in the glacial deposits. This aquifer is under unconfined conditions with water levels representing the water table. Layer 2 represents the lower aquifer in the glacial deposits. The lower aquifer represented in layer 2 is simulated as convertible and can be either confined or unconfined depending on the position of the water table. The intermediate confining unit that separates the upper and lower aquifers was included only in the leakance term between layers 1 and 2. The bottom model layer (layer 3) represents either metasedimentary rocks (in the vicinity of Pine Mountain, Millie Hill, and Trader Hill) or lodgement till (over the remaining area) and also is represented as convertible.

Land-surface information for most of the model area is available from a digital 1:100,000 topographic quadrangle map for Iron Mountain. This information was interpolated to a grid with a 50 ft spacing. Additional point-altitude data were taken from the topographic quadrangle map for Iron Mountain SW. This information was contoured over the model area to create land-surface elevations, which constitute the top of glacial deposits (fig. 9). The land surface is characterized by steep gradients in the vicinities of Pine Mountain, Millie Hill, and Trader Hill.
Layer surfaces were determined based on logs of domestic and public-supply wells obtained from Geological Survey Division, Michigan Department of Environmental Quality; Statewide Ground-Water Database, MIRIS-Geologic Resource Mapping Unit, Michigan Department of Environmental Quality; and Wisconsin Geologic and Natural History Survey (figs. 10, 11,12). The top of the silt unit formed the bottom of layer 1 and the silt unit thickness was assigned to layer 2. In areas where the geologic materials composing a layer are absent, the layer is assigned a minimal thickness and the hydraulic properties of the overlying layer. Information for the Wisconsin part of the model area did not differentiate the glacial deposits into separate vertical units; therefore, on the western side of the Menominee River, glacial deposits are modeled as the upper layer with layer 2 assigned a 1 ft thickness and the hydraulic conductivity of the overlying layer. Additional delineation of layer surfaces in areas of sparse data was done by extrapolating from known nearby points, primarily within Wisconsin where few well logs are available. In the vicinity of Pine Mountain, Millie Hill, and Trader Hill, till overlies the bedrock units and the sand-and-gravel units are missing. Outcrops provide additional control for the delineation of the bedrock surface, especially in Wisconsin where outcrops indicate that glacial deposits are thin.

Figure 8. Model grid for simulating the ground-water-flow system underlying the Kingsford-Iron Mountain area, Michigan

Figure 9. Estimated land surface elevations in the Kingsford-Iron Mountain area, Michigan.

Figure 10. Estimated configuration of the bottom of the upper aquifer model (model layer 1) in the Kingsford-Iron Mountain area, Michigan.

Figure 11. Estimated configuration of the bottom of the lower aquifer (model layer 2), Kingsford-Iron Mountain area, Michigan

Boundary Conditions

Model boundaries extend approximately 1 mi north of the Iron Mountain well field, approximately 1 mi south of the Kingsford well field, and approximately 2 mi east and 3 mi west of the well fields (fig. 8). These boundaries are based on surface-water elevations and topographic information, and are represented by specified head, specified flux, and head-dependent flux conditions. Specified head boundaries (also referred to as constant head boundaries) are modeled by specifying head values that do not change during numeric simulation. Specified flux boundaries can be modeled by specifying the flux equal to zero to represent ground-water-flow divide boundaries and streamlines. Head-dependent flux boundaries are modeled by specifying a conductance, which limits the amount of water that can enter or leave any cell.
Specified head boundaries are used to represent hydraulic conditions along parts of the east and west model boundaries and are based on surface-water elevations (fig. 13). Specified head boundaries also were used to represent the southernmost part of the Menominee River. Specified flux boundaries with flux equal to zero are used to represent natural physical boundaries formed by outcropping or shallow bedrock along parts of the north and south model boundaries. Specified flux boundaries with flux equal to zero also are used to represent ground-water-flow divides and streamlines along parts of the north and west boundaries. Finally, head-dependent flux boundaries are used to represent Lake Antoine along part of the east boundary and three small lakes along the west boundary.
In areas where a natural physical boundary does not exist, specified head boundaries located far from the well fields were used, based on available water-level data. After model development, sensitivity of model results to boundary conditions were analyzed. Fluxes through cells forming specified head boundaries were compared for simulations in which public-supply wells were being pumped and when wells were idle. As an additional check on the ground-water-flow solution, specified head boundaries were changed to specified flow boundaries and the solutions compared. The effect of stresses induced by pumping (on the assumption of no-flow across streamlines) also was investigated.
All boundaries of layer 2 are the same as those of layer 1, except that no-flow boundaries are located in the area of Lake Antoine, because the lake does not penetrate the lower sand-and-gravel unit. Layer 3 is bounded by no-flow boundaries. In layers 1 and 2, the Chapin Mine area of the model is represented by inactive (no-flow) cells and by head-dependent flux cells in layer 3. In this area, glacial deposits are relatively thin and unsaturated, and ground water is present only in the bedrock unit below the glacial deposits.

Figure 12. Estimated configuration of the bottom of the lodgement till and metasedimentary rocks (model layer 3), Kingsford-Iron Mountain area, Michigan

Figure 13. Boundary of model area in the Kingsford-Iron Mountain area, Michigan.

Hydraulic Properties

Hydraulic properties used in model simulations include layer hydraulic conductivities and leakances, recharge rates, and streambed conductance. Aquifer hydraulic properties and leakances affect ground-water flow through and between model layers. Recharge rates indicate the amount of water movement through the upper surface of the model. Streambed conductance affects vertical flow of ground water from an aquifer to a stream or from a stream to an aquifer.
Based on interpretation of well and geophysical logs, eleven zones were used to describe conditions within the three layers (table 2, figs. 14, 15, 16). The upper and lower aquifers were separated into different zones based on calculated horizontal hydraulic conductivities from aquifer-test analyses. Transition zones were added to reduce the difference in hydraulic conductivities between zones in some areas. Initial estimates of hydraulic conductivity for the zones within each layer were based on calculated hydraulic conductivities from aquifer-test analyses and on representative values estimated by Freeze and Cherry (1979) for the geologic materials present in the area.
Leakances as described by McDonald and Harbaugh (1988) are based on the layer thicknesses and the vertical hydraulic conductivity for each layer. The initial estimate of leakance between layers 1 and 2 was based on the thickness and assigned values for the hydraulic conductivity of the confining unit that separates the upper and lower aquifers. Vertical hydraulic conductivities were assigned based on the type of geologic materials present and consisted of zones representing till in the vicinity of Pine Mountain and Millie Hill, silt and clay in the south Kingsford area and to the north near Iron Mountain Well #4, and silt and silty sand over the rest of the area (fig. 17). These vertical hydraulic conductivities then were divided at each grid cell location by the layer thickness. The resulting values for leakance between layers 1 and 2 ranged from 3.8x10-6 to 2.5x10-3 /d (fig. 18). The initial estimate of leakance between layers 2 and 3 was based on the vertical hydraulic conductivity of the geologic materials present and the thickness of the over- and underlying layers. The resulting values for leakance between layers 2 and 3 ranged from 1.8x10-7 to 4x10-3 /d (fig. 19). Initial estimates of vertical hydraulic conductivity were assumed to be 0.01 times the horizontal hydraulic conductivity for each zone.

Figure 14. Hydraulic conductivity zones in the upper aquifer (model layer 1) for the ground-water-flow model, Kingsford-Iron Mountain area, Michigan.

Figure 15. Hydraulic conductivity zones in the lower aquifer (model layer 2) for the ground-water-flow model, Kingsford-Iron Mountain area, Michigan.

Figure 16. Hydraulic conductivity zones in the lodgement till and metasedimentary rocks (model layer 3) for the ground-water-flow model, Kingsford-Iron Mountain area, Michigan.

Figure 17. Hydraulic conductivity zones in theconfining unit between the upper and lower aquifers (model layers 1 and 2) Kingsford-Iron Mountain area, Michigan.

Figure 18. Estimated vertical leakance between the upper and lower aquifers (model layers 1 and 2), Kingsford-Iron Mountain area, Michigan.

Figure 19. Estimated vertical leakance between the lower aquifer and the lodgement till-metasedimentary rocks (model layers 2 and 3), Kingsford-Iron Mountain area, Michigan.

Recharge was assumed to differ in areas depending on the type and thickness of glacial deposits and was applied to the top model layer. In areas of surficial till, and in areas where the bedrock is located near the land surface, recharge was assumed to be zero. Over the rest of the area, recharge initially was assigned a value of 10 in/yr. Estimates of average discharge obtained from three nearby U.S. Geological Survey (USGS) gaging stations equaled 18.3 in/yr from 191to 1996 for the Menominee River near Florence, Wisc. (USGS gaging station 04063000); 15.1 in/yr from 1914 to 1997 for the Menominee River at Twin Falls near Iron Mountain, Mich. (USGS gaging station 04063500); and 13.8 in/yr from 1993 to 1997 for the Menominee River at Niagara, Wisc. (USGS gaging station 04065106, Water Resources Data for Michigan, Water Years 1996 and 1997). These discharges represent an upper limit for recharge because they include both baseflow and runoff. Actual estimates of baseflow could not be determined because of the presence of numerous control structures along the Menominee River.
In MODFLOW (McDonald and Harbaugh, 1988), streambed conductance is calculated as the product of the hydraulic conductivity of the streambed materials, stream length, and stream width, divided by the streambed thickness. Stream lengths for cells representing the Menominee River were equal to the length of the cell (fig. 8). Stream widths were assumed to be equal to approximately 100 ft. Stream lengths for cells representing the Pine River, Pine Creek, and Little Popple River in Wisconsin were equal to the length of the cell. Stream widths were assumed to equal 10 ft except along one section of the Pine River where the width was assumed to equal 100 ft. Stream lengths and widths for cells representing Antoine Creek were assumed to be 10 ft. Streambed thicknesses were assumed to be 1 ft. Hydraulic conductivities of the streambed materials were initially assigned a value of 10 ft/d, except for the section of Antoine Creek in the area immediately north of Iron Mountain Well #3 which was assigned an initial value of 1 ft/d. In this area, Coleman Engineering Company (1980) observed that the creek bottom consisted of silty sand which limits the direct flow of water from the creek to the well at this location. Horizontal hydraulic conductivities for the lake cells were initially estimated to be 10 ft/d.
Ground-water withdrawals were simulated from the center of the cell containing a pumping well. Because of the close proximity of Iron Mountain wells 1 and 2, these wells were simulated as pumping from the same cell. Withdrawals from Kingsford and Iron Mountain wells totaled:

Model Calibration

Model calibration is the process of reducing the difference between measured and simulated water levels and flows by adjusting model parameters. Multiple calibration trials were used to determine the set of parameter values that produced the best fit between simulated and measured values, whereas keeping parameters and the overall water budget within estimated values.
Head data used in the calibration of the ground-water-flow model include water-level measurements made during this study, measurements made during previous investigations of the Kingsford area, and measurements recorded on domestic well logs. Of the measurements used for this study, 24 are in layer 1, 44 are in layer 2, and 5 are in layer 3. Long-term drawdown measurements for the public-supply wells were included for comparison purposes. Because of the way in which pumping is represented in the model (continuous pumping at a constant rate, rather than intermittent pumping at varying rates); however, the water-level measurements at pumping wells will not be represented exactly.
Flow data at USGS gaging stations were examined to determine upper and lower bounds on flow within the Menominee River. Long-term data (greater than 80 years of record) are available for two stations (04063000 and 04063500) on the Menominee River upstream of the model area. Flow at these stations averaged 1,797 ft3/s for 04063000 and 1,812 ft3/s for 04063500 over the available years of record for similar size drainage basins (Water Resources Data for Michigan, Water Years 1996 and 1997). The length of river reach between these two stations is approximately the same length as that in the model area and the flow is assumed to be representative. The difference between these flows is small and is probably within the expected error range of the measurements; however, this comparison serves to illustrate that flow to river cells within the model is expected to be small (less than 15 ft3/s).
Model sensitivity was investigated by varying model parameters within expected ranges (scenarios 1 - 9, Appendix) and comparing the resulting root mean square errors (RMSE). Different scenarios were considered because of the complexity of the ground-water system being represented, differing geologic materials and water levels observed at adjacent wells, the steep gradients present, and the unknown location of the water table in some areas of the model. The root mean square error is the average of the squared differences in measured and simulated water levels. RMSE values and simulated model flow budgets were compared for scenarios in which recharge to ground water; hydraulic conductivity of the streambed materials, zone 1 in layer 1, Cowboy Lake area, Iron Mountain Well #4 area, and zone 5 in layer 2, and all zones in layer 3; and the leakances between layers were varied. The best fit between measured and simulated water levels was achieved when the Menominee River was represented as specified head cells in both layers 1 and 2 (scenario 1); however, this representation of the river does not allow for the possibility of flow beneath the river which was investigated in other scenarios (appendix).
The final parameter estimates (scenario 4) for each zone are within reasonable ranges (table 3).

The horizontal hydraulic conductivities of layer zones are within the ranges determined from aquifer tests and flow to river cells is within the expected order of magnitude. Flow to river cells increased and water levels in layer 1 were higher with an increase in recharge. Recharge of 2 in/yr to till areas did not change the extent of the unsaturated areas and water levels remained about the same in all layers. Water levels in layer 1 were lower with increases in hydraulic conductivity for layer 1 or 2. Water levels in layer 2 remained about the same with an increase in the hydraulic conductivity of layer 1 and decreased with an increase in the hydraulic conductivity of layer 2. Water levels remained about the same with changes in the hydraulic conductivity of layer 3. Water levels were higher in layer 1 and remained about the same in layer 2 with changes in streambed hydraulic conductivity; however, flow to river cells decreased with lower streambed conductivity and increased with higher streambed conductivity. Changes in leakance resulted in small changes in water levels.
Results of model simulations indicate close agreement between water flowing into and out of the model area; percent discrepancy equalled -0.49 (table 4). Generally simulated water levels are close to measured water levels with the residuals about evenly split above and below the line representing simulated equal to measured water levels (fig. 20a, b, and c). Seventy-one percent of the simulated water levels are within 15 ft of the measured values. Differences between simulated and measured values are shown in figure 21 and do not indicate any spatial trends. Simulated water levels for layers 1 and 2 are shown in figures 22 and 23, respectively. Part of the difficulty in matching measured values is due to large observed differences in nearby wells. The simulated streamflow of 9.6 ft3/s is less than 15 ft3/s as expected. The simulated values obtained represent the best match to these different values.

Figure 20a. Relation between observed and simulated water levels in the upper aquifer (model layer 1), Kingsford-Iron Mountain area, Michigan.

Figure 20b. Relation between observed and simulated water levels in the lower aquifer (model layer 2), Kingsford-Iron Mountain area, Michigan.

Figure 20c. Relation between observed and simulated water levels in the metasedimentary rocks (model layer 3), Kingsford-Iron Mountain area, Michigan.

Investigation of Boundary Effects

Choice of hydraulic boundaries (such as the specified heads used in this model) does not allow conditions along the boundaries to change from initial values; thus, one must investigate whether the boundaries are affecting the solution in the interior portion of the model (Anderson and Woessner, 1992). After selection of a final model, changes in conditions at boundaries were investigated by comparing flow through specified head cells forming boundaries for simulations in which wells were pumping and in which wells were idle. Flow through the individual specified head cells when the public-supply wells were pumping differed by an average of 2.6 percent (standard deviation equaled 23.8) from the flow through these same cells when the public-supply wells were simulated as idle. As an additional check on the boundaries, specified head cells were changed to specified flow cells using the well package to represent the specified flows. The ground-water-flow solution differed slightly near the edge of the dry cell area when the specified heads were changed to specified flows. Water-level differences between the two solutions at measurement locations were within 2 ft or less. Stresses induced from pumping by the public-supply wells only affected water levels locally and did not affect flow across streamlines.

Model Limitations

The ground-water-flow model was developed to simulate the regional steady-state response of the flow system in the Kingsford-Iron Mountain area to ground-water withdrawals by the public-supply wells. Hydraulic properties in the aquifers were assumed to be isotropic. (Within a cell hydraulic properties are the same in the north-south direction as in the east-west direction; hydraulic properties vary from location to location). Vertical variations in aquifer properties within layers and any variations in head or flow within the aquifers are not represented in the model. Model simulations are restricted to steady-state conditions. All stresses within and inputs to the system, including well withdrawals and recharge rates, remain constant throughout the simulation. No net gain or loss of flow is simulated in the system and no changes in ground-water storage occur. This model, in its current (2000) form, cannot be used to simulate transient-flow conditions.
Each grid cell represents the average hydrologic and hydraulic properties in the volume of aquifer represented by the cell, and any variations in properties within the volume represented by the grid cell cannot be represented. Local flows over distances smaller than the dimensions of the grid cell also cannot be accurately represented. Additional geologic and hydrologic data, as well as finer discretization of the model, would be needed to simulate smaller-scale flow systems. Hydraulic properties of the individual layers are represented by zones and local variations are not represented. This representation by zones generalizes the glacial and bedrock properties and small-scale differences in ground-water flows are not represented. This generalization was necessary because nearby wells often penetrate different lithologic units with varying thicknesses and detailed geologic information is unavailable for large portions of the model area. The accuracy of layer surfaces and hydraulic conductivity estimates are limited by the available data at well and boring locations. Additional control and accuracy could be achieved by inclusion of more data points.
Recharge to ground water was assumed to vary over the model area depending on the presence or absence of surficial till deposits and poorly permeable bedrock; thus, local variations in recharge rates, such as those associated with impermeable surfaces or differences in other types of surficial materials, are not represented in the model. Recharge in areas of surficial till deposits is assumed to be negligible. Simulated well withdrawals are assumed to come from the centers of grid cells. Small withdrawals from domestic wells were not included because of the difficulty in obtaining reliable data and the limitations in representing small-scale flow systems. However, domestic ground-water withdrawals are probably small at the scale of the model.
Other than flow within the Precambrian metasedimentary bedrock units in the vicinity of Pine Mountain and Millie Hill, the model does not simulate flow within the bedrock units. Thus, the bottom of the model is considered to be impermeable. This assumption is considered adequate for model development because of the limited flow available from fractures in bedrock. The model does not simulate flow through fractures or fissures in the rocks. Primary and secondary porosities are combined to represent average flow conditions. Thus, the presence of fractures and fissures could create local heterogeneities in hydraulic properties, which can not be represented by the model.
External boundary conditions were based on natural hydrologic conditions and were located distant from Kingsford's and Iron Mountain's well fields, and, based on the investigation of flow through the specified head cells, are assumed to have minimal effect on the solution in the interior portion of the model. Enlargement of the model area to natural physical hydrologic boundaries might improve model simulations and convergence and improve the match between simulated and observed ground water flows.
Different scenarios (in which recharge, hydraulic conductivity, and leakance parameters were varied) were considered because of the complexity of the system, different lithologic units and water levels observed at nearby wells, the steep gradients present, and the unknown location of the water table in some parts of the model area. Cells that went dry during model simulations were permitted to rewet; however, most remained dry. Most of the dry cells were in the vicinity of Pine Mountain, Millie Hill, and Trader Hill. In these areas the glacial deposits are thin and unsaturated. The measured water level in a well drilled on top of Pine Mountain shows that glacial deposits and bedrock are dry to the base of the mountain.
The accuracy of particle-tracking simulations is limited by the accuracy of the numerical model on which the simulations are based, the estimates of the effective porosity of the flow system, and the approximation of the cell flow velocities to the local ground-water-flow velocities. Additionally, the particle-tracking program considers ground-water flow by advection only. If the effects of dispersion were included, the contributing areas could be larger. Because the model does not specifically describe flow through fractures, ground-water flow and travel times in areas where fractures are present may not be accurately represented. Particle tracking was performed for public-supply wells near the center of the model area. Detailed analysis of contributing areas for small localized areas is inappropriate because the model was not designed to represent local flow systems. The areas that are represented as being dry in model simulations may contribute water and contaminants, if present, to the public-supply wells.

DELINEATION OF CONTRIBUTING AREAS

The particle-tracking program MODPATH (Pollock, 1989) can be combined with MODFLOW-calculated flow in each cell to determine areas contributing water to public-supply wells. In this study, the particle-tracking program was used by specifying the locations of hypothetical water particles along the faces of each cell containing a public-supply well. These particles then were tracked backward in time through the ground-water-flow system until they reached the top cell face in layer 1, which represents the starting location of the particle. This collection of starting locations, projected up to the land surface, represents the contributing areas for the public-supply wells. The subsurface zones through which the particles travel are known as the zone of transport. Particles also can be tracked for a specified amount of time, such as 10 years. This 10 year time-of-travel area is identified by those particles whose travel times are 10 years or less. An estimated effective porosity of 15 percent was used in the delineation for all layers. This value is less than the estimated total porosity of 25 percent, which is representative of these types of materials (Freeze and Cherry, 1996) and provides a conservative estimate of the contributing areas. Particle tracking results were compared for several of the scenarios and found to delineate similar areas.
Particle tracking describes the advective movement of ground water and does not incorporate the effects of adsorption, diffusion, dispersion, and degradation. Therefore, particle tracking is not intended as a substitute for modeling transport of dissolved chemicals in a ground-water system.
The 10-year contributing areas for Kingsford's public-supply wells encompass about 0.11 mi2, or 302 model cells, and consist of elongated areas to the east of the well fields (fig. 24). The 10-year contributing areas for Iron Mountain's public-supply wells encompass about 0.09 mi2, or 30 model cells, and consist of elongated areas to the east of the well field. The 10-year zone-of-transport consist of elongated zones east and west of the well fields. Actual contributing areas may be larger because the areas that are represented as being dry in model simulations may contribute water and contaminants, if present, to the public-supply wells. Forty year time-of-travel contributing areas are shown in figure 25.

Figure 21. Distribution of simulated minus observed water levels in the Kingsford-Iron Mountain area, Michigan.

Figure 22. Simulated ground-water levels in the upper aquifer (model layer 1), Kingsford-Iron Mountain area, Michigan.

Figure 23. Simulated ground-water levels in the lower aquifer (model layer 2), Kingsford-Iron Mountain area, Michigan.

Figure 24. Estimated 10-year contributing and zone-of-transport areas projected to the land surface in the Kingsford-Iron Mountain area, Michigan.

Figure 25. Estimated 40-year contributing and zone-of-transport areas projected to the land surface in the Kingsford-Iron Mountain area, Michigan.

SUMMARY

A ground-water-flow model was developed to simulate the regional, steady-state ground-water flow system in the Kingsford-Iron Mountain area. Kingsford and Iron Mountain are located in Dickinson County in the Upper Peninsula of Michigan. These two cities rely entirely on ground water withdrawn for public supply from the glacial deposits. The glacial deposits overlie Precambrian metasedimentary (Michigamme Slate, Vulcan Iron Formation, and Randville Dolomite) and metavolcanic (Badwater Greenstone) rocks, except along the Menominee River and in isolated areas where bedrock is exposed. Glacial deposits range from 0 to about 380 ft in thickness.
Two glacial sand-and-gravel aquifers are present in the study area. The upper aquifer ranges from 0 to 160 ft in thickness and consists of medium to fine grained sand with varying amounts of silt and gravel. The confining unit ranges from 0 to 150 ft in thickness and consists of well sorted and well compacted fine silt with clay. In some areas the silt and clay are absent due to partial burial and protrusion of stranded ice blocks that prevented silt deposition. The intermediate silt layer is a fairly effective barrier to hydraulic communication between the two aquifers. However, where the silt is absent in the Kingsford area, the two aquifers appear to be hydraulically connected. The lower aquifer ranges from 0 to 237 ft in thickness and consists of moderately well sorted fine to coarse sand with layers and lenses of silt. The lodgement till unit ranges from 0 to 25 ft in thickness.
Ground-water flow within glacial deposits in the Kingsford-Iron Mountain area is primarily towards surface-water bodies to the west and southwest and follows the regional topographic slope and bedrock surface. Ground-water flow is probably confined primarily within the permeable sand-and-gravel units because the Precambrian bedrock units appear to be largely unfractured and limit the downward migration of ground water. However, bedrock consists of metasedimentary bedrock units in the vicinity of Pine Mountain, and ground-water flow is locally affected by these units where they are fractured.
The ground-water-flow model developed for the Kingsford-Iron Mountain area encompasses part of Dickinson County, Michigan and Florence County, Wisconsin. This area is horizontally discretized into a variably spaced grid of cells in 213 columns and 204 rows. The model area is vertically discretized into three model layers. Layer 1 represents the upper aquifer in the glacial deposits; layer 2 represents the lower aquifer in the glacial deposits; and layer 3 represents either lodgement till or metasedimentary bedrock units where they are present in the model area. The intermediate confining unit between the upper and lower glacial aquifers is included only in the leakance term between layers 1 and 2. Specified-head, specified flux, and head-dependent flux cells form the model boundaries.
Conceptualized estimates of hydraulic conductivity for each layer were based on aquifer tests and on representative values taken from Freeze and Cherry (1979). These estimates were modified during model calibration. Each of the three layers was divided into zones containing similar geologic materials based on interpretation of well and geophysical logs.
The model was calibrated using multiple trials where parameters were changed to reduce differences between simulated and observed values, and to keep parameter estimates and the overall model budget within reasonable limits. Differences between observed and simulated values are attributed to the complex geology and hydrology of this area, the steep gradients in the vicinities of Pine Mountain, Millie Hill, and Trader Hill, and the large areas in which detailed hydrogeologic information was unavailable. Different scenarios (in which recharge, hydraulic conductivity, and leakance parameters were varied) were simulated in order to verify the reproducibility of regional ground-water-flow directions and contributing areas.
Contributing areas were delineated for Kingsford's and Iron Mountain's public-supply wells using particle-tracking analysis. Results of flow simulations based on 1997-1998 pumping conditions indicate that 10-year contributing areas, delineated by those particles with travel times of 10 years or less, cover 0.2 mi2.

REFERENCES CITED

APPENDIX

Citation:

Back to USGS, WRD Michigan Home Page