U.S. Department of the Interior

U.S. Geological Survey

Scientific Investigation Report 2004-5159

Simulation of Ground-Water Flow in 

Glaciofluvial Aquifers in the Grand Rapids 

Area, Minnesota
Simulation of Ground-Water Flow in 

Glaciofluvial Aquifers in the Grand Rapids 

Area, Minnesota

By Perry M. Jones

Prepared in cooperation with the Minnesota Department of Health

Scientific Investigation Report 2004-5159

U.S. Department of the Interior

U.S. Geological Survey
U.S. Department of the Interior

Gale A. Norton, Secretary

U.S. Geological Survey

Charles G. Groat, Director

U.S. Geological Survey, Reston, Virginia: 2004

For sale by U.S. Geological Survey, Information Services

Box 25286, Denver Federal Center

Denver, CO 80225

For more information about the USGS and its products:

Telephone: 1-888-ASK-USGS

World Wide Web: http://www.usgs.gov/

Any use of trade, product, or firm names in this publication is for descriptive purposes only and does not imply 

endorsement by the U.S. Government.

Although this report is in the public domain, permission must be secured from the individual copyright owners to 

reproduce any copyrighted materials contained within this report.

Jones, P.M., 2004, Simulatioin of ground-water flow in glaciofluvial aquifers in the Grand Rapids area, Minnesota: 

U.S. Geological Survey Scientific Investigations Report 2004-5159, 23 p.
iii

Contents

Abstract .......................................................................................................................................................... 1

Introduction ................................................................................................................................................... 1

Physical setting of the Grand Rapids area ...................................................................................... 3

Hydrogeology ....................................................................................................................................... 3

Simulation of ground-water flow ............................................................................................................... 5

Ground-water-flow model of the Grand Rapids study area .......................................................... 5

Data used ..................................................................................................................................... 8

Model description and assumptions ....................................................................................... 8

Model calibration ............................................................................................................................... 13

Simulation results .............................................................................................................................. 13

Model sensitivity to hydraulic conductivity, areal recharge, river-bed conductance,

and general-head-boundary conductance ...................................................................... 15

Model limitations and accuracy ...................................................................................................... 21

Summary ...................................................................................................................................................... 22

Acknowledgments ...................................................................................................................................... 23

References ................................................................................................................................................... 23

Figures

Figure 1. Location of modeled area, mining features, wetlands, public supply wells, 

gaging stations, lakes, and rivers in the study area. ............................................................ 2

Figure 2. Geology, geologic sections, locations of wells, and reaches of Mississippi and 

Prairie Rivers where base flow calibrations Model Layer five, and (f) Model 

Layer six—Lower aquiferwere done in (a) Model Layer one, (b) Model Layer 

two—Upper aquifer, (c) Model Layer three, (d) Model Layer four—Middle 

aquifer, (e) Model Layer five, and (f) Model Layer six—Lower aquifer .............................. 6

Figure 3. Model segmentation, showing model boundaries, and locations of ground-water 

withdrawals, for the MODFLOW simulation of ground-water flow in glaciofluvial 

aquifers in the Grand Rapids area, Minnesota ..................................................................... 12

Figure 4. Measured and simulated water-level altitude in wells (a) completed in upper, 

middle, and lower aquifers, (b) completed in the upper aquifer, (c) completed 

in the middle aquifer, and (d) completed in the lower aquifer for the calibrated 

MODFLOW simulation of ground-water flow in glaciofluvial aquifers in the 

Grand Rapids area, Minnesota ................................................................................................ 14

Figure 5. Difference between simulated and measured water-level contours for groundwater flow for the (a) Model Layer two—Upper aquifer, (b) Model Layer four—

Middle aquifer, and (c) Model Layer six—Lower aquifer, in the Grand Rapids 

area, Minnesota ........................................................................................................................ 16

Figure 6. Simulated potentiometric surfaces for ground-water flow and ground-water 

flow direction for the (a) Model Layer two—Upper aquifer, (b) Model Layer four—

Middle aquifer, and (c) Model Layer six—Lower aquifer, in the Grand Rapids 

area, Minnesota ......................................................................................................................... 18

Figure 7. Model sensitivity to selected input variables based on (a) root mean square error 

between simulated and measured ground-water levels, and (b) percentage error 

between simulated and estimated base flow for the calibrated MODFLOW

simulation of ground-water flow in glaciofluvial aquifers in the Grand Rapids 

area, Minnesota ......................................................................................................................... 22
iv

 Tables

Table 1. Geologic and hydrogeologic units in the Grand Rapids area, Minnesota and 

their water-bearing characteristics ......................................................................................... 4

Table 2. Data sources used in the MODFLOW steady-state simulation of ground-water 

flow in glaciofluvial aquifers in the Grand Rapids area, Minnesota ................................... 9

Table 3. Ground-water recharge and hydraulic conductivity values for the best-fit 

calibration of the MODFLOW steady-state simulation of ground-water flow in 

glaciofluvial aquifers in the Grand Rapids area, Minnesota .............................................. 10

Table 4. Regression statistics for measured and simulated water levels in wells and 

simulated and estimated base flow to reaches of the Mississippi and Prairie 

Rivers for the best-fit calibration of the MODFLOW steady-state simulation of 

ground-water flow in glaciofluvial aquifers in the Grand Rapids area, Minnesota ....... 14

Table 5. Relation among model variables used in sensitivity analysis for the calibrated 

MODFLOW simulation of ground-water flow in the upper, middle, and lower 

aquifers in the Grand Rapids area, Minnesota ..................................................................... 20

Conversion Factors, Abbreviations, and Vertical and 

Horizontal Datums, 

Multiply By To obtain

 

inch (in.) 2.54 centimeter (cm)

foot (ft) 0.3048 meter

mile (mi) 1.609 kilometer

square mile (mi2) 2.590 square kilometer 

foot per day (ft/day) 0.3048 meter per day

square foot per day (ft2/day) 0.09290 square meter per day

cubic foot (ft3) 0.02832 cubic meter

gallon (gal) 0.003785 cubic meter

cubic foot per second (ft3/s) 0.02832 cubic meter per second

degrees Fahrenheit  ̊C = (°F – 32) / 1.8 degrees Celsius

Vertical coordinate information to the North American Vertical Datum of 1988 (NAVD88): 

horizontal coordinate informatioin is referred to the North American Datum of 1927 (NAVD27)

Water year: The water year is October 1 through September 30 and is named for the calendar 

year in which it ends. 
ABSTRACT

A calibrated steady-state, finite-difference, ground-waterflow model was constructed to simulate ground-water flow in 

three glaciofluvial aquifers, defined in this report as the upper, 

middle, and lower aquifers, in an area of about 114 mi2

 surrounding the city of Grand Rapids in north-central Minnesota. 

The calibrated model will be used by Minnesota Department 

of Health and communities in the Grand Rapids area in the 

development of wellhead protection plans for their water 

supplies. The model was calibrated through comparison of 

simulated ground-water levels to measured static water levels 

in 351 wells, and comparison of simulated base-flow rates to 

estimated base-flow rates for reaches of the Mississippi and 

Prairie Rivers. Model statistics indicate that the model tends 

to overestimate ground-water levels. The root mean square 

errors ranged from +12.83 ft in wells completed in the upper 

aquifer to +19.10 ft in wells completed in the middle aquifer. 

Mean absolute differences between simulated and measured 

water levels ranged from +4.43 ft for wells completed in the 

upper aquifer to +9.25 ft for wells completed in the middle 

aquifer. Mean algebraic differences ranged from +9.35 ft for 

wells completed in the upper aquifer to +14.44 ft for wells 

completed in the middle aquifer, with the positive differences 

indicating that the simulated water levels were higher than the 

measured water levels. Percentage errors between simulated 

and estimated base-flow rates for the three monitored reaches 

all were less than 10 percent, indicating good agreement. 

Simulated ground-water levels were most sensitive to changes 

in general-head boundary conductance, indicating that this 

characteristic is the predominant model input variable controlling steady-state water-level conditions. Simulated groundwater flow to stream reaches was most sensitive to changes in 

horizontal hydraulic conductivity, indicating that this characteristic is the predominant model input variable controlling 

steady-state flow conditions.

INTRODUCTION

About two-thirds of the population of Minnesota uses 

ground water for their drinking water. About 2,600 community water-supply wells are present in Minnesota, with about 

one-half of these wells potentially vulnerable to contamination 

(Minnesota Department of Health, 2003).

Public water supplies in Minnesota are defined as sources 

of water to at least 15 service connections, or as sources to 26 

or more people per day for at least 60 days per year (Minnesota Department of Health, 2003). Public water supplies are 

divided into two categories: community and non-community 

water supplies. Community water supplies serve residents 

year-round, and include municipalities, mobile home parks, 

and apartment complexes. Non-community water supplies 

provide water to people in places other than their homes, such 

as hotels, restaurants, campgrounds, and schools. Non-community water supplies consist of transient and non-transient 

systems, where transient systems serve the traveling or transient public, and non-transient systems serve at least 25 nontransient individuals (Minnesota Department of Health, 2003). 

Communities in Minnesota and throughout the United 

States are being required to develop and implement wellhead 

protection plans to safeguard their public water supplies from 

contamination. Few tools are currently (2004) available to 

communities to effectively delineate drinking-water-supply 

management (or protection) areas surrounding their public 

supply wells. 

The Minnesota Department of Health (MDH) is working to assist communities to develop and implement wellhead 

protection plans through delineation of drinking-water-supply 

management (or protection) areas surrounding public supply 

wells. MDH hydrologists are developing tools and techniques 

for determining capture zones around public supply wells. 

One tool is a set of regional ground-water flow models in areas 

where groups of community water supply wells are located. 

The Grand Rapids area in north-central Minnesota is one of 

the areas where a group of public wells are located (fig. 1). 

Eighteen community supply wells and seven non-community, 

non-transient wells are present in a 45-mi2 area around Grand 

Rapids. Most of the water supplied to these wells is obtained 

from glaciofluvial aquifers. Glaciofluvial aquifers consist of 

boulders, gravels, sands, and silts deposited by running waters 

from melting glaciers. In the Grand Rapids area, these aquifers lie between glacial till units.

The U.S. Geological Survey (USGS), in cooperation 

with the MDH, conducted a 5-year study to simulate groundwater-flow conditions in aquifers in Minnesota. As part of 

this 5-year study, a 2-year study (2002-2003) was conducted 

Simulation of Ground-Water Flow in Glaciofluvial Aquifers 

in the Grand Rapids Area, Minnesota

By Perry M. Jones
2 Simulation of ground-water flow in glaciofluvial aquifers in the Grand Rapids area, Minnesota

Figure 1. Location of modeled area, mining features, wetlands, public supply wells, gaging stations, lakes, and rivers in the study area. 

�� ��

� �

�� ��

� �

�� ��

� �

�� ��

� �

�� ��

� �

�� ��

� �

�� �� �� �� �� ��

�

��

�

�

��

�

�

��

�

������ ������ ������ ������

��������

����

����������

����

����

����������

�������

����

����

����

�����������

����

����������

�������

���

��������

����������

����

��� �����

����

����������

����������

�������

����

������

����

����

������

����

����

�����������

���������� ����

������������� ����

���������

���������

�����������������

��������

��������

����������

�������� �����������

�����

�������

�����

����������� �����

�����������

�����

�����

������

��������

����������

��������

�����

������ ������ ������ ������

������

������

������

������������������������������������������������

���������������������������

����������������������

� � � � � �������

� � � � � ������������

������������

����������

������������

���������

����������������

�������

��������

����������

��������

������������������������������

�������������������������������������

����������������������������������������

��������������������������������������

���������������

�����������

B2

�����������������������

�����������������������

������������������

���������������

��������������

�������������

������������

� �

2

2

� �

�
to simulate ground-water flow conditions in glaciofluvial 

aquifers in the Grand Rapids area (fig. 1). A steady-state, 

regional, ground-water-flow model of glaciofluvial aquifers 

was constructed to simulate ground-water flow in the area. 

This regional ground-water flow model will be used by MDH 

and communities in the Grand Rapids area in the development 

of drinking-water-supply management areas surrounding public supply wells. This report presents the data sets, calibration 

results, and sensitivity analyses from the ground-water-flow 

model of the Grand Rapids area. Available geologic and 

hydrologic data collected from 1975 through 2001 were used 

for model construction and calibration. Data from scientific 

literature and other sources also were used.

Physical Setting of the Grand Rapids Area

The study area is about 114 mi2 surrounding the city 

of Grand Rapids (population 7,764 in 2000) (U.S. Census 

Bureau, 2002) in Itasca County, north-central Minnesota (fig. 

1). Cohasset and La Prairie are the only other cities in the 

study area. The watershed of the study area is part of the 

upper Mississippi River Basin, near the headwaters. Land-surface altitude in the study area ranges from 1,240 ft along the 

Mississippi River in the southeast portion of the study area to 

1,480 ft east of Pokegama Lake and south of Grand Rapids. 

Climate in the watershed is continental, with cold winters 

and hot summers. The mean annual temperature (1961-90) 

at Grand Rapids, Minnesota is 38.6°F, and the mean annual 

precipitation is 27.54 in. (Minnesota State Climatologist, 

2001). January is the coldest month, and July is the warmest 

month. February is the driest month and June is the wettest 

month. Mean January temperature is 3.8°F, and mean July 

temperature is 67.4°F. Mean February precipitation is 0.54 in., 

and mean June precipitation is 4.11 in.

Vegetation in the study area consists of a mix of northern 

hardwood forest and grasslands. Thirty-seven percent of the 

study area is forest, which consists of 18 percent deciduous 

forest, 11 percent mixed-wood forest, 5 percent young forest, 

and 3 percent coniferous forest (Minnesota Department of 

Natural Resources, 2002a). Grasslands cover 25 percent of 

the study area (Minnesota Department of Natural Resources, 

2002a).

Open water and wetlands cover 29 percent of the study 

area (Minnesota Department of Natural Resources, 2002a). 

Farmsteads, cultivated lands, and other rural development 

cover 3 percent of the study area, whereas 4 percent of the 

land cover consists of urban and industrial areas, mainly in the 

cities of Grand Rapids and Cohasset (Minnesota Department 

of Natural Resources, 2002a). 

Two percent of the study area consists of gravel pits, 

abandoned iron-ore open mines, and lands affected by mining activities (Minnesota Department of Natural Resources, 

2002a). Local topography and hydrology south of Lower 

Prairie Lake and west of Pokegema Lake have been affected 

by previous iron-ore mining activities at the Tioga No. 2 and 

Greenway Mines, and gravel mines (fig. 1). Abandoned open 

mine pits, tailings, stockpiles, and tailings and settling ponds 

are present in these areas (fig. 1) (Minnesota Department of 

Natural Resources, 2002a). No active iron-ore mining currently (2004) occurs in the study area. Tailings and stockpiles 

are as high as 150 ft above the land surface and extend over 

an area of several square miles. Because no active iron ore 

mining is present in the study area, the heights and extents of 

these piles are relatively static. Water levels in iron-ore mine 

pit lakes are relatively stable.

Surface drainage and ground-water flow through the 

study area is generally to the south-southeast (Oakes, 1970). 

Surface drainage north of the Mississippi River flows south to 

the river. The Mississippi River enters the study area from the 

west, flows east-southeast through a series of lakes (Cut-off, 

Blackwater, and Jay Gould Lakes) and the Paper Mill Reservoir, through the city of Grand Rapids, and southeast out of 

the study area (fig. 1). The Prairie River drains Lower Prairie 

Lake, located north of Grand Rapids, and flows south into the 

Mississippi River. Surface drainage south of the Mississippi 

River flows either north to the Mississippi River or south to 

Pokegama Lake (fig. 1). Pokegama Lake drains to the northwest through Little Jay Gould and Jay Gould Lakes to the 

Mississippi River. 

Hydrogeology

Glaciofluvial and bedrock aquifers are present in the 

study area (table 1). Three glaciofluvial aquifers are defined 

for the purposes of this report as the upper, middle, and lower 

aquifers. These aquifers are composed of glaciofluvial sediments that are separated or overlain by glacial tills or discontinuous deposits of glaciofluvial sediments. The primary 

bedrock aquifer in the study area is the Biwabik Iron Formation, which is composed of chert, shale, and iron minerals 

(table 1). The geologic units and water-bearing characteristics 

in the study area are presented in table 1.

Precambrian-age metavolcanic, sedimentary, and igneous 

bedrock, and Cretaceous-age sedimentary bedrock underlie 

glacial sediments in the study area (Morey, 1972). The Ely 

Greenstone of Archean age underlies glacial sediments in the 

northwestern part of the study area (fig. 2), and the Giants 

Range Granite of lower Precambrian age lie to the southsoutheast of the Ely Greenstone. A series of Precambrian-age 

and Cretaceous-age bedrock units lie above the Giants Range 

Granite in the central and southeastern part of the study area, 

striking to the east-northeast and dipping generally between 5° 

and 15° to the south-southeast (Oakes, 1970). The Pokegama 

Quartzite of Precambrian age overlies the Giants Range Granite in the central part of the study area (fig. 2). The Biwabik 

Iron Formation, composed of chert, shale, and iron minerals, 

overlies the Pokegama Quartzite and is overlain and bounded 

to the south by the Virginia Formation, consisting of argillites, siltstones, and graywackes (Morey, 1972) (fig. 2). The 

Biwabik Iron Formation is mined for iron ore and the main 

Introduction 3
4 Simulation of ground-water flow in glaciofluvial aquifers in the Grand Rapids area, Minnesota Table 1. Geologic and hydrogeologic units in the Grand Rapids area, Minnesota and their water-bearing characteristics [gal/min, gallons per minute] System Series Geologic unit Lithology Water-bearing characteristics Hydrogeologic unit in model Equivalent layer in the digital ground-water flow model Quaternary Pleistocene Glaciofluvial sediments fluvial sands, gravels, and boulders well yields more than 1,000 gal/min in middle and lower aquifers, upper, middle, and lower aquifers; glaciofluvial sediments separating aquifers layers 1, 2, 3, 4, 5, and 6 Quaternary Pleistocene Glacial tills clay, sands, silts, pebbles, cobbles, and boulders limited permeability glacial tills layers 1, 2, 3, 4, 5, and 6, basal (no-flow) boundary below lower aquifer Cretaceous Upper Coleraine Formation sandstones, iron-forma- tion, and shales limited, occurs in isolated patches above the Virginia Formation not represented in model basal (no-flow) boundary Precambrian Middle Virginia Formation argillite, siltstones, and graywackes limited permeability, well yields as much as 10 gal/mn not represented in model basal (no-flow) boundary Precambrian Middle Biwabik Iron-Formation ferruginous chert, shale, and iron minerals, also called taconite. well yields as much as 500 gal/min not represented in model basal (no-flow) boundary Precambrian Middle Pokegama Quartzite quartzite, dense, hard, conglomeratic at base not used as a source of water not represented in model basal (no-flow) boundary Precambrian Middle Giants Range Granite granite intrusive well yields are low (less than 10 gal/min), only used for domestic water supplies where glacial drift thick- nesses are small not represented in model basal (no-flow) boundary Precambrian Middle Ely Greenstone greenstone, green schist not used as a source of water not represented in model basal (no-flow) boundary
bedrock aquifer in the area. Isolated patches of sandstones, 

iron-formation, and shales of the Cretaceous-age Coleraine 

Formation overlie the Virginia Formation in the southeastern 

part of the study area. Bedrock valleys have been identified 

in the study area (Bruce A. Bloomgren, Minnesota Geological Survey, oral commun., April 17, 2001), but the extent and 

depth of these valleys is poorly defined by the small number 

and distribution of wells and boreholes. 

Glacial sediments cover the entire study area, with 

thicknesses ranging from less than 50 ft in the central and 

northwestern part of the study area to greater than 400 ft in 

the southeastern part of the study area (Oakes, 1970; Winter, 

1973a; Jirsa, and others, 2002). Three major morainal tills and 

associated glaciofluvial sediments exist, which were formed 

during the Wisconsin ice advances from the north and west 

of the study area (Winter, 1971). All three of the tills are 

discontinuous in the study area, with the extents and thicknesses varying over short distances, forming areas of hydrologic confinement. The stratigraphically lowest till, the basal 

till, is a dark-greenish and brownish-gray till that ranges from 

clayey to sandy in texture, and is calcareous (Winter, 1971). 

The middle boulder till ranges widely in color from gray to 

yellow, and consists of sands and silts, with abundant cobbles 

and boulders (Winter, 1971). This till tends to be thickest in 

the southeastern part of the study area, and greater than 100 

ft thick. The surficial till is brown in color; sandy, silty, and 

calcareous; and is generally less than 30 ft thick in the study 

area. Sand, silt, and clay deposits reworked and redeposited 

by Glacial Lake Aitkin overlie surficial till and glaciofluvial 

sediments southeast and northwest of Grand Rapids (Winter, 

1973b).

Glaciofluvial sediments of the three glaciofluvial aquifers 

lie stratigraphically above the surficial till, between surficial 

and boulder tills, and commonly lie between the boulder and 

basal till or bedrock (Winter, 1973a). These deposits consist 

largely of sands, gravels, and boulders. Glaciofluvial sediments above the surficial and between the surficial and boulder 

tills, the upper and middle aquifers, respectively, are the most 

continuous outwash deposits in the study area, commonly 

greater than 50 ft thick and sometimes greater than 100 ft thick 

in portions of buried valleys (Winter, 1973a). These sediments 

consist of fine-grained sands throughout much of the study 

area. However, outwash deposits can be highly transmissive 

where coarse-grained sands, gravels, and boulders were deposited in buried valleys and other locations where the bedrock 

surface was low. The glaciofluvial sediments below the boulder till, the lower aquifer, are fairly continuous in most of the 

study area, but are absent in the northwestern part of the study 

area. The glaciofluvial sediments below the boulder till are 

poorly sorted and are generally less than 50 ft thick, but are 

greater than 100 ft thick at some locations in buried bedrock 

valleys (Winter, 1973a). Some glaciofluvial sediments are in 

direct contact with glaciofluvial sediments of earlier or later 

glacial flow events.

Glaciofluvial aquifers and the Biwabik Iron Formation 

are the main sources of ground-water in the Grand Rapids 

area. Well yields from glaciofluvial aquifers can be more than 

1,000 gal/min where glaciofluvial sediment thickness and the 

horizontal extent are large (table 1). These high well yields 

are found in the middle or lower aquifers. The Biwabik Iron 

Formation ranges in thickness between 350 to 500 ft (Oakes, 

1970), producing water mainly within fractures and bedding 

planes. Wells opened to the Biwabik Iron Formation can 

produce yields as high as 500 gal/min (table 1). These wells 

generally are located near the subcrop of the formation in the 

central portion of the study area (fig. 2). The Virginia Formation and Giants Range Granite are used as a source of water 

for some domestic water supplies. However, well yields typically are low, generally less than 10 gal/min (table 1).

Ground-water withdrawals mainly are from glaciofluvial 

aquifers through municipal, small industrial, and domestic 

wells. In 2001, the city of Grand Rapids withdrew 342.4 

million gallons of water from four wells completed in 

glaciofluvial sands and gravels, 256.6 million gallons of water 

from a well completed in the Biwabik Iron Formation, and 

15.9 million gallons of water from Pokegama Lake (fig. 1) 

(Minnesota Department of Natural Resources, 2002b). The 

city of La Prairie (population 605 in 2000) (U.S. Census 

Bureau, 2002) purchases water from Grand Rapids. In 2001, 

the city of Cohasset (population 2,481 in 2000) (U.S. Census 

Bureau, 2002) withdrew 22.5 million gallons of water from 

two wells completed in buried glaciofluvial sands and gravels to serve 755 people (Minnesota Department of Natural 

Resources, 2002b). The remaining population of Cohasset 

withdrew water from private wells (Michael Stejskal, City of 

Cohasset, Minnesota, oral commun., March 4, 2004). In 2001, 

Minnesota Power withdrew 357.4 million gallons of water 

from three wells completed in glaciofluvial sands and gravels 

for generation of electricity from their coal-fire power plant 

in Cohasset (Minnesota Department of Natural Resources, 

2002b). Within the study area, 14 other public supply wells 

and 2 irrigation wells annually pumped about 39.6 million 

gallons of water from wells completed in glaciofluvial aquifers 

in 2001 (Minnesota Department of Natural Resources, 2002b). 

Domestic wells in the study area extract water mostly from 

glaciofluvial aquifers, with a few households using the aquifer 

in the Biwabik Iron Formation for a source of water. Only 

withdrawals from wells completed in glaciofluvial aquifers 

were simulated in the model described in this report.

SIMULATION OF GROUND-WATER 

FLOW

Ground-water flow in the Grand Rapids area was 

simulated using a calibrated, steady-state, finite-difference, 

ground-water-flow model. The model was calibrated using 

water levels in wells and estimated base flows in reaches of 

the Mississippi and Prairie Rivers. Once the steady-state simulation was calibrated, a series of sensitivity analyses were run 

with the model to determine the effect of changes in hydraulic 

Simulation of ground-water flow 5
6 Simulation of ground-water flow in glaciofluvial aquifers in the Grand Rapids area, Minnesota

������ ������ ������ ������

������

������

������

�������������������������������������

�����������������

���������

�������

���� ����������

����

����

��������������������������������������������������������

�������������������������������

����������������������������������������������������������

���������������������������������

���������������������������������������������������������

�������������������������������

����������������������������������������������

�������

�����

�������������������

�����������

�������

�����

����

B

B

B B

B

B

B BBB

B

B

B

B

BBB

B

B

B

B

B

BB

BBBB

B

BB

B

B

BB

B

B

B

B

B

B

BB

B

B B

B

B

BB

B

B

B B BB

BB

BB

B

B

BB

B

BB

B B

B B

BB

B

B

B B BB

B

BB

B BB

BB

BB

BBB

BB

B

B

B

B B

B

B

B

B

B

������ ������ ������ ������

������

������

������

�����������������������������������

������������

����������

�����������������

���������

�������

���� ����������

����

����

�����������

����

�����������������

�������������

�����������

�����

�������

�����

�������

�����

�����

�����������������

�������������

�����������

�����

�������

�����

�������

�����

�����

�����������

�������

�����

���� ������ ������ ������ ������

������

������

������

������������

����������

�����������������

���������

�������

���� ����������

����

����

�����������

����

�����������

�����

�������

����� �������

�����

�����

���������������������

�����������

�������

�����

����

B

B

B

B

B

BB

B

B BB B B B B BB B B

B B

BB B

BB B

B

BBB B B

B B

B B

B B BB B

B B B

B B B

B B B B B B B

B B B B BBB

BB B B B B BBB B B

B B B

B B B BBBB B BBB

B B B BB B

B

BB B

BB

B B B B BB

B BB B

B

B B

BBBBBBB

B

BB B

BB

B

BB

BBB B B B B

B B BBB B BBB B B BB B

B B B BB B

B

B

B B BB

B

B

B BB

BB BBB

B B B B

B

BB B B BBBB B B BB B BB

B B B BB BB BB BB BB B B BB B

B B BB B B B B

B

B

B

B

B

B

B

B

B

B

B

B

B

���� ��������

����������

�����������

����

�����������

�������

�����

�������� ����

����������

���� ����

���������� �������

����

���� ����

���������� ����

�����������������

������ ������ ������ ������

������

������

������

� ��

�����

���������������������������������������������������

��������������������������������������������

�����������������������������

� ��

��������������������

��������������

Pokegama Quartzite

�������� ��������������

���������������������������������

����������������������������������������

Figure 2. Geology, geologic sections, locations of wells, and reaches of Mississippi and Prairie Rivers where base flow calibrations 

Model Layer five, and (f) Model Layer six—Lower aquifer
� � � � � �������

� � � � � ������������

�������������������� ������������������������������������

�����������

�������

�����

���� ������ ������ ������ ������

������

������

������

������������

����������

�����������������

���������

�������

���� ����������

����

����

�����������

����

�����������

�������

�����

����

B

B

B

B B

B

B BB

B B B

B BB

BBB

B

B BB B

BBBBB

B B

BB

B

B B

BB BBBB

BBB

B

B

B

������ ������ ������ ������

������

������

������

������������

����������

�����������������

���������

�������

���� ����������

����

����

�����������

����

B

�����������

������������������������������������������

��������������������������������������������������������������

���������������������������������������������������������������

����������������������������������������������������������

������������������������������������������������������������

����������������������������������������������������������

����������������������������������������������������������

����������������������������������������������������������

���������������������������������������������������������

����������������������������������������������������������

����������������������������������������������������������

����������������������������������������������������������

����������������������������������������������������������

�������������������

�����������������������������

����

������������������

�������������������������������������������

����

�����

�����

�����

�����

�����

�����

�����

�����

�����

�

����

��

����

�����������������

�����������������

�������������

����������������

����������������������������

��������������������

���������

���������

�������

��������������

�������������������������������������������������������������������������

�

�

�

� �

�

������������

������� �����������������������������

��������������������������������

������

���������

��������

����������

�������

���������

�����

�����

�����

�����

�����

�����

�����

�����

����

� ��

����� ����� ���������� ������������� �����������������

��������������������

�������

������������������ ��������������

����������������

����������������������������

�����������������������������

��������������������������������

���������������

�������������������������������

�����������������

���������������������������������

����������������

�������������������������������

��������

���������

were done in (a) Model Layer one, (b) Model Layer two—Upper aquifer, (c) Model Layer three, (d) Model Layer four—Middle aquifer, (e) 

Simulation of ground-water flow 7
8 Simulation of ground-water flow in glaciofluvial aquifers in the Grand Rapids area, Minnesota

conductivity, areal recharge, and river-bed conductance on 

simulated water levels and streamflows in the model. 

Ground-Water-Flow Model of the Grand Rapids 

Study Area

A three-dimensional, finite-difference, ground-waterflow model was developed, incorporating an area of about 

114 mi2 surrounding the city of Grand Rapids (fig. 1). The 

USGS MODular ground-water-FLOW model, MODFLOW, 

was used to simulate ground-water flow conditions in the 

Grand Rapids area. MODFLOW ground-water-flow model 

that simulates steady-state and transient ground-water flow 

in multiple aquifers (Harbaugh and McDonald, 1996) and 

simultaneously solves a series of mathematical equations that 

describe saturated ground-water flow where Darcy’s Law 

applies. In MODFLOW, ground-water flow can be simulated in unconfined aquifers, confined aquifers, and confining 

units. A variety of hydrologic features and processes, such as 

rivers, streams, wells, evapotranspiration, and recharge from 

precipitation, also can be simulated. The steady-state representation of the study area was done using the BASIC, BCF, 

General-Head, River, Well, Recharge, and Preconditioned 

Conjugate Gradient Solver (PCG2) packages of MODFLOW. 

The MODFLOW simulation was developed, run, and analyzed 

using the Groundwater Modeling System (GMS) graphic-user 

interface (U.S. Department of Defense, 1998). Hydrologic and 

geologic data available in GIS coverages were imported into 

the model using GMS.

Simulation of ground-water flow conditions in the study 

area was done using the following four-step approach: (1) 

compile existing hydrologic and geologic data; (2) discretize 

the compiled data; (3) calibrate the model through the comparison of simulated and measured ground-water levels and 

stream base-flow rates; and (4) perform sensitivity analyses on 

the calibrated model to assess the effects of hydraulic conductivity, areal recharge, and river-bed conductance on simulated 

ground-water levels and base-flow rates.

Data Used

A variety of data sources were used to represent various hydrologic features in the model (table 2). The model 

properties required to represent these features include extent, 

thicknesses, and boundaries of aquifers; areal recharge rates, 

water-level altitudes, depths of surface-water bodies that affect 

ground-water flow, and withdrawal rates from wells completed 

in aquifers. The types of data sources consisted of geologic 

descriptions, Geographic Information System (GIS) data, 

water-level altitudes, and municipal well pumping records. 

A three-dimensional representation of the glacial 

hydrogeologic units in the study area was created by MDH 

hydrologists based on interpolation of borehole data and geologic maps. Geologic logs of existing municipal, domestic, 

and monitoring wells were used to develop representations of 

hydrologic units. Logs for wells were obtained from waterwell records in the Minnesota Geological Survey’s County 

Well Index System (table 2). Geologic descriptions on the 

logs were not detailed enough to distinguish between the surficial, middle boulder, and basal till units, so these units were 

simulated as a single hydrogeologic unit named glacial tills 

(table 1). The inverse-distance-weighting method was used to 

interpolate between well logs, using the 12 nearest neighbors, 

a power of 6, and outlined aquifer boundaries for barriers. The 

interpolation was refined by referencing geologic maps and 

publications by Oakes (1970) and Winter (1971, 1973a, and 

1973b), and bedrock-depth maps developed by the Minnesota 

Department of Natural Resources (MNDNR). 

A series of GIS coverages were used to identify and 

represent hydrologic features and processes in the study area 

(table 2). Data sets for perennial wetlands, lakes, and rivers were obtained from the National Wetlands Inventory data 

base (U.S. Fish & Wildlife Service, 1994). Included were 

natural and man-made lakes, such as tailings ponds, settling 

ponds, and mine pit lakes. Water-level altitude data for the 

perennial wetlands, rivers, and lakes used in the model were 

obtained directly from the National Wetlands Inventory data 

base and from the MNDNR Lake Level data base (Minnesota 

Department of Natural Resources, 2002c). A 98-ft (30-m) 

digital elevation model (DEM) based on the hypsography of 

1:24,000-scale USGS quadrangle maps was interpolated to 

develop a 328-ft (100-m) land-surface altitude grid used for 

the top of the upper layer of the model. 

Annual pumping records for 1997-2001 for public supply wells in the study area were obtained from the MNDNR 

water-appropriation-permit data base (Minnesota Department 

of Natural Resources, 2002b). The maximum annual withdrawal rates for 1997-2001 were used in the model to represent steady-state pumping rates from public supply wells.

Model Description and Assumptions

A three-dimensional, numerical model of ground-water 

flow was constructed based on a conceptual model of the 

hydrogeology in the study area. The conceptual model was 

created based on knowledge of the hydrogeologic setting, 

aquifer characteristics, distributions and amounts of groundwater recharge and discharge, and aquifer boundaries. No 

storage terms were included in the model. Static water levels 

measured between 1975 and 1999 and the maximum annual 

withdrawal rates for 1997-2001 were considered an acceptable 

estimate of the ground-water system at equilibrium. Withdrawals from individual wells vary annually, but the maximum 

annual withdrawal rate for 1997-2001 approximates the average withdrawals through the 1975-1999 water-level measurement period (James Walsh, Minnesota Department of Health, 

oral commun., March 7, 2003).

A “true-layer” approach was used to define layering 

represented by the model, explicitly defining altitudes and 

aquifer hydraulic properties of cells in each layer based on the 

three-dimensional geologic model representation. Using this 
approach, MODFLOW geometry-dependent parameters such 

as transmissivity and Vcont (“leakance”, see McDonald and 

Harbaugh, 1988) are automatically computed when the model 

files are saved. The study area was discretized into rectangular finite-difference grid cells within which the properties of 

the hydrogeologic unit represented were assumed to be uniform and isotropic. Ground-water flow was simulated within 

the aquifers using a block-centered approach, in which flow is 

calculated between discretized cells based on head conditions 

at the central nodes of the cells (McDonald and Harbaugh, 

1988). Hydrogeologic properties and stresses were applied 

assuming that the assigned properties and stresses represent 

average conditions within that cell. Starting hydraulic head 

values in the model were extrapolated from a water-level 

surface obtained from interpolation of measured water levels 

in wells. Those head values were later reset to head values 

calculated by early iterations of the model to improve model 

convergence.

The three-dimensional, finite-difference grid used in the 

model was evenly spaced, consisting of 239 rows and 263 

columns. The dimensions of the grid cells were 328.08 ft 

(100 m) along rows and columns. The model was vertically 

divided into six layers, based on the hydrogeologic units (fig. 

2, table 1). The layers were represented as either confined or 

unconfined, with their transmissivities varying with saturated 

thicknesses under unconfined conditions and being constant 

under confined conditions. Simulation of flow in the model 

was dependent on the thicknesses and hydraulic conductivities 

between adjacent cells and layers. A detailed discussion of 

the simulation of flow in the model can be found in McDonald 

and Harbaugh (1988).

The three principal glaciofluvial aquifers used in the 

study area for ground-water supplies and represented in the 

model are herein defined as the upper, middle, and lower 

aquifers. In the model, these aquifers are simulated in layer 

2 (upper aquifer), layer 4 (middle aquifer), and layer 6 (lower 

aquifer) (fig. 2). Glaciofluvial sediments and glacial tills 

Table 2. Data sources used in the MODFLOW steady-state simulation of ground-water flow in glaciofluvial aquifers in the Grand Rapids 

area, Minnesota

Hydrologic feature Types of data Sources of data Types of information Date compiled

Aquifer and model 

boundaries and geology

GIS data sets (shapefiles) and Digital 

Elevation Model 

(DEM)

U.S. Geological Survey, Minnesota Department of Health, 

which used Minnesota Geological Survey’s County Well Index 

record, Minnesota Department 

of Natural Resources bedrockdepth maps, and U.S. Geological Survey publications

location, geologic logs, 

geologic description, 

stratigraphy, land-surface altitude, bedrockdepth maps

February - July 

2002

Wetlands GIS data sets U.S. Fish & Wildlife Service, National Wetlands Inventory data 

base (1994)

location, area, water-level 

altitude

February 2002

Lakes GIS data sets U.S. Fish & Wildlife Service, National Wetlands Inventory data 

base (1994)

location, area, water-level 

altitude

February 2002

Rivers (polygons) GIS data sets U.S. Fish & Wildlife Service, National Wetlands Inventory data 

base (1994)

location, area, water-level 

altitude (stage)

February 2002

Rivers (Arc and Node) GIS data sets U.S. Fish & Wildlife Service, National Wetlands Inventory data 

base (1994)

location, area, water-level 

altitude (stage)

February 2002

Lake water-level altitude U.S. Fish & Wildlife Service, National Wetlands Inventory data 

base (1994), Minnesota Department of Natural Resources, 

Division of Waters, Lake Level 

data base (2002c)

altitude for lakes February and

 March 2002

Areal Recharge GIS data sets Developed from U.S. Geological 

Survey data

location, recharge rate June 2002

Public supply and other 

large-capacity wells

well pumping records Minnesota Department of Natural 

Resources, Division of Waters, 

Water-Appropriations-Permit 

data base (2002b)

annual pumping rates June 2002

Simulation of ground-water flow 9
10 Simulation of ground-water flow in glaciofluvial aquifers in the Grand Rapids area, Minnesota

in layers 1, 3, and 5 represent glaciofluvial sediments and 

glacial tills lying above or between, and hydraulically interconnected with the principal aquifers. Glacial clays and tills 

were present in parts of each model layer, and were simulated 

in the model as a single hydrogeologic unit; herein termed 

glacial tills. Hydraulic conductivities were specified for the 

glaciofluvial sediments in each layer and a single hydraulic 

conductivity value was specified for the glacial tills (table 

3). Initial hydraulic conductivity values for each of the 

hydrogeologic units were obtained from MODFLOW simulations in the Coleraine/Bovey area east of the Grand Rapids 

modeled area (Jones, 2002). In GMS, horizontal and vertical hydraulic conductivity values are entered directly into the 

program. Model layer elevations are used by GMS with the 

entered vertical hydraulic conductivity values to get the Vcont 

(“leakance”) values used in MODFLOW (McDonald and Harbaugh, 1988). Precambrian and Cretaceous bedrock, and basal 

till below the lower aquifer, represented the basal boundary of 

the model (table 1).

Perennial rivers, streams, lakes, and wetlands were specified as head-dependent boundaries (general-head or river) in 

layers 1-5 of the model (fig. 3). The model area was bounded 

to the north by Stevens Lake, Bass Lake, Shoal Lake, the Prairie Lake-Lower Prairie Lake System, and a series of wetlands 

(figs. 1 and 3). The eastern part of the model is bounded by 

the Prairie River, Mud Lake, Blackberry Lake, Mississippi 

River, and a series of wetlands and tributaries (fig. 1). Pokegama Lake, Rice Lake, and a series of small lakes, ponds, and 

wetlands bound the model to the south. The western boundary of the model consists of wetlands and two lakes, Little 

Drum and Little Rice Lakes. Cells located inside the model 

boundary were active, and cells outside of the model boundaries were inactive, where no flow into or out of the cells was 

simulated. 

 The bottom of the lowest layer in the model was simulated as a no-flow boundary. The volume of water moving 

downward across the bottom of the lowest layer in the model 

was assumed to be small relative to the amount of horizontal 

flow. With this assumption, ground-water withdrawals from 

the Biwabik Iron Formation were assumed to have minimal 

effect on steady-state water levels in the glaciofluvial aquifers.

A total of 14,314 general-head boundary cells, representing lakes and perennial wetlands, were specified in the 

model. Flow into and out of these cells was calculated by the 

model based on the difference between the head in the cell and 

the assigned general head, multiplied by a conductance term 

(McDonald and Harbaugh, 1988). In each of the general-head 

boundary cells, the conductance was defined as the hydraulic 

conductivity of the lake or wetland bed material divided by 

the vertical thickness of that material, multiplied by the area 

of the lake or wetland in that cell. Both hydraulic conductivity and vertical thickness of lake and wetland bed material 

can be highly variable for lakes and wetlands. The selected 

hydraulic conductivity value of 0.066 ft/day (0.021 m/day) 

represents a small value typical of the range of hydraulic 

conductivity values for glacial tills obtained from single-well 

hydraulic conductivity tests done by Jones (2002) and for 

hydraulic conductivity values determined for lake-bed material 

of nearby Shingobee Lake (Kishel and Gerla, 2002). Many of 

the lakes and wetlands in the study area lie above glacial till, 

so a hydraulic conductivity value for the glacial till is relevant 

for the lake-bed material. An assumed value of 3.28 ft (1.0 m) 

was used to represent an average thickness for the bed material 

for lakes and wetlands in the area, which was used by Jones 

(2002) for general-head boundaries. Lake and wetland areas 

and altitudes from the National Wetland Inventory GIS coverages were used for the areas and head values, respectively, for 

the general-head boundaries representing the lake and perennial wetland segments. 

Streams and rivers were simulated using the river package of MODFLOW (fig. 3). The river package was used to 

simulate flow between the surface-water features and groundwater systems (McDonald and Harbaugh, 1988). Streams 

and rivers were segmented into reaches based on underlying 

geology and altitudes of the streams and rivers for assigning hydrologic variables to model cells. These reaches were 

segmented into cells (fig. 3). A total of 1,840 river cells were 

simulated in the model. Flow between the reaches and the 

ground-water flow systems was calculated for each cell based 

on the head difference between the river and the aquifer, 

multiplied by a conductance term (McDonald and Harbaugh, 

1988). In the river package, the conductance term is defined 

as the hydraulic conductivity of the river-bed materials divided 

by the vertical thickness of the river-bed materials, multiplied 

by the surficial area of the river bed in the cell (McDonald and 

Harbaugh, 1988). Hydraulic conductivities of 0.066 ft/day 

(0.021 m/day) and 13 ft/day (4.0 m/day) were used for the 

conductance terms for river segments where the underlying 

geology was till or glaciofluvial sediments, respectively. A 

vertical thickness of 3.28 ft (1.0 m) was used for the thickness 

of each river segment. Jones (2002) also used these hydraulic 

conductivity and vertical thickness values in simulation of 

ground-water flow surrounding the nearby Canisteo Mine (fig. 

1). Surficial areas and altitudes for the river segments were 

obtained from the National Wetland Inventory GIS coverages. 

The altitude of the river bed was assumed to be 6.56 ft (2.0 m) 

below the altitudes for the river segments. Jones (2002) also 

used this assumption.

A specified-flux boundary was used to represent areal 

recharge to the upper layer (layer 1) of the model using the 

recharge package in MODFLOW. Areal recharge was represented as the net difference between precipitation and evapotranspiration losses occurring above the water table. Initial 

recharge rates were proportioned based on results from stream 

hydrograph analyses from Jones (2002) and the geology of 

the upper two layers of the model. The largest areal recharge 

rate (1.1 x 10-3 ft/day (3.4 x 10-4 m/day)) was assumed to occur 

where glaciofluvial sediments (sands, gravels and boulders) 

were present in layer 1 of the model and where glacial tills in 

model layer 1 were less than 3.28 ft (1.0 m) thick and overlied 

glaciofluvial sediments in layer 2 (fig. 2; table 3). The smallest recharge rate (4.6 x 10-6 ft/day (1.4 x 10-6 m/day)) occurred 
Simulation of ground-water flow 11

Table 3. Ground-water recharge and hydraulic conductivity values for the best-fit calibration of the MODFLOW steady-state simulation of ground-water flow in glaciofluvial aquifers in the Grand Rapids area, Minnesota[All values are in feet per day: ---, no value]

Geologic unit, model layer, and hydrogeologic units Ground-water recharge Horizontal hydraulic conductivity Vertical hydraulic conductivity

Initial Final Initial Final Initial Final

Glacial tills - Layer 1, where till thickness is equal to or 

greater than 5 meters 

4.6 x 10-6 4.6 x 10-6 6.6 x 10-2 2.0 x 10-2 6.6 x 10-2 2.0 x 10-3

Glacial tills - Layer 1, where till thickness is equal to or 

greater than 1 meter and less than 5 meters 

2.3 x 10-4 2.3 x 10-4 6.6 x 10-2 2.0 x 10-2 6.6 x 10-2 2.0 x 10-3

Glacial tills - Layer 1, where till thickness is less than 1 meter and glaciofluvial sediments are present in Layer 2 1.1 x 10-3 1.1 x 10-3 6.6 x 10-2 2.0 x 10-2 6.6 x 10-2 2.0 x 10-3Glaciofluvial sediments - Layer 1 1.1 x 10-3 1.1 x 10-3 13 160 13 16Glacial tills - Layers 2-6 ---- ---- 6.6 x 10-2 2.0 x 10-2 6.6 x 10-2 2.0 x 10-3Glaciofluvial sediments - Layer 2 (upper aquifer) ---- ---- 13 98, 100, 170, 

200, 220

 13 9.8, 10, 17, 20, 22

Glaciofluvial sediments - Layer 3 ---- ---- 13 98 13 9.8Glaciofluvial sediments - Layer 4 (middle aquifer) ---- ---- 13 98 13 9.8Glaciofluvial sediments - Layer 5 ---- ---- 13 98 13 9.8Glaciofluvial sediments - Layer 6 (lower aquifer) ---- ---- 13 110, 270 13 11, 27
12 Simulation of ground-water flow in glaciofluvial aquifers in the Grand Rapids area, Minnesota

Figure 3. Model segmentation, showing model boundaries, and locations of ground-water withdrawals, for the MODFLOW simulation of 

ground-water flow in glaciofluvial aquifers in the Grand Rapids area, Minnesota

��������

����

����������

����

�������

����

����

����������

�������

����

����

����

�����������

����

�����������������

������ ������ ������ ������

������

������

������

� �

� �

2

2 �

������������������������������������������������

���������������������������

���������������������� � � � � � �������

� � � � � ������������

�����������

�����

�������

�����

�����������

��

�����������

�����

�������

����� �������

�����

�����

�����������

�����

�������

�����

���������������

���������������������

������������

���������������������

�������������������

���������������������

�����������������������

���������������������

��������������������������

���������������������

�����������������������������

���������������

�����������������������

����������������������������

�����������

���������������������

�������������������

���������������������������������������������

����������������������������������������

���������������������������������������������

����������������������������������������

�������������������������������������

�������������������������������������������

�����������������������������������������������

�

�
where glacial tills were present in layer 1 of the model and 

were at least 16.4 ft (5.0 m) thick (fig. 2; table 3). Recharge 

rates were 2.3 x 10-4 ft/day (7.0 x 10-5 m/day) for glacial tills 

present in model layer 1 that were greater than 3.28 ft (1.0 m) 

and less than 16.4 ft (5.0 m) thick (fig. 2; table 3).

Ground-water withdrawals from the 23 public supply 

wells and 2 irrigation wells were simulated in the model using 

the well package in MODFLOW (fig. 3). These withdrawals included withdrawals from four municipal wells for the 

city of Grand Rapids and two municipal wells for the city of 

Cohasset. In the model, water was withdrawn from either the 

middle or lower aquifers at a specified rate, where the rate was 

independent of the cell area and head (McDonald and Harbaugh, 1988). Little water is used from the upper aquifer for 

public supply or irrigation use in the study area, and therefore was not simulated in the model. Pumping rates for the 

wells were based on the highest annual pumping rate between 

1997 and 2001. Annual pumping rates for the municipal 

wells were obtained from the MNDNR Water Appropriations 

Permit Records (Minnesota Department of Natural Resources, 

2002b). A total of 383,500 ft3/day of pumping was simulated 

from the lower aquifer, whereas a total of 3,700 ft3

/day of 

pumping was simulated from the middle aquifer. Withdrawals from the Biwabik Iron Formation were not simulated in 

the model because it was assumed that the volume of water 

moving downward across the bottom of the lowest layer in 

the model was small relative to the amount of horizontal 

flow. The actual recharge could be greater than the simulated 

recharge in localized areas of the model where water is leaving the model through the bottom of the lowest model layer, 

recharging the Biwabik Iron Formation.

The Preconditioned Conjugate Gradient (PCG2) solver 

package was used with the modified incomplete Cholesky preconditioning option (relaxation parameter = 1.0) to solve the 

matrix equations produced by the model (Hill, 1990). Convergence criteria of 0.03 ft (0.01 m) of both head change and 

residual criterion were used during model calibration. Within 

the model boundaries, no-flow model cells were not allowed 

to rewet.

Model Calibration

Model calibration was done by adjusting initial estimates 

of aquifer properties and boundary conditions within a plausible range of values until simulated water levels and groundwater discharge to streams acceptably match measured water 

levels and estimated base flows, respectively. Recharge rates 

and hydraulic conductivities were varied during calibration 

because values for these variables were uncertain and initial 

model runs were most sensitive to changes in these variables.

Water levels in 351 wells were used to calibrate the 

model. The 351 wells were completed in the upper (60 wells), 

middle (240 wells), and lower (51 wells) aquifers (fig. 2). In 

general, static water levels for these wells were reported by 

well drillers following installation, and were obtained from 

water well records in the Minnesota Geological Survey’s 

County Well Index System. Measurement dates for these 

static water levels ranged from 1975 to 1999. A water-table 

contour map published by Oakes (1970) also was used as a 

reference for assessing the simulated water-table contours.

Base-flow rates into two reaches of the Mississippi River 

(eastern and western) and a reach of the Prairie River (fig. 2) 

were used in model calibration. The base-flow rate for the 

western Mississippi River reach was determined by multiplying the length of the reach by the base flow per river mile 

value determined by Payne (1995) for the subreach between 

the USGS gaging station at Lake Winnibigoshish Dam (67 

miles upstream from Grand Rapids) and the USGS gaging station at Grand Rapids (fig. 1). In a similar manner, base-flow 

rates for the eastern Mississippi River reach and the Prairie 

River reach were determined by multiplying the length of the 

reaches by the base flow per river mile value determined by 

Payne (1995) for the subreach between the USGS gaging station at Grand Rapids and the USGS gaging station near Libby, 

Minnesota (75 miles downstream from Grand Rapids) (fig.1). 

Simulation Results

The final calibration values for hydraulic conductivity 

and ground-water recharge rates are listed in table 3. The 

match between simulated and measured water levels was 

improved by (1) decreasing the horizontal hydraulic conductivity for the glacial tills from 6.6 x 10-2 to 2.0 x 10-2 ft/day, 

(2) increasing the horizontal hydraulic conductivity for the 

glaciofluvial sediments from 13 to values ranging from 98 to 

270 ft/day, (3) decreasing the vertical hydraulic conductivity for the glacial tills from 6.6 x 10-2 to 2.0 x 10-3 ft/day, (4) 

increasing the vertical hydraulic conductivity for glaciofluvial 

sediments in layers 1, 2, and 6 from 13 to values ranging 

from 16 to 27 ft/day, and (5) decreasing the vertical hydraulic 

conductivity for glaciofluvial sediments in layers 2–6 from 13 

to values ranging from 9.8 to 11 ft/day. The above changes are 

within the range of values calculated or measured in existing 

aquifer tests or reported in previous studies cited in this report.

Plots of measured and best-fit simulated water levels and 

model statistics indicate that the model tends to overestimate 

ground-water levels. Plots of measured and best-fit simulated 

water levels for the wells completed in all glaciofluvial aquifers (fig. 4a), the upper (fig. 4b), the middle (fig. 4c), and the 

lower (fig. 4d) glaciofluvial aquifers are shown with the 1:1 

linear relation between measured and simulated water levels in 

figure 4. Root mean square error, mean absolute differences, 

and mean algebraic differences for the simulated and measured water levels indicate that the simulated water levels are 

biased high compared to the measured water levels (table 4). 

The root mean square errors ranged from +12.83 ft in wells 

completed in the upper aquifer to +19.10 ft in wells completed 

in the middle aquifer. The mean absolute difference between 

simulated and measured water levels, computed as the sum of 

the absolute values of the differences divided by the number of 

Simulation of ground-water flow 13
14 Simulation of ground-water flow in glaciofluvial aquifers in the Grand Rapids area, Minnesota

����������������������������������������������������

����� ����� ����� ����� ����� ����� ����� ����� ����� ����� ����� ����� ����� ����� ����� ����� ����� �����

�������������������������������������������������������

�����

�����

�����

�����

�����

�����

�����

�����

�����

�����

�����

�����

�����

�����

�����

�����

�����

�����

�������������������������������������������������������� ����������������������������������������

����������������������������������������� ����������������������������������������

��������������������������������

��������������������������������

��������������������������������

��

���������������� ��������������

Figure 4. Measured and simulated water-level altitude in wells (a) completed in upper, middle, and lower aquifers, (b) completed in the 

upper aquifer, (c) completed in the middle aquifer, and (d) completed in the lower aquifer for the calibrated MODFLOW simulation of 

ground-water flow in glaciofluvial aquifers in the Grand Rapids area, Minnesota

Table 4. Regression statistics for measured and simulated water levels in wells and simulated and estimated base flow to reaches of the 

Mississippi and Prairie Rivers for the best-fit calibration of the MODFLOW steady-state simulation of ground-water flow in glaciofluvial 

aquifers in the Grand Rapids area, Minnesota

[ft3/day, cubic feet per day; ft, foot; values in parentheses are percentage error]

Groups of wells Well measured versus simulated water-level regression statistics (ft)

Root mean squared error Mean absolute difference Mean algebraic difference

All wells +17.55 +7.61 +12.99

Wells completed in upper aquifer +12.83 +4.43 +9.35

Wells completed in middle aquifer +19.10 +9.25 +14.44

Wells completed in lower aquifer +18.18 +6.17 +13.35

Monitored reaches Base flow to reach (ft3/day)

Simulated Estimated Simulated - estimated (ft3/day)

Eastern reach of Mississippi River 700,000 730,000 - 30,000 (-4.1)

Western reach of Mississippi River 132,000 133,000 - 1,000 (-0.8)

Reach of Prairie River 324,000 297,000 + 27,000 (+9.1)

Total reach 1,160,000 1,160,000 0 (0.0)
wells, ranged from +4.43 ft for wells completed in the upper 

aquifer to +9.25 ft in wells completed in the middle aquifer. 

These mean absolute differences are less than 10 percent of 

the differences between maximum and minimum groundwater altitudes in each aquifer. The mean algebraic difference 

between simulated and measured water levels, computed as 

the algebraic sum of the differences divided by the number of 

wells, ranged from +9.35 ft in wells completed in the upper 

aquifer to +14.44 ft in wells completed in the middle aquifer, 

indicating the positive differences were not balanced by the 

negative differences. 

Differences between simulated and measured water levels 

at the well locations ranged between –74.25 to +56.10 ft, with 

a standard deviation of 15.87 ft for all of the wells. The largest positive differences occurred in isolated wells at relatively 

low land-surface altitudes near the Mississippi River or Pokegama Lake, or at relatively high land-surface altitudes in the 

central part of study area (fig. 5). The largest negative differences occurred in the middle and lower aquifers in the vicinity 

of the Grand Rapids well field (fig. 5). One possible explanation for these large differences is extreme weather conditions 

(droughts or periods of high precipitation) during the time of 

water-level measurements, which would make the assumption 

of steady-state conditions inaccurate. Extreme weather conditions are unlikely factors because no correlation was found 

between these large differences and local precipitation records 

or the Palmer Hydrological Drought Index for Central Minnesota. Other possible explanations for the large differences 

include substantial differences between simulated land-surface 

altitude for the model grid and actual land-surface altitude at 

the wells and local areas of high or low hydraulic conductivity. Differences between simulated land-surface altitude for 

the model grid and actual land-surface altitude at some wells 

were as great as 20 ft. Positive differences at some wells in 

the central part of the study area may result from ground-water 

recharge to the Biwabik Iron Formation where the formation 

subcrops in the modeled area (figs. 2 and 5). Recharge to the 

formation from the overlying glaciofluvial sediments results 

from withdrawals from the Biwabik Iron Formation aquifer.

Percentage errors between simulated and estimated baseflow rates for the three monitored reaches all were less than 

10 percent (table 4), indicating good agreement. Percentage 

errors between simulated base-flow rates and estimated baseflow rates were –4.1 for the eastern reach of the Mississippi 

River, –0.8 for the western reach of the Mississippi River, and 

+9.1 for the reach of the Prairie River (table 4). The percentage error was 0.0 between total simulated and total estimated 

base-flow rates for the three reaches (table 4).

Simulated ground-water flow in the upper aquifer is 

discontinuous due to the discontinuous areal extent of the 

aquifer and the presence of dry cells that are mostly found in 

the higher altitudes of the aquifer (fig. 6). In general, ground 

water in this aquifer flows toward the Mississippi River and 

the chain of lakes connected with the river in the central part 

of the study area (fig. 6a). A local cone of induced drawdown 

averaging about 8 ft is present beneath the city of Grand 

Rapids well field, resulting from pumping from the lower 

aquifer. Because steady-state conditions are assumed for the 

simulation, this cone of induced drawdown is only an estimate 

of long-term average conditions. The highest altitudes in the 

upper aquifer are found in the south-central part of the study 

area east of Pokegama Lake, an area of low hydraulic conductivity, and the lowest altitudes are in the southeastern part of 

the study area near the Mississippi River (fig. 6a). 

Pumping from the lower aquifer by the city of Grand 

Rapids public-supply wells lowers the potentiometric surface 

of both the middle and lower aquifers by about 60 ft (fig. 6b 

and 6c). A radial cone of depression of greater than 1 mile in 

diameter exists in both aquifers beneath the well field. The 

highest altitudes in the middle aquifer are found in the southcentral part of the study area east of Pokegama Lake, in the 

north-central part of the study area near Shoal Lake, and in 

the eastern part of the study area near the Prairie River (fig. 

6b). The lowest altitudes are present beneath the city of Grand 

Rapids well field and in the southeastern part of the study area 

(fig. 6b). In the lower aquifer, the highest altitudes are found 

in the southeastern part of the aquifer and the lowest altitudes 

are at the city of Grand Rapids well field (fig. 6c).

A water budget computed by a steady-state model is an 

accounting of inflow to, and outflow from, the model. Inflow 

(sources of water) to the model should equal outflow (discharges). In the calibrated model, total inflow was 51.38 ft3

/s, 

whereas total outflow was 51.37 ft3/s. Of the sources of water, 

inflow from wetlands and lakes accounted for 63.0 percent 

(32.35 ft3/s), simulated recharge from the river boundaries 

accounted for 19.3 percent (9.93 ft3/s), and simulated areal 

recharge (precipitation – evaporation) accounted for 17.7 

percent of the total recharge to the model (9.10 ft3

/s). Discharge from rivers accounted for 45.7 percent (23.46 ft3

/s) and 

discharge from wetlands and lakes accounted for 45.6 percent 

(23.42 ft3/s) of the total discharges from the model. Pumping 

wells accounted for the remaining 8.7 percent (4.48 ft3/s) of 

the total simulated discharge.

Model Sensitivity to Hydraulic Conductivity, 

Areal Recharge, River-Bed Conductance, and 

General-Head-Boundary Conductance

Sensitivity analyses were conducted to determine the 

effect of changes in areal recharge, hydraulic conductivity, 

river-bed conductance, and general-head-boundary conductivity on simulated water levels and ground-water discharge 

to streams. Sensitivity analyses provides an understanding 

of the importance of these changes on simulated results, and 

identifies areas for possible model improvement. Thirty-nine 

simulations were run to assess the sensitivity of the model to 

changes in areal recharge rates, aquifer hydraulic conductivity, 

river-bed conductance, and general-head-boundary conductance. Eight of the 39 simulations were conducted varying the 

recharge rates by factors of 0.6, 0.7, 0.8, 0.9, 1.1, 1.2, 1.3, and 

1.4 times the calibrated rates (table 5). Twenty-four simulaSimulation of ground-water flow 15
16 Simulation of ground-water flow in glaciofluvial aquifers in the Grand Rapids area, Minnesota

��������������������

��������������������

��������

����

����������

����

����

����������

�������

����

����

����

�����������

����

�����������

�����

�������

�����

����������� ��

�����������������

������ ������ ������ ������

������

������

������ B

B

B

B

B

BB

B

B BB B B B B BB B B

B B

BB B

BB B

B

BBB B B

B B

B B

B B BB BB

B B B

B B B

B B B B B B B

B B B B BBB

BB B B B B BBB B B

B B B

B B B BBBB B BBB

B B B BB B

B

BBB B

BB

B B B BBB BB

B BB B

B BB

B B

BB

BB

BBB B B BBB

B

BB BB

B B B B

B B BBB B BBB B B BB B

B B B BB B

B

B

B B BB

B

B

B BB

BB BBB

B B B B

B

BB B B BBBB B B BB B BB

B B B BB BB BB BB BB B B BB B

B B BB BB B B BBB

B

B

B

B

BB

B

B

B

B

B

B

B

B �

�

�

��

��

��

��

��

��

�� ���

���

��������������������

��������������������

��������

����

����������

����

����

����������

�������

����

����

����

�����������

����

�����������

�������

�����

��

�����������������

������ ������ ������ ������

������

������

������ B

B

B B

B

B

B BBB

B

B

B

B

BBB

B

B

B

B

B

BB

BBBB

B

BB

B

B

BB

B

B

B

B

B

B

BB

B

B B

B

B

BB

B

B

B B BB

BB

BB

B

B

BB

B

BB

B B

B B

BB

B

B

B B BB

B

BB

B BB

BB

BB

BBB

BB

B

B

B

B B

B

B

B

B

B

�

�

��

��

��

��

��

��

��

�� �

���

���

��

�

�

�

��

�

�� ��

�

Figure 5. Difference between simulated and measured water-level contours for ground-water flow for the (a) Model Layer two—
������������

���������

�������������������

����

�����������

B

����������

����������

����������

����������

����������

����������

����������

�������

������

��������

��������

��������

��������

��������

��������

�

���������������������������

�����������������������

������������������

� � � � � �������

� � � � � ������������

�������������������

�������������������

��������

����

����������

����

����

����������

�������

����

����

����

�����������

����

�����������

�������

�����

����������� ��

���

��

����

��

��

�

�

�

��

���

�����������������

������ ������ ������ ������

������

������

������

B

B

B

B B

B

B BB

B B B

B BB

BBB

B

B BB B

BBBBB

B B

BBB

B

B B

BB BBBBB

BBB

B

B

B ���

��� ���

���

��� ���

��� ���

���

Upper aquifer, (b) Model Layer four—Middle aquifer, and (c) Model Layer six—Lower aquifer, in the Grand Rapids area, Minnesota

Simulation of ground-water flow 17
18 Simulation of ground-water flow in glaciofluvial aquifers in the Grand Rapids area, Minnesota

��������������������

��������������������

��������

����

����������

����

����

����������

�������

����

����

����

�����������

����

�����������

�����

�������

�����

����������� ��

�����������������

������ ������ ������ ������

������

������

������

����

����

���� ����

����

����

����

����

����

����

����

����

����

����

����

����

����

����

����

����

����

��������

����

����

����

����

����

����

�

�

�

��������������������

��������������������

��������

����

����������

����

����

����������

�������

����

����

����

����������� ����

�����������

������� �����

����������� ��

�����������������

������ ������ ������ ������

������

������

������

����

����

����

����

����

����

����

����

����

����

����

����

����

���� ����

����

����

����

����

����

����

����

���� ����

����

����

����

Figure 6. Simulated potentiometric surfaces for ground-water flow and ground-water flow direction for the (a) Model Layer two—
������������

���������

�������������������

����������������������������������������������

����������������������������������������

���������������������������������������������

����������������������������������������

�������������������������������������

�������������������������������������������

�����������������������������������������������

������������������������

��������������������

�����������

��������������

��������������

��������������

��������������

��������������

��������������

��������������

��������������

��������������

��������������

��������������

����������������������

�����������������������

������������������

�

�

�������������������

�������������������

��������

����

����������

����

����

����������

�������

����

����

����

�����������

����

�����������

�������

�����

����������� ��

�����������������

������ ������ ������ ������

������

������

������

����

����

����

����

����

���� ����

����

����

����

����

����

�

�

�

�

� � � � � �������

� � � � � ������������

Upper aquifer, (b) Model Layer four—Middle aquifer, and (c) Model Layer six—Lower aquifer, in the Grand Rapids area, Minnesota

Simulation of ground-water flow 19
20 Simulation of ground-water flow in glaciofluvial aquifers in the Grand Rapids area, Minnesota

 Table 5. Relation among model variables used in sensitivity analysis for the calibrated MODFLOW simulation of ground-water flow in 

the upper, middle, and lower aquifers in the Grand Rapids area, Minnesota

[---, no value]

Groups of wells Multiplication 

factor

Root mean square error in ground-water levels (ft)

Areal recharge Horizontal 

hydraulic 

conductivity

Vertical 

hydraulic 

conductivity

River-bed 

conductance

General-head 

boundary 

conductance

Wells completed in upper 

aquifer

0.6 12.70 14.01 12.99 12.83 12.21

0.7 12.73 13.65 12.89 12.83 12.30

0.8 12.76 13.35 12.83 12.83 12.50

0.9 12.80 13.06 12.83 12.83 12.66

1.0 12.83 12.83 12.83 12.83 12.83

1.1 12.86 12.66 12.80 12.83 12.96

1.2 12.89 12.50 12.80 12.83 13.12

1.3 12.89 12.57 12.80 12.83 13.26

1.4 12.93 --- 12.80 12.83 13.39

Wells completed in middle 

aquifer

0.6 19.03 19.00 18.83 19.10 18.18

0.7 18.96 18.93 18.73 19.10 18.28

0.8 19.00 18.83 18.67 19.10 18.93

0.9 19.03 18.64 18.64 19.10 19.03

1.0 19.10 19.10 19.10 19.10 19.10

1.1 19.13 18.96 19.03 19.10 19.16

1.2 19.16 18.87 18.96 19.10 19.19

1.3 19.16 18.83 18.93 19.10 19.26

1.4 19.03 --- 18.90 19.10 19.33

Wells completed in lower 

aquifer

0.6 18.05 17.95 18.18 18.21 17.88

0.7 18.05 17.91 18.57 18.21 17.98

0.8 18.08 17.68 18.37 18.18 18.08

0.9 18.14 17.88 18.28 18.18 18.14

1.0 18.18 18.18 18.18 18.18 18.18

1.1 18.24 18.01 18.11 18.18 18.24

1.2 18.28 17.88 18.08 18.18 18.31

1.3 18.28 17.75 18.01 18.18 18.31

1.4 18.24 --- 17.98 18.18 18.41

Monitored reaches Multiplication 

factor

Percentage error in base flow

Areal Recharge Horizontal 

hydraulic 

conductivity

Vertical 

hydraulic 

conductivity

River-bed 

conductance

General-head 

boundary 

conductance

Eastern reach of Mississippi 

River

0.6 -9.1 -27.4 -9.6 -5.0 -15.0

0.7 -7.0 -20.1 -8.4 -4.7 -12.0

0.8 -6.1 -14.8 -6.9 -4.5 -9.5

0.9 -5.2 -9.7 -5.6 -4.4 -6.4

1.0 -4.1 -4.1 -4.1 -4.1 -4.1

1.1 -3.4 1.0 -3.2 -4.3 -2.4

1.2 -2.5 5.3 -2.1 -4.2 -0.7

1.3 -1.0 8.0 -1.1 -4.0 -1.1

1.4 -0.9 --- -0.1 -4.0 -2.7

Western reach of Mississippi 

River

0.6 -5.3 -24.1 -5.0 -0.5 -12.8

0.7 -3.9 -17.3 -4.1 -0.4 -9.1

0.8 -2.4 -11.1 -2.2 -0.2 -5.8

0.9 -1.6 -6.0 -1.1 -0.1 -2.8

1.0 -0.8 -0.8 -0.8 -0.1 -0.8

1.1 0.7 4.9 0.9 -0.1 1.7

1.2 2.1 9.5 1.9 -0.1 4.8

1.3 3.2 13.7 2.8 0.0 7.0

1.4 4.2 --- 3.7 -0.1 8.3

Reach of Prairie River 0.6 4.8 -5.6 -7.1 8.6 -1.6

0.7 5.8 1.2 -3.1 8.8 1.7

0.8 5.5 -1.8 0.1 8.9 4.7

0.9 8.0 3.8 6.0 9.0 7.0

1.0 9.1 9.1 9.1 9.1 9.1

1.1 10.0 13.7 12.0 9.0 10.8

1.2 10.9 15.0 14.8 9.2 12.2

1.3 11.7 16.6 20.8 9.2 13.6

1.4 10.0 --- 23.5 9.2 14.8
tions were run using various vertical hydraulic conductivity, 

river-bed conductance, and general-head-boundary values, 

varying the values by the same multiplication factors as in the 

areal recharge simulations. Horizontal hydraulic conductivity 

values were varied by factors of 0.6, 0.7, 0.8, 0.9, 1.1, 1.2, and 

1.3 times the calibrated values in seven simulations. Sensitivity simulations were compared to the calibrated model simulations by observing changes in the root mean square error in 

water levels in the wells and percentage error in the base-flow 

rates to the Mississippi and Prairie River reaches.

Simulated water levels were most sensitive to changes in 

general-head-boundary conductance (fig. 7 and table 5), indicating that it is the predominant variable controlling steadystate water levels. Simulated base-flow for the model was 

most sensitive to changes in horizontal hydraulic conductivity 

(fig. 7, table 5), indicating that it is the predominant variable controlling steady-state flow conditions. Compared to 

general-head-boundary conductance and horizontal hydraulic 

conductivity, the model was less sensitive to changes in vertical hydraulic conductivity and areal recharge values (table 5). 

The model was least sensitive to the changes in river-bed 

conductance (table 5 and fig. 7). This low sensitivity may be 

due in part to the smaller values for river-bed conductance 

relative to the horizontal hydraulic conductivity values used in 

the model and the smaller areal extent of the river beds relative 

to the extent of the aquifers, general-head boundaries, and 

areal recharge.

The model sensitivity varied among the three aquifers and 

among the different reaches. Simulated ground-water levels 

in the upper aquifer were most sensitive to changes in horizontal hydraulic conductivity (table 5), whereas ground-water 

levels in the middle aquifer were most sensitive to changes in 

general-head boundary conductance. Simulated ground-water 

levels in the lower aquifer were equally sensitive to changes in 

horizontal hydraulic conductivity, vertical hydraulic conductivity, and general-head boundary conductance (table 5). Base 

flow to the simulated reaches for the Mississippi River were 

most sensitive to changes in horizontal hydraulic conductivity, 

whereas base flow to the Prairie River reach was most sensitive to changes in vertical hydraulic conductivity of the aquifer 

(table 5). 

�����������������������������������������

�����������������������������������������

��������������������������������

�����������������������

��

����

����

����

����

����

��

����

��

��

�

�

��

���

���

���

���

��� ��� ��� ��� ��� ��� ��� ��� ��� ��� ��� ��� ��� ��� ��� ��� ��� ��� ��� ��� ��� ���

��������������������� ���������������������

���������������������������������������������������������

�������������������

��������������������������������������������������������������

��������������

���������������������������������

�������������������������������

���������������������

���������������������������������

�����������

Figure 7. Model sensitivity to selected input variables based on (a) root mean square error between simulated and measured groundwater levels, and (b) percentage error between simulated and estimated base flow for the calibrated MODFLOW simulation of groundwater flow in glaciofluvial aquifers in the Grand Rapids area, Minnesota

Simulation of ground-water flow 21
22 Simulation of ground-water flow in glaciofluvial aquifers in the Grand Rapids area, Minnesota

Model Limitations and Accuracy

The numerical flow model is a simplification of a 

complex glaciated-terrain, ground-water flow system. The 

accuracy of the simulation is limited to the accuracy, amount, 

and distribution of the data used to describe the hydrologic 

variables of the flow system and the accuracy of the assumption of steady-state. These variables include the hydrologic 

properties of the aquifers and confining units, areal recharge 

rates, and hydrologic boundary conditions. Variables determined for the model during calibration may not be unique. 

Different combinations of variables for the model could 

produce similar results from the model. The model grid was 

designed to simulate hydrologic conditions at the scale of 

the study area. For example, the model can be modified to 

include storage terms and used with particle-tracking software 

(Pollock, 1994) to determine approximate areal extents of the 

zones of contribution for discharging wells in the study area 

on a regional scale. However, the model grid could be refined 

to address hydrologic issues on a smaller, more accurate scale, 

such as the estimation of the areal extent of the zone of contribution for individual public supply wells.

The accuracy of a model depends upon the accuracy 

of the measured data used to calibrate a model. The model 

was calibrated using water-level data recorded, for the most 

part, within a week following well installation; therefore, it is 

unlikely that static conditions were present for all of the wells 

following installation. These wells were installed between 

1975 and 1999, when hydrologic conditions varied between 

periods of droughts to periods of above-normal precipitation. 

Thus, these water levels do not necessarily represent steadystate conditions. Also, the base-flow rates used for calibration 

were estimated from values obtained by Payne (1995) over 

large reaches of the Mississippi River. Any local variations 

in the base-flow values were not taken into account in the 

estimated values. Therefore, the model may not accurately 

represent small-scale flow conditions where local variations 

from more regional patterns may exist.

The calibrated, steady-state model is a tool for waterresources management based on the assumption that future 

hydrologic conditions will be similar to historical, or assumed 

steady-state, conditions. Furthermore, it would be assumed 

that variations in annual recharge and discharge for a hypothetical simulation using the calibrated model are similar to 

hydrologic conditions of the calibrated model. Accuracy of 

the hypothetical simulation becomes more uncertain if the 

variation in annual recharge or discharge exceeds the range 

used in the calibrated model. Because the numerical model 

was calibrated under the assumption of steady-state flow 

conditions, the model will most accurately reflect the effects 

of long-term, annual or multiple-year stresses.

SUMMARY

A calibrated steady-state, finite-difference, ground-water 

flow model was constructed to simulate a system of three 

glaciofluvial aquifers, defined as the upper, middle, and lower 

aquifers, in an area of about 114 mi2 surrounding the city of 

Grand Rapids in north-central Minnesota. This model was 

constructed as part of a 5-year cooperative study between the 

U.S. Geological Survey (USGS) and the Minnesota Department of Health (MDH) to simulate ground-water-flow conditions in aquifers. The calibrated model will be used by MDH 

and communities in the Grand Rapids area in the development 

of wellhead protection plans for their water supplies. 

The USGS MODular ground-water-FLOW model, MODFLOW, was used to simulate ground-water flow conditions 

in the Grand Rapids area. Available geologic and hydrologic 

data collected from 1975 through 2001 were used for model 

construction and calibration. Data from scientific literature 

and other sources also were used. A series of GIS coverages 

were used to identify and represent hydrologic features and 

processes in the model. Data sets for perennial wetlands, 

lakes, and rivers were obtained from the National Wetlands 

Inventory data base. A “true-layer” approach was used to 

define layering represented by the model, explicitly defining altitudes and aquifer hydraulic properties of cells in each 

model layer based on the three-dimensional geologic model 

representation. Perennial rivers, streams, lakes, and wetlands 

were specified as head-dependent boundaries in five of the 

six layers in the model. Lakes and wetlands were simulated 

as general-head boundaries, and streams and were simulated 

using the river package of MODFLOW. A specified-flux 

boundary was used to represent areal recharge to the upper 

layer (layer 1) of the model using the recharge package in 

MODFLOW. Ground-water withdrawals from the 23 public supply wells and 2 irrigation wells were simulated in the 

model using the well package in MODFLOW. 

The model was calibrated through comparison of simulated ground-water levels to measured static water levels in 

351 wells, and comparison of simulated base-flow rates to 

estimated base-flow rates for two reaches of the Mississippi 

River (eastern and western) and a reach of the Prairie River. 

Measurement dates for the static water levels in the wells 

ranged from 1975 to 1999. Model statistics indicate that the 

model tends to overestimate ground-water levels. The root 

mean square errors ranged from +12.83 ft in wells completed 

in the upper aquifer to +19.10 ft in wells completed in the 

middle aquifer. Mean absolute differences between simulated 

and measured water levels ranged from +4.43 ft for wells completed in the upper aquifer to +9.25 ft for wells completed in 

the middle aquifer. Mean algebraic differences ranged from 

+9.35 ft for wells completed in the upper aquifer to +14.44 ft 

for wells completed in the middle aquifer, with the positive 

differences indicating that positive water levels were not balanced by the negative differences. Percentage errors between 

simulated and estimated base-flow rates for the three moni-
tored reaches were less that 10 percent, indicating good agreement. Simulated ground-water levels were most sensitive to 

changes in general-head boundary conductance, indicating 

that this characteristic is the predominant model input variable controlling steady-state water-level conditions. Simulated 

ground-water flow to stream reaches was most sensitive to 

changes in horizontal hydraulic conductivity, indicating that 

this characteristic is the predominant model input variable 

controlling steady-state flow conditions.

The accuracy of the model is limited to the accuracy, 

amount, and distribution of the data used to describe the 

hydrologic variables of the flow system and the accuracy of 

the assumption of steady-state. The model grid was designed 

to simulate hydrologic conditions at the scale of the study 

area. The model may not accurately represent small-scale 

flow conditions where local variations from more regional patterns may exist.

ACKNOWLEDGMENTS

The author thanks the Minnesota Department of Health, 

in particular Mr. James Walsh and Mr. Bruce Olsen, for providing assistance needed to complete this study. The author 

thanks the Minnesota Department of Natural Resources, in 

particular Mr. Sean Hunt and Mr. Kevin Zytkovicz, for providing annual pumping records for 1997-2001 for public supply 

wells and bathymetric data for lakes in the study area. The 

author thanks Minnesota Power, in particular Lowell Neudahl, 

for providing bathymetric data for Blackwater Lake.

REFERENCES

Harbaugh, A.W., and McDonald, M.G., 1996, User’s documentation for MODFLOW-96, an update to the U.S. Geological Survey Modular Finite-Difference Ground-Water 

Flow Model: U.S. Geological Survey Open-File Report 

96-485, 63 p.

Hill, M.C., 1990, Preconditioned conjugate-gradient 2 

(PCG2)—A computer program for solving ground-water 

flow equations: U.S. Geological Survey Water-Resources 

Investigations Report 90-4048, 43 p.

Jirsa, M.A., Setterholm, D.R., Bloomgren, B.A., Bauer, E.J., 

and Lively, R.S., 2002, Bedrock topographic and depth to 

bedrock maps of the western half of the Mesabi Iron Range, 

Northern Minnesota: Minnesota Geological Survey, Miscellaneous Map Series, Map M-126.

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, 30 p.

Kishel, H.F., and Gerla, P.J., 2002, Effects of small-scale 

hydraulic variation on ground-water discharge to Shingobee Lake, Minnesota: The University of North Dakota, 

Department of Geology web page, accessed April 23, 2002, 

at URL http://www.geology.und.nodak.edu/gerla/home/

cur_res.html

McDonald, M.G., and Harbaugh, A.W., 1988, A modular 

three-dimensional finite-difference ground-water flow 

model: U.S. Geological Survey Techniques of WaterResources Investigations, book 6, chap A-1, p. 2-1 – 2-35.

Minnesota Department of Health, 2003, Community Public 

Water Supply General Information: accessed March 18, 

2003, at URL http://www.health.state.mn.us/divs/eh/water/

index.htm

Minnesota Department of Natural Resources, 2002a, LandSatbased land use-land cover (Raster), 1995 (Minnesota): Minnesota Department of Natural Resources online metadata, 

accessed February 11, 2002, at URL http://deli.dnr.state.

mn.us/metadata/full/lusatra3.html

Minnesota Department of Natural Resources, 2002b, Minnesota Water Use Information (Raster), 1996 (Minnesota): 

Minnesota Department of Natural Resources online metadata, accessed February 14, 2002, at URL http://www.dnr.

state.mn.us/waters/programs/water_mgt_section/appropriations/wateruse.html

Minnesota Department of Natural Resources, 2002c, Lake 

finder, 1996: accessed March 25, 2002, at URL http://www.

dnr.state.mn.us/lakefind/index.html

Minnesota State Climatologist, 2001, Normal precipitation and 

temperature data: accessed July 11, 2001, at URL http://climate.umn.edu/doc/historical/normals.htm

Morey, G.B., 1972, Mesabi Range, in Sims, P.K., and Morey, 

G.B., eds., Geology of Minnesota—A Centennial Volume: 

Minnesota Geological Survey, p. 204-217.

Oakes, E.L., 1970, Geology and ground-water resources of the 

Grand Rapids area, north-central Minnesota: U.S. Geological Survey Hydrologic Investigations Atlas, HA-322, scale 

1:48,000.

Payne, G.A., 1995, Ground-water base flow to the Upper 

Mississippi River upstream of the Minneapolis-St. Paul 

Area, Minnesota during July 1988: U.S. Geological Survey 

Open-File Report, 94-478, 28 p.

Pollock, D.W., 1994, User’s guide for MODPATH/MODPATH-PLOT, Version 3—A particle tracking post-processing package for MODFLOW, the U. S. Geological Survey 

finite-difference ground-water flow model: U.S. Geological 

Survey Open-File Report 94-464, 249 p.

Winter, T.C., 1971, Sequence of glaciation in the MesabiVermilion Iron Range area, northeastern Minnesota: U.S. 

Geological Survey Professional Paper 750-C, p. C82-C88.

Winter, T.C., 1973a, Hydrogeology of glacial drift, Mesabi 

Iron Range, northeastern Minnesota: U.S. Geological Survey Water-Supply Paper 2029-A, 23 p.

Winter, T.C., 1973b, Petrography and stratigraphy of glacial 

drift, Mesabi-Vermilion Iron Range area, northeastern Minnesota: U.S. Geological Survey Bulletin 1331-C, 41 p.

U.S. Census Bureau, 2002, Population data: accessed February 

6, 2002, at URL http://factfinder.census.gov/servlet/BasicFactsServlet

U.S. Department of Defense, 1998, Groundwater modeling 

system, user’s manual: Brigham Young University, 620 p.

U.S. Fish & Wildlife Service, 1994, National wetlands inventory digital data, October 1991 to 1994: accessed February 

11, 2002, at URL http://www.nwi.fws.gov

References 23