U.S. Department of the Interior

U.S. Geological Survey

Scientific Investigations Report 2010–5117

National Water Availability and Use Pilot Program

Implementation of Local Grid Refinement (LGR) for the 

Lake Michigan Basin Regional Groundwater-Flow Model

Implementation of Local Grid Refinement 

(LGR) for the Lake Michigan Basin Regional 

Groundwater-Flow Model

By C.J. Hoard

National Water Availability and Use Pilot Program 

Scientific Investigations Report 2010–5117

U.S. Department of the Interior

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

KEN SALAZAR, Secretary

U.S. Geological Survey

Marcia K. McNutt, Director

U.S. Geological Survey, Reston, Virginia 2010

For product and ordering information:

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

Telephone: 1-888-ASK-USGS

For more information on the USGS—the Federal source for science about the Earth,

its natural and living resources, natural hazards, and the environment:

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

Telephone: 1-888-ASK-USGS

Suggested citation:

Hoard, C.J., 2010, Implementation of local grid refinement (LGR) for the Lake Michigan Basin regional 

groundwater-flow model: U.S. Geological Survey Scientific Investigations Report 2010–5117, 25 p.

Any use of trade, product, or firm names 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 material contained within this report.

ISBN 978-1-4113-2941-0
iii

Contents

Abbreviations and Acronyms.......................................................................................................................v

Acknowledgements.......................................................................................................................................v

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

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

Purpose and Scope ..............................................................................................................................2

Regional Groundwater-Flow Model...................................................................................................3

Model-Grid Refinement.................................................................................................................................3

Refined-Model Development..............................................................................................................3

Construction of Intermediate and Local Model Grids ...........................................................7

Boundary Conditions of the Refined Grids ..............................................................................7

Recharge................................................................................................................................................7

Surface-Water Network ......................................................................................................................9

Refined-Model Scenarios...........................................................................................................................13

Streamflow Depletion ........................................................................................................................13

Climate-Change Scenario .................................................................................................................16

Assumptions and Limitations of the Models ...........................................................................................19

Summary........................................................................................................................................................19

References Cited..........................................................................................................................................21

Appendix 1. Surface-Water Network....................................................................................................23

Figures

1. Location and features of the regional, intermediate, and local models..............................4

2. Geologic cross sections for the regional and local models..................................................5

3. Example of shared-node distribution between the intermediate and local models .........8

4. Seasonal distribution of recharge for the adjusted and non-adjusted recharge 

values............................................................................................................................................10

5. Surface-water network for the area of the local model shown in comparison to 

how it was discretized at the levels of grid refinement used in the regional 

and intermediate models...........................................................................................................11

6. Stream network for the intermediate- and local-model areas.................................................12

7. Location of test watershed, test stream, stream gage, and pumped well for 

refined-model scenarios ...........................................................................................................14

8. Hydrographs of test stream for streamflow-depletion scenarios......................................15

9. Distribution of the source of water to the pumping well for the stream-depletion 

scenarios......................................................................................................................................17

10. Recharge and streamflow distribution shown for three time periods for the test 

stream in the climate-change scenario..................................................................................18
iv

Tables

1. Summary of physical properties used for the intermediate model......................................6

2. Summary of physical properties used for the local model....................................................6

3. Summary of the overall water budget in the test watershed with a well pumping 

in layer 1 .......................................................................................................................................16

Appendix figures

1–1. Stream-segment construction across the intermediate and local-model boundary.............25

Conversion Factors

Multiply By To obtain

Length

foot (ft) 0.3048 meter (m)

meter (m) 3.281 foot (ft)

mile (mi) 1.609 kilometer (km)

Area

square foot (ft2) 0.09290 square meter (m2)

square mile (mi2) 259.0 hectare (ha)

square mile (mi2) 2.590 square kilometer (km2) 

Volume

gallon (gal) 3.785 liter (L) 

gallon (gal) 0.003785 cubic meter (m3) 

million gallons (Mgal) 3,785 cubic meter (m3)

cubic foot (ft3) 0.02832 cubic meter (m3) 

Flow rate

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

cubic foot per day (ft3/d) 0.02832 cubic meter per day (m3/d)

gallon per minute (gal/min) 0.06309 liter per second (L/s)

million gallons per day (Mgal/d) 0.04381 cubic meter per second (m3/s)

Hydraulic conductivity

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

Hydraulic gradient

foot per mile (ft/mi) 0.1894 meter per kilometer (m/km)

Vertical coordinate information is referenced to the North American Vertical Datum of 1988 

(NAVD 88).

Horizontal coordinate information is referenced to the North American Datum of 1983 (NAD 83).

Elevation, as used in this report, refers to distance above the vertical datum.

Concentrations of chemical constituents in water are given in milligrams per liter (mg/L).
v

Abbreviations and Acronyms

AFINCH Analysis of Flows in Networks of Channels

AOGCM Atmosphere Ocean Coupled General Circulation Model

CMIP3 Coupled Model Intercomparison Project 3

DEM Digital Elevation Model

GHB General Head Boundary Package

GIS Geographic Information System

GLBP Great Lakes Basin Pilot

HUC Hydrologic Unit Code

LGR Local Grid Refinement

LPF Layer Property Flow Package

NHD National Hydrography Dataset

PCMDI Program for Climate Model Diagnosis and Intercomparison

SFR1 Streamflow Routing Package

SWB Soil Water Balance

TMR Telescopic Mesh Refinement

USGS United States Geological Survey

WCRP World Climate Research Programs

WGCM Working Group on Coupled Modelling

Acknowledgments

The author wishes to acknowledge the modeling groups, the Program for Climate Model 

Diagnosis and Intercomparison (PCMDI) and the World Climate Research Programs (WCRP) 

Working Group on Coupled Modelling (WGCM) for their roles in making available the WCRP 

Coupled Model Intercomparison Project 3 (CMIP3) multi-model dataset. Support of this dataset is provided by the Office of Science, U.S. Department of Energy. Additional simulations 

provided courtesy of Katherine Hayhoe, Edwin Maurer, and Jeff VanDorn of ATMOS Research 

and Consulting.

The author also wishes to thank Daniel Feinstein, Mike Fienen, Randy Hunt, Steffen 

Mehl, Howard Reeves, and Steve Westenbroek from the USGS for the assistance in planning 

and developing the various simulations described in this report. Jesse Dickinson and Connor 

Haugh, also from the USGS, provided thorough technical reviews of this report. A complete 

editorial review of this report was provided by Bonnie Fink (USGS ).

Implementation of Local Grid Refinement (LGR) for the 

Lake Michigan Basin Regional Groundwater-Flow Model

By C.J. Hoard

Abstract

The U.S. Geological Survey is evaluating water availability and use within the Great Lakes Basin. This is a pilot 

effort to develop new techniques and methods to aid in the 

assessment of water availability. As part of the pilot program, 

a regional groundwater-flow model for the Lake Michigan 

Basin was developed using SEAWAT-2000. The regional 

model was used as a framework for assessing local-scale water 

availability through grid-refinement techniques. Two gridrefinement techniques, telescopic mesh refinement and local 

grid refinement, were used to illustrate the capability of the 

regional model to evaluate local-scale problems. An intermediate model was developed in central Michigan spanning an area 

of 454 square miles (mi2) using telescopic mesh refinement. 

Within the intermediate model, a smaller local model covering 

an area of 21.7 mi2 was developed and simulated using local 

grid refinement. Recharge was distributed in space and time 

using a daily output from a modified Thornthwaite-Mather 

soil-water-balance method. The soil-water-balance method 

derived recharge estimates from temperature and precipitation data output from an atmosphere-ocean coupled generalcirculation model. The particular atmosphere-ocean coupled 

general-circulation model used, simulated climate change 

caused by high global greenhouse-gas emissions to the atmosphere. The surface-water network simulated in the regional 

model was refined and simulated using a streamflow-routing 

package for MODFLOW. 

The refined models were used to demonstrate streamflow 

depletion and potential climate change using five scenarios. 

The streamflow-depletion scenarios include (1) natural conditions (no pumping), (2) a pumping well near a stream; the well 

is screened in surficial glacial deposits, (3) a pumping well 

near a stream; the well is screened in deeper glacial deposits, 

and (4) a pumping well near a stream; the well is open to a 

deep bedrock aquifer. Results indicated that a range of 59 to 

50 percent of the water pumped originated from the stream 

for the shallow glacial and deep bedrock pumping scenarios, 

respectively. The difference in streamflow reduction between 

the shallow and deep pumping scenarios was compensated for 

in the deep well by deriving more water from regional sources. 

The climate-change scenario only simulated natural conditions 

from 1991–2044, so there was no pumping stress simulated. 

Streamflows were calculated for the simulated period and 

indicated that recharge over the period generally increased 

from the start of the simulation until approximately 2017, and 

decreased from then to the end of the simulation. Streamflow 

was highly correlated with recharge so that the lowest streamflows occurred in the later stress periods of the model when 

recharge was lowest.

Introduction

In 2005, the U.S. Congress requested the U.S. Geological Survey (USGS) to initiate an assessment of the availability and use of water resources and the factors that affect 

water resources throughout the United States. The goal of the 

National Assessment of Water Availability and Use Program 

is to improve methods for predicting water availability for 

future economic and environmental uses (Grannemann and 

Reeves, 2005). This program began as a pilot study within the 

Great Lakes Basin to develop the tools and methods necessary to study water availability and to disseminate the results. 

The techniques developed can be used in future studies as the 

National Assessment of Water Availability and Use Program 

expands into other areas of the Nation.

As demands on water resources increase, it becomes 

more important to quantify the water resources available to 

effectively manage their use. Groundwater and surface water 

are best treated as a single resource—an important consideration as water-resource managers assess environmental flows, 

or those flows required to sustain aquatic ecosystems while 

maintaining anthropogenic uses of the resource (The Nature 

Conservancy, 2008). The environmental flows in streams 

reflect the need for a better understanding of stream-aquifer 

interaction, especially when groundwater withdrawals may 

affect streamflow- and groundwater-dependent ecosystems. 

As part of the Great Lakes Basin Pilot (GLBP) study, 

a regional groundwater-flow model was constructed for the 

Lake Michigan Basin (Feinstein and others, in press). The 

objectives of the regional model were to provide (1) a forecasting tool to assess regional effects of water use and climate 

in the western portion of the Great Lakes Basin; (2) a platform 
2 Implementation of Local Grid Refinement (LGR) for the Lake Michigan Basin Regional Groundwater-Flow Model

model and a less accurate solution by the parent model. Head 

and flux conditions taken from the parent and used initially 

along the child boundaries are allowed to change, and as 

changes occur they are propagated back through to the parent 

model. The disadvantage of TMR coupling is that changes to 

the child model are not propagated out to the parent model 

without user intervention. The primary disadvantage of the 

LGR approach is the longer runtimes and relatively more 

complex input requirement. Both TMR and LGR approaches 

are able to use the streamflow-routing (SFR1) package (Prudic 

and others, 2004) of MODFLOW to improve the evaluation of 

stream-aquifer interactions within the model.

As part of the effort to demonstrate the tools and methods available for assessing water availability, grid-refinement 

techniques were applied to the Lake Michigan Basin regional 

model. These grid-refinement techniques were used at multiple 

scales to achieve a model cell size capable of accurately 

modeling stream-aquifer interactions. Stream-aquifer interactions are of particular importance in evaluating how groundwater withdrawals potentially affect environmental flows. 

To that end, hypothetical model scenarios of refined portions 

of the regional model were developed to evaluate the effects 

of groundwater withdrawals on streams. It should be noted 

that streamflows required to support specific aquatic habitat 

were not evaluated as part of this investigation. In addition, 

hypothetical model scenarios examining the potential effect 

of climate change on groundwater recharge and the resultant 

effect on the water resources of the area of interest also were 

developed.

Purpose and Scope

The purpose of this report is to describe the use of gridrefinement techniques for MODFLOW to examine local-scale 

hydrology within a regional model. This modeling effort 

builds on the regional model developed for the Lake Michigan Basin as part of the GLBP project of the USGS national 

assessment of water availability and use. Grid refinement was 

used to demonstrate (1) local-scale groundwater and surfacewater interactions—specifically the effects of groundwater 

withdrawals on streamflow, and (2) potential climate change 

effects on water resources at a local scale. Two stages of grid 

refinement were performed on the regional model: one on an 

area of 453.8 mi2 and another on a smaller subset of the initial 

refined area of 22.4 mi2. The National Hydrography Dataset 

(NHD) NHDPlus database (Horizon Systems Corporation, 

2006) was used to develop the surface-water networks used 

within the refined model areas. Climate projections based on 

Hayhoe and others (2008) and Maurer and others (2002) were 

used to develop recharge estimates for evaluation of climate 

effects on water resources for a 53-year period (1991–2044).

for developing embedded, higher-resolution models that can 

be used to address water-management issues at local scales; 

(3) a means of documenting and archiving information from 

a wide variety of sources on the hydrogeology and water use 

in the region; and (4) a basis for developing indicators of sustainability of water resources (Feinstein and others, in press). 

The results of this modeling effort have yielded (1) improved 

estimates of the various components of the water budget for 

the region; (2) synthesis and improved estimates of the various 

hydraulic properties of the geologic units in the region; (3) 

a better understanding of groundwater flow throughout the 

region; and (4) a tool for water managers in the region to aid 

in evaluating the sustainability of water resources. 

Although the Lake Michigan Basin groundwater-flow 

model has provided a tool to assess water availability for the 

region, the large grid-cell structure needed in the regional 

model (5,000 by 5,000 ft) likely is too coarse to address many 

smaller-scale problems. Regional models are designed to provide answers to questions on a large scale, such as identifying 

and quantifying the various components of the regional water 

budget. Often the issues related to groundwater and surfacewater interaction occur at a smaller scale than the size of a 

model cell in the regional model. 

Situations that benefit from improved simulation accuracy obtained from fine grid resolution include, but are not 

limited to, the assessment of drawdown near pumping centers 

and the assessment of surface-water and groundwater interaction near streams. Regional models can be refined by variablegrid spacing to achieve small grid sizes in the area of interest, or by global refinement using a small cell size over the 

whole region, or by locally refining the model through various 

grid-refinement techniques. Model refinement of the regional 

model by either of the first two options generates a vast 

number of cells in the overall model. Although this improves 

model accuracy, it also leads to lengthy computation time. 

Using grid-refinement techniques limits the cell refinement to 

the area of interest, which allows for improved model accuracy only in the area of interest.

Two grid-refinement techniques developed for MODFLOW (Harbaugh, 2005)—telescopic mesh refinement (TMR) 

(Leake and Claar, 1999) and local grid refinement (LGR) 

(Mehl and Hill, 2006)—allow for grid refinement using a 

fine grid “child” model within a coarse grid “parent” model. 

However, the difference between the two methods is in how 

the child model is linked to the parent model. With TMR, the 

coupling between the parent and child model is one way, such 

that each model is simulated once, and the boundary conditions for the refined model are established using head and flux 

boundary conditions determined from the parent model; the 

child model is then run independently. The LGR technique 

couples the parent and child model in such a way that the 

models are run sequentially and iterate to a final solution that 

is a compromise between a more accurate solution by the child 
Model-Grid Refinement 3

Regional Groundwater-Flow Model

The Lake Michigan Basin regional groundwater-flow 

model encompasses nearly 180,963 mi2 of the land contributing drainage to Lake Michigan and the surrounding area, 

which includes parts of Wisconsin, Illinois, Indiana, Ohio, 

Michigan, and Ontario, Canada (fig. 1). Laterally, the model 

grid is divided into a series of 391 rows and 261 columns of 

variably spaced cells. Model cells range in area from 0.9 mi2

(5,000 ft by 5,000 ft) in the center of the model to 288 mi2

(68,930 ft by 116,490 ft) at the external boundaries of the 

model. The variably spaced grid cells were established so that 

the refined portion in the center of the model, or near-field, 

(fig. 1) overlies the major pumping centers for greater modelsolution accuracy. 

The model simulates groundwater flow in geologic units 

within the region, ranging in age from the glacial deposits of 

Quaternary age at the land surface to the bedrock of Precambrian age at depth. The model is vertically discretized into 20 

variable-thickness layers. The top three layers of the model 

represent the unconsolidated glacial and alluvial deposits 

within the basin. The remaining 17 layers simulate the stratified bedrock units of the Wisconsin, Kankakee, Findlay, and 

Algonquin Arches and the Michigan Basin (Feinstein and others, in press). In total, the maximum vertical thickness of the 

model reaches approximately 15,000 ft. Initial aquifer properties used in the model were gathered from previous studies 

and a variety of databases in the region (Feinstein and others, 

in press; Lampe, 2009; Arihood, 2009). Groundwater-flow 

conditions were simulated over the time period 1864–2005. 

A version of the regional model that simulated all geologic 

units as initially unconfined was used as the base simulation 

for all TMR and LGR simulations developed during this study. 

The unconfined version of the regional model was used as the 

platform because it does a better job than the confined version 

in representing the conditions of the upper part of the aquifer 

system where groundwater and surface- water exchange is 

important. 

Lastly, one of the objectives for developing the regional 

model was to provide a platform for development of embedded, higher-resolution models that can be used to address 

water-management issues at smaller (local) scales (Feinstein 

and others, in press). The current study is a demonstration of 

how the regional model may be used to that end.

Model-Grid Refinement

The location for the grid refinement investigation was 

chosen to illustrate a setting typical of the Great Lakes Basin 

and was based on three criteria: the inset model should be in 

the near-field portion of the regional model, it should include 

the glacial- and bedrock-aquifer units, and the inset area 

should include some headwater stream locations. The selected 

area was in the central part of Michigan’s Lower Peninsula 

(fig. 1) where there is appreciable glacial material (layers 1–3) 

as well as a thick bedrock aquifer (sandstone of Pennsylvanian 

age, layer 5) (fig. 2). 

Refined-Model Development

Prior to refining the grid, the regional groundwaterflow model was transferred from MODFLOW-2000 (SEAWAT-2000) to MODFLOW-2005 using Groundwater Vistas©.

Grid refinement for the area of interest was done as a twostage process to achieve a cell size suitable for stream-aquifer 

interactions. This was done to balance the improvement in 

model accuracy in the area of interest while minimizing errors 

along the boundaries of the refined model. When very large 

refinement ratios are used, errors in the model solution at the 

boundaries of the refined model may increase (Mehl and Hill, 

2006). To avoid that phenomenon, the refinement was done 

in two stages. In addition, downsizing to an intermediatescale model made the process more computationally efficient. 

Rather than having to solve the entire regional model in conjunction with the target local area, a smaller intermediate area 

made simulation times more manageable. Lastly, the regional 

model simulated the surface-water network using the river 

package. By refining to an intermediate-scale model, it became 

easier to simulate surface-water features using the SFR1 

package in MODFLOW-2005. This was important because 

it allowed for the use of a beta version of MODFLOW-LGR, 

which enabled routing of streamflow between the coupled 

intermediate and local models.

An intermediate-scale model was initially extracted 

from the regional model using traditional TMR techniques. A 

second refinement was done on the intermediate-scale model 

to produce a local-scale model using the LGR technique. 

All physical properties from the calibrated regional model 

except recharge, land-surface elevation, and stream geometry 

were retained. A summary of the physical properties used in 

the intermediate and local models is found in tables 1 and 2, 

respectively. A complete description of the development and 

properties of the regional Lake Michigan Basin groundwaterflow model is documented in Feinstein and others (in press).

The regional model was simulated and calibrated using 

the codes MODFLOW-2000 (Harbaugh and others, 2000) and 

SEAWAT-2000 (Langevin and others, 2003). SEAWAT-2000 

simulates flow of variable-density fluids, which was necessary owing to the presence of saline water and brines in some 

of the deeper bedrock units in the regional model. For the 

intermediate and local models developed during this study, 

the code MODFLOW-LGR (Mehl and Hill, 2006) was used 

for all simulations, although this code does not simulate flow 

of variable-density fluids. However, this was not a problem 

because the salinity of fluids in the bedrock simulated in the 

area of grid refinement is well below the 10,000 milligrams 

per liter (mg/L) threshold for saline water.
4 Implementation of Local Grid Refinement (LGR) for the Lake Michigan Basin Regional Groundwater-Flow Model

Wisconsin

Illinois

Michigan

Ohio Indiana

MekaL hci nagi

Lake Superior

Lake

Erie

Lake

St. Clair Michigan

84°0'0" 87°0'0" 90°0'0"

48°0'0"45°0'0"42°0'0"

0 5 10 2.5 Kilometers

0 5 10 2.5 Miles

0 100 200 50 Miles

0 100 200 50 Kilometers

Lake

Huron

E

X

P

L

A

N

ATIO

N

Stre

a

m

n

etw

ork

L

a

k

e

Mic

hig

a

n

B

a

sin

Local

m

o

d

el e

xte

nt

Intermediate

m

o

d

el e

xte

nt

R

e

gio

n

al m

o

d

el e

xte

nt

Location of geologic section in figure 2

Pumping center

Figure 1. Location and features of the regional, intermediate, and local models.
Model-Grid Refinement 

5

West East

Vertical Exaggeration = 150x

EXPLANATION

Model near-field

Lake Michigan

00

50 100 MILES 25

25 50 100 KILOMETERS

Unit Layer Color

Upper Quaternary 1Middle Quaternary 2Lower Quaternary 3Jurassic Red Beds 4Upper Pennsylvanian 5

Lower Pennsylvanian 6Michigan Formation 7Marshall Sandstone 8Mississippian-Devonian 9Silurian-Devonian 10Silurian-Devonian 11Silurian-Devonian 12Maquoketa Shale 13Sinnipee 14St. Peter 15Prarie-du-Chien/Franconia 16Ironton-Galesville 17Eau Claire 18Upper Mount Simon 19Lower Mount Simon 20

0 1 2 0.5 Miles

5. 0 2 1 0 Kilometers

West East

Regional model geologic cross section through row 260(modified from Feinstein and others, in press) 

Vertical Exaggeration = 150x

Local model geologic cross section through row 196

Figure 2. Geologic cross sections for the regional and local models.
6 Implementation of Local Grid Refinement (LGR) for the Lake Michigan Basin Regional Groundwater-Flow Model

Table 1. Summary of physical properties used for the intermediate model.

[ft, feet; d, day; E, the number is multiplied by 10 to the denoted power]

Thickness 

(ft)

Horizontal hydraulic 

conductivity (ft/d)

Vertical hydraulic 

conductivity (ft/d)

Confined specific 

storage (ft-1)

Specific yield

Layer 1 Minimum 0.2 1.28E+0 2.17E-2 5.49E-6 0.15

Mean 92.7 1.63E+1 4.04E-1 5.49E-6 0.15

Maximum 171.6 1.92E+2 5.52E+0 5.49E-6 0.15

Layer 2 Minimum 0.2 7.94E-1 4.90E-2 5.49E-6 0.15

Mean 67.5 5.85E+0 2.88E-1 5.49E-6 0.15

Maximum 200.0 9.73E+1 3.07E+0 5.49E-6 0.15

1

Layer 3 Minimum 0.2 7.94E-1 4.90E-2 5.49E-6 0.15

Mean 2.6 6.02E+0 3.01E-1 5.49E-6 0.15

Maximum 95.0 9.73E+1 3.07E+0 5.49E-6 0.15

1

Layer 4 Minimum 0.2 5.00E-1 6.59E-4 2.87E-7 0.07

Mean 4.0 4.22E+0 2.35E-1 4.65E-6 0.14

Maximum 64.3 9.63E+1 3.07E+0 5.49E-6 0.15

Layer 5 Minimum 0.2 1.01E+0 8.16E-2 2.87E-7 0.07

Mean 213.6 9.76E+0 1.03E-1 3.39E-7 0.07

Maximum 371.4 2.60E+1 1.67E+0 5.49E-6 0.15

Layer 6 Minimum 41.6 3.59E-1 1.82E-4 2.87E-7 0.07

Mean 107.8 4.63E+0 1.79E-3 2.87E-7 0.07

Maximum 190.6 5.76E+0 2.85E-3 2.87E-7 0.07

1

Layers 3 and 4 are absent in areas of the model so the physical properties reflect properties of the layer above.

Table 2. Summary of physical properties used for the local model.

[ft, feet; d, day; E, the number is multiplied by 10 to the denoted power]

Thickness 

(ft)

Horizontal hydraulic 

conductivity (ft/d)

Vertical hydraulic 

conductivity (ft/d)

Confined specific 

storage (ft-1)

Specific yield

Layer 1 Minimum 39.0 3.11E+0 4.78E-2 5.49E-6 0.15

Mean 101.7 1.19E+1 2.98E-1 5.49E-6 0.15

Maximum 135.5 6.17E+1 1.78E+0 5.49E-6 0.15

Layer 2 Minimum 50.2 1.06E+0 8.57E-2 5.49E-6 0.15

Mean 154.7 4.73E+0 3.84E-1 5.49E-6 0.15

Maximum 200.0 3.78E+1 3.07E+0 5.49E-6 0.15

1

Layer 3 Minimum 0.2 1.06E+0 8.57E-2 5.49E-6 0.15

Mean 7.6 5.16E+0 4.19E-1 5.49E-6 0.15

Maximum 58.7 3.78E+1 3.07E+0 5.49E-6 0.15

1

Layer 4 Minimum 0.2 1.06E+0 8.57E-2 5.49E-6 0.15

Mean 0.3 5.16E+0 4.19E-1 5.49E-6 0.15

Maximum 5.9 3.78E+1 3.07E+0 5.49E-6 0.15

Layer 5 Minimum 44.5 4.98E+0 1.00E-1 2.87E-7 0.07

Mean 140.0 4.98E+0 1.00E-1 2.87E-7 0.07

Maximum 234.6 4.98E+0 1.00E-1 2.87E-7 0.07

Layer 6 Minimum 73.3 4.00E+0 1.21E-3 2.87E-7 0.07

Mean 96.9 5.06E+0 2.10E-3 2.87E-7 0.07

Maximum 113.1 5.76E+0 2.81E-3 2.87E-7 0.07

1

Layers 3 and 4 are absent in areas of the model so the physical properties reflect properties of the layer above.
Model-Grid Refinement 7

Construction of Intermediate and Local Model 

Grids

An initial grid-refinement routine was used to develop 

an intermediate-scale model from the Lake Michigan Basin 

regional model. To assist in modifying the files associated with 

MODFLOW, the graphical pre- and post-processing program, 

Groundwater Vistas©, was used. Within Groundwater Vistas©, 

an automated TMR process is available (Rumbaugh and Rumbaugh, 2007) for extracting a refined portion of the regional 

model. A 453.8 mi2 area was extracted from the regional 

model and discretized into 230 rows and 220 columns. The 

resulting intermediate-scale model had cell dimensions of 

500 ft per side, or 1/10th the side dimensions of the cells in the 

regional model. Only the top six layers and properties in the 

regional model were translated directly to the intermediate 

(and local) models. All layers were simulated using the convertible layer type 3 option in the layer property flow (LPF) 

package (Harbaugh, 2005). However, the layer top and bottom 

elevations were interpolated between regional model cells so 

that the cell boundaries did not have sharp, blocky transitions 

at the locations of the original regional cell boundaries.

The elevation of the top of layer 1 of both the intermediate and local models was handled differently than the elevations of other layer surfaces. This top surface was recalculated 

for the refined grids by taking the average 30-m Digital Elevation Model (DEM) (U.S. Geological Survey, 2006) elevation 

for each cell and assigning that as the top elevation within the 

geospatial data-processing program, ArcMap (Environmental Systems Research Institute, Inc., 2008). This was done 

because the top of layer 1 elevation in the regional model used 

the average DEM elevation for the 5,000 ft cell, which appreciably smoothed the surface topography. The refined cell size 

of the intermediate model allowed for improved resolution of 

the land-surface features.

Only layers 1–6 were used in the simulations for the 

intermediate- and local-scale models because the hypothetical 

scenarios in these models tested only near-surface processes. 

Keeping the deeper layers in the regional model (layers 7–20) 

would make the evaluation of the head solution computationally expensive, leading to longer simulation run times. 

Additionally, in those deeper layers, the salinity of the groundwater at depth adds complexity to the simulation of flow for 

this area. There was no additional vertical refinement of the 

regional model layering for the intermediate- and local-scale 

models.

The intermediate model that was generated from the 

TMR of the regional model had a cell size of 500 ft per side. 

Although this increased the resolution of the model, an even 

smaller cell size was required for investigation of groundwater/surface-water interactions. Therefore, a local-scale 

model with a cell size of approximately 71.43 ft per side was 

generated using Groundwater Vistas©. The local-scale model 

covered an area (fig. 1) of approximately 21.7 mi2 and was 

composed of 350 rows by 350 columns. Again, in this case, 

the top and bottom elevation of the layers of the model were 

interpolated from elevations used in the intermediate-scale 

model to avoid coarse transitions in cell geometry between 

locations of the former intermediate cell boundaries. As 

before, the 30-m DEM was used to assign elevation to the top 

of layer 1. However, in this case, a nearest-neighbor algorithm 

within ArcMap (Environmental Systems Research Institute, 

Inc., 2008) was used to assign the DEM elevation to the model 

cell in layer 1. This algorithm was used because the model cell 

size was smaller than the DEM cell size. 

The local-model cell-size of 71.43 ft per side was chosen 

to be a 7-to-1 ratio from the local to intermediate model-cell 

dimensions. This ratio of cell dimensions was chosen so that 

the LGR process for MODFLOW-2005 (Harbaugh, 2005) 

could be used. An odd number ratio is required for a cell center to be shared correctly between the parent and child models 

(Mehl and Hill, 2006); in this case, the intermediate and 

the local model. The dimensions of the local grid were then 

adjusted to 344 rows by 344 columns so that only the shared 

centers of the local model align with the local model region of 

the intermediate model (fig. 3).

Boundary Conditions of the Refined Grids

A specified-head-boundary condition was used on the 

external cells of the intermediate model. A no-flow boundary was used at the bottom of the model. This boundary was 

selected at the base of layer 6 because layers 6 and 7 serve as 

confining units in the regional model. Likewise, the groundwater begins to get saline below layer 6. Because of the density 

difference in the water, as well as the low hydraulic conductivity of layers 6 and 7, very little flow was expected upward 

into layer 6. The values assigned to each boundary cell were 

taken from the head solution of the calibrated regional model. 

Specifically, the head solution used was from the final stress 

period of the regional model, representing the period 2001–05. 

For the TMR process, a linear interpolation of the regional 

model heads was used in assigning the specified-head boundary to the intermediate model cells along the outer boundary. 

This was done to avoid a head distribution with sharp contrasts 

along the cell boundaries that overlie the previous regional 

model cell boundaries. 

Recharge

Recharge estimates used for the intermediate- and localscale models were derived from a modified ThornthwaiteMather soil water balance (SWB) code (Westenbroek and others, 2010). The SWB code estimates recharge for an area using 

precipitation, temperature, land-use classification, hydrologic 

soil type, land-surface slope, and soil-water capacity data. The 

initial estimates of these inputs were used to generate an estimate of recharge used in the regional flow model (Feinstein 

and others, in press). The recharge estimates developed from 

the SWB code were compared to recharge estimates developed from base-flow analysis of gaged watersheds within the 
8 Implementation of Local Grid Refinement (LGR) for the Lake Michigan Basin Regional Groundwater-Flow Model

Figure 3. Example of shared-node distribution between the intermediate and local models.

E

X

P

L

A

N

ATIO

N

Intermediate grid

Intermediate nodeLocal node

Shared node

0 5 10 2.5 Miles0 5 10 2.5 Kilometers

Local grid
Model-Grid Refinement 9

regional study area. For the regional model, the SWB inputs 

were then adjusted to produce a recharge result that better 

matched the base-flow-derived recharge estimates. With the 

exception of temperature and precipitation, input parameters 

to the SWB code for the intermediate and local models were 

identical to those used in calculating recharge for the regional 

model. Temperature and precipitation data sets for the intermediate- and local-scale models were obtained from an atmosphere-ocean coupled general circulation model (AOGCM) 

(Hayhoe and others, 2008; Maurer and others, 2002). 

The AOGCM used provided potential temperature and 

precipitation realizations on a daily basis for the period 1960–

2099. The data output from the AOGCM reflects potential 

conditions assuming a high, global greenhouse-gas-emission 

scenario (A1fi scenario) (Hayhoe and others, 2008). Using 

the climatic data derived from the AOGCM as the input to the 

SWB model provided a means to estimate future recharge and 

the effects that global climate change may have on recharge. 

Likewise, the effect of potential climate change on a groundwater system can be examined by using the estimated recharge 

for simulations on the local and intermediate models. To that 

end, one of many possible long-term climate scenarios was 

developed using AOGCM-derived recharge estimates from 

the period 1991–2044. Thus, climate-change results should be 

viewed as illustrative of the approach taken and what could 

happen in a temperate watershed, but should not be considered 

a definitive forecast.

The recharge results obtained by using the AOGCM input 

to the SWB model were interpolated to both the intermediate- and local-model scales using a b-spline method (Stephen 

M. Westenbroek, U.S. Geological Survey, personal commun., 2009). A series of monthly recharge realizations were 

developed from the SWB model for the period 1991–2044. 

The recharge results were output as a series of grids that 

reflected the estimated recharge for a given month. These 

monthly recharge grids then were combined and translated 

into a MODFLOW-2005 recharge package (Harbaugh, 2005). 

The initial calculation of recharge estimates, based on the 

AOGCM, was appreciably lower than current-condition average recharge used in the calibrated regional model for the 

interval over which the two sets of estimates overlap. For the 

time period 1991–2005, the AOGCM-derived recharge estimates for the local model averaged 4.82 in/yr. Recharge values 

estimated in the calibrated regional model over the same 

period for the local model area averaged 7.81 in/yr. To maintain continuity between the regional, intermediate, and local 

models, recharge estimates derived from the AOGCM model 

for the entire period (1991–2044) were scaled up by a factor of 

1.622. This adjustment raised recharge estimates derived from 

the AOGCM inputs to an average of 7.81 in/yr over the local 

model area. Although this adjustment changes the magnitude 

of the recharge, it did not affect the seasonal distribution of the 

recharge given by the AOGCM model (fig. 4). 

Surface-Water Network

Refinement of the regional model allowed for a more 

accurate representation of the surface-water network in the 

intermediate and local models. A comparison of the same 

surface-water network over a common area among the three 

scales of models (fig. 5) illustrates the effect of cell size on 

the surface-water representation. For the local-model area, 

roughly 64 percent of the cells in the regional model include 

surface-water features; in the intermediate model, approximately 14 percent of the cells include surface-water features; 

and for the local model, 2.3 percent of the cells include 

surface-water features. For these reasons, it was necessary 

to reexamine the surface-water network developed for the 

regional model and adapt it for use with the refined model 

grids.

There are many surface-water features within the 

intermediate- and local-model areas. The intermediate model 

includes surface-water features from four different 8-digit 

hydrologic unit codes (HUCs) that drain into Lake Michigan (fig. 6). (Surface-water drainage generally is to the west 

toward Lake Michigan.) There are approximately 170 streams 

and 48 surface-water bodies (lakes and wetlands) within the 

intermediate-model area; there are16 streams but no other 

surface-water bodies within the local-model area.

The source of the surface-water features used for the 

refined models was the 1:100,000 NHDPlus dataset (Horizon 

Systems Corporation, 2006). The NHDPlus dataset is a series 

of geospatial datasets built on the framework of the National 

Hydrography Dataset (NHD). Results of this compilation and 

analysis yielded many useful products, some of which include 

an improved 1:100,000 NHD; a set of value-added attributes 

to enhance stream-network navigation and analysis, flow 

direction, flow accumulation, and elevation grids; and flow 

line min/max elevations and slopes (Horizon Systems Corporation, 2006). These enhancements to the NHD (provided by 

the NHDPlus effort) provided a routed-streamflow network, 

water-feature elevations, and water-feature geometries, which 

were needed to develop the model stress packages that simulate surface-water features in the intermediate- and local-scale 

models.

Simulation of surface-water features other than streams 

(i.e., lakes and wetlands) within the intermediate model (none 

are present in the local model) was done using the generalhead boundary (GHB) package (Harbaugh, 2005). Required 

data to simulate the surface-water features other than streams 

within the intermediate-model area were extracted from the 

NHDPlus dataset within ArcMap (Environmental Systems 

Research Institute, Inc., 2008). The stages for the individual 

water features were based on the average DEM elevation 

within the water body. A hydraulic-conductivity value of 

0.015 ft/d was assigned to the bed material of surface-water 

bodies other than streams based on values from the literature 

(Calver, 2001). An arbitrary bed thickness of 0.5 ft was used 

for all water bodies. This information was exported to a text 

file and edited in a text editor to create the GHB package.
10 Implementation of Local Grid Refinement (LGR) for the Lake Michigan Basin Regional Groundwater-Flow Model

x

x

ox

o

o

o

o

o o

oo o o xx

o

o

o

o

x

x

x

x

o

0 2 4 6 8

56 56 56 56 56 56 56 56 56 56 56 56

x

x

ox

o

o

o

o

o

o o

ooo o oo

xx

o

o

o

o

x

x

x

x

o

Jan

Feb

Mar

Apr

May

Jun

Jul

AugSep

Oct

Nov

Dec

0 2 4 6 8

56 56 56 56 56 56 56 56 56 56 56 56

1990−2045

INCHES

EXPLANATION 

x

o

Data value between 

1.5 times the interquartile 

range and the last quartile

75th percentile

median

25th percentile

Outlier data value 

greater than 3 times 

the interquartile range 

Data value between 

1.5 and 3 times the

 interquartile range

INCHES

56 Number of samples

NON-ADJUSTED

ADJUSTED

Figure 4. Seasonal distribution of recharge for the adjusted and non-adjusted recharge values.
Model-Grid Refinement 11

Figure 5. Surface-water network for the area of the local model shown in comparison to how it was discretized at the levels of grid 

refinement used in the regional and intermediate models. 

42 ̊55’30”

42 ̊58’30”

42 ̊57’00”

84 ̊52’30” 84 ̊51’00” 84 ̊49’30”

0 1 2 Miles

0 1 2 Kilometers

42 ̊55’30”

42 ̊58’30”

42 ̊57’00”

84 ̊52’30” 84 ̊51’00” 84 ̊49’30”

0 1 2 Miles

0 1 2 Kilometers

42 ̊55’30”

42 ̊58’30”

42 ̊57’00”

84 ̊52’30” 84 ̊51’00” 84 ̊49’30”

0 1 2 Miles

0 1 2 Kilometers

Intermediate model 

Local model 

Regional model 
12 Implementation of Local Grid Refinement (LGR) for the Lake Michigan Basin Regional Groundwater-Flow Model

Figure 6. Stream network for the intermediate- and local-model areas.

E

X

PLAN

ATIO

N

8-digit hydrologic unit code (H

U

C)

040500060405000504050007

04050004

LakeLocal model area

Streams

0 5 10 2.5 Miles0 5 10 2.5 Kilometers

84°40'0" 84°50'0" 85°0'0"

43°0'0"

42°50'0"

Intermediate

model area
Refined-Model Scenarios 13

Simulation of streams within the intermediate and local 

models was done using the SFR1 package (Prudic and others, 

2004). This package was chosen because it provides a more 

detailed account of surface-water/groundwater interactions 

than the river package (Harbaugh, 2005). The SFR1 package 

limits water lost from a stream to that which is carried by the 

stream, rather than using a set stream stage to estimate water 

lost from a losing reach, and the inputs needed to generate 

the package generally were available within the NHDPlus 

dataset. These required inputs included an ordered stream 

network, location, stream elevation, streambed thickness, 

streambed hydraulic conductivity, and stream geometry. After 

these inputs were compiled, a program add-on for ArcMap 

(Steven Predmore, ESRI, personal commun., 2009) was used 

to process the inputs and generate a formatted SFR1 file. A 

hydraulic-conductivity value of 0.283 ft/d was assigned to 

the bed material based on values from the literature (Calver, 

2001). As above, an arbitrary bed thickness of 0.5 ft was used 

for all streams. Further detail on the construction of the stream 

network used for the beta version of MODFLOW–LGR is 

described in appendix 1.

Refined-Model Scenarios

Five model scenarios were developed to demonstrate the 

utility of grid refinement on modeling the effect of groundwater withdrawals and climate change on streamflow for 

the regional LMB model. All scenarios were hypothetical 

and built on the hydrogeologic framework of the calibrated 

regional model. There were no attempts to calibrate or execute 

parameter estimation on any of the following scenarios. 

Streamflow Depletion

A series of scenarios was developed to demonstrate the 

effect of groundwater withdrawals on a headwater stream. 

Conceptually, this was done by placing a well near a headwater stream and evaluating the effect of pumping on flow in the 

stream (fig. 7). All scenarios were transient simulations using 

monthly stress periods with five time steps per stress period. 

Current (2010) condition simulations were run for 9 years 

(1991–99); these simulations used a 4-year run-up (steadystate) period to allow conditions to equilibrate. Each scenario 

used the modified recharge estimates developed from the SWB 

code discussed previously. 

To achieve steady-state conditions prior to the transient 

simulation, the model was simulated for 3 years (36 stress 

periods) using constant recharge rates applicable to January 1990. Because the model had not fully equilibrated after 

the 3-year run-up, the model run-up period was extended to 

include 1990, and the evaluation of model results was done 

only for the period 1991–99. 

In addition, a monthly flow value had to be estimated for 

each stream that originated outside the intermediate-model 

area. For those streams, the precipitation data used for the 

SWB (Westenbroek and others, 2010) code to calculate 

recharge was used as the input to the AFINCH (Holtschlag, 

2009) code to derive the appropriate monthly streamflow for 

that segment.

Four scenarios were developed to demonstrate aquiferstream interaction and response to different screen depths 

of the pumped well near the stream gage (fig. 7). The first 

scenario was a control case with no pumping simulated near 

the stream. The second scenario simulated a well that was 

screened in the glacial deposits of layer 1 at approximately 

95 ft below ground surface. The third scenario simulated a 

pumping well that was screened in the deeper glacial deposits of layer 2, approximately 280 ft below land surface. The 

fourth scenario placed the pumping well open to the Pennsylvanian sandstone of layer 5 at approximately 405 ft below 

land surface. In each pumping scenario, the well was assigned 

a discharge rate of 71.4 gal/min and was located approximately 460 ft east of the stream (fig. 7). Pumping was initiated 

in the first stress period and remained constant throughout the 

simulation.

To demonstrate the effects of pumping on streamflow, the 

Gage package (Prudic and others, 2004) was used to record 

the streamflow at a point near the well (fig. 8). The streamflow 

for the non-pumping control case was treated as an unstressed 

base case. Streamflows from the three other scenarios were 

compared to the streamflow from the non-pumping case. The 

hydrograph in figure 8 indicates that the stress on the stream 

is greater in the low-flow months during the late summer, fall 

and winter seasons owing to a combination of low recharge 

and continued pumping from the nearby well. During those 

periods, streamflow was depleted up to 100 percent, causing 

the stream to dry up. It should be noted that in some dry years 

the stream dried up in the control case, which was not affected 

by pumping, so during some periods there was no streamflow 

available for depletion.

In all cases, pumping had the effect of decreasing streamflow in the test stream, which is expected from conservation 

of mass (fig. 8). Although the amount of depletion among the 

various scenarios was similar, pumping from layer 1 had the 

largest effect on streamflow in the test stream. Pumping from 

layer 2 and pumping from layer 5 had slightly less effect on 

the streamflow than pumping from layer 1 because the wells 

in these scenarios were placed in deeper layers farther from 

the stream. Table 3 provides a summary of the cumulative 

water budget for pumping from layer 1. This table illustrates 

the relative magnitude of the various components of the water 

budget in the test watershed. To better discern the effect 

of relocating the pumping well to deeper parts of the flow 

system, a mass-balance approach was used to determine the 

source of the water to the pumping well. This was done by 

compiling the water-budget information for the test watershed 

for all four scenarios. The sources of water to the pumping 

well were then derived by subtracting the water budget for 

the pumping scenario from the water budget of the naturalconditions scenario. Recharge was the same for all scenarios, 
14 Implementation of Local Grid Refinement (LGR) for the Lake Michigan Basin Regional Groundwater-Flow Model

Streamgage Well

0 1 2 0.5 Miles

0 1 2 0.5 Kilom eters

Test watershedStream sTest stream

EXPLAN

ATION

Figure 7. Location of test watershed, test stream, stream gage, and pumped well for refined-model scenarios.
Refined-Model Scenarios 15 Figure 8. Hydrographs of test stream for streamflow-depletion scenarios.

Jan91 Jan92 Jan93 Jan94 Jan95 Jan96 Jan97 Jan98 Jan99 Jan00

0

0.5

1

1.5

2

Flow, in cubic 

feet per second

Pumping in layer 1

Streamflow

Streamflow with pumping

Jan91 Jan92 Jan93 Jan94 Jan95 Jan96 Jan97 Jan98 Jan99 Jan00

DATE

Reduction of flow, 

in cubic feet per second

Jan91 Jan92 Jan93 Jan94 Jan95 Jan96 Jan97 Jan98 Jan99 Jan00

0

0.5

1

1.5

2

Flow, in cubic 

feet per second

Pumping in layer 2

Jan91 Jan92 Jan93 Jan94 Jan95 Jan96 Jan97 Jan98 Jan99 Jan00

DATE

StreamflowStreamflow with pumping

Jan91 Jan92 Jan93 Jan94 Jan95 Jan96 Jan97 Jan98 Jan99 Jan00

0

0.5

1

1.5

2

Flow, in cubic 

feet per second

Pumping in layer 5

Jan91 Jan92 Jan93 Jan94 Jan95 Jan96 Jan97 Jan98 Jan99 Jan00

DATE

StreamflowStreamflow with pumping

0

0.05

0.1

0.15

0.2

0

0.05

0.1

0.15

0.2

Reduction of flow, 

in cubic feet per second

0.05

0.1

0.15

0.2

0

Reduction of flow, 

in cubic feet per second
16 Implementation of Local Grid Refinement (LGR) for the Lake Michigan Basin Regional Groundwater-Flow Model

so that term cancelled out in the analysis. The water-budget 

components of interest became streamflow, storage, and lateral 

flow (water derived from outside the test watershed). Results 

of this analysis are illustrated in figure 9. 

In all scenarios, the source of water to the well varied 

annually. During dry years, such as 1993, more water pumped 

from the well was derived from lateral groundwater flow and 

not from the stream and storage (fig. 9) because there was 

not as much water available in the stream (fig. 8). During a 

wet year, like 1997, much of the water entering the well was 

coming from the stream. The annual variation in sources of 

water to the well was similar in all cases, yet there were small 

differences in the magnitude of the various sources depending on the layer in which the well was located. For instance, 

in 1997 the well pumping from layer 1 captured 62 gal/min 

from streamflow, the well pumping from layer 2 captured 55 

gal/min from streamflow, and the well pumping from layer 5 

captured 52 gal/min from streamflow. These results indicate 

that deepening the well changes the distribution of the sources 

of water to the well. 

The effect of deepening the well was more apparent when 

comparing the net source of water to the well for the entire 

period (1991–99) (fig. 9). In all cases, the majority of the 

water to the well came from streamflow. As the pumping well 

was placed into deeper geologic units, the effect on streamflow 

was reduced, largely owing to the increased contribution of 

lateral flow of water from induced recharge from surrounding 

watersheds. This finding implies that for this hydrogeologic 

setting, to lessen the effect of pumping on a stream, a well can 

be placed in a deeper geologic unit, causing the well to capture 

the lateral flow and, perhaps, to spread its effect to other 

streams and (or) other watersheds.

Climate-Change Scenario

The last scenario developed for the refined models was 

a scenario examining the effects of potential climate change 

on streamflow (scenario 5). Similar to the streamflow-depletion scenarios, the climate-change scenario was a transient 

simulation that used a monthly stress period with five time 

steps in each stress period. This scenario simulated a period 

of 54 years (1991–2044), again with a 4-year run-up (steadystate) period to allow the model to equilibrate. 

The Gage package (Prudic and others, 2004) was used to 

monitor streamflow at the same location as in the streamflowdepletion scenarios (fig. 7). In this scenario, climate change 

was simulated by variable areal recharge. The annual distribution of recharge for the simulated period is illustrated in figure 

10. SWB-derived recharge was quite variable over the period 

simulated, ranging from 17.62 in/yr in 2011 to 3.38 in/yr in 

2042. There was a slight increase in recharge until 2017 and 

then it generally decreased until the end of the simulation in 

2044 (illustrated by the LOESS-smoothed curve in figure 10). 

This variation in recharge was reflected in streamflow during 

the simulation; however, this is only one realization of possible climate and should be considered illustrative of the type 

of system response that could potentially occur with climate 

change. 

Simulated streamflow was greater during periods of high 

recharge and less during periods of low recharge (fig. 10). 

The period 2013–2022 had the greatest streamflows, which 

correlated to a period of higher recharge (fig. 10). The period 

2035–2044 generally had the lowest recharge values, and 

that was reflected in low streamflows for that period (fig. 10). 

During this lower recharge period, the median flow for several 

months was zero, indicating that this stream was frequently 

dry.

Generally, the monthly distribution of streamflow was 

similar for all periods, with the highest flows during March 

and April and the lowest flows during October–December. 

However, the recharge distribution by month is slightly different (fig. 4). The highest recharge occurred during December–April, with little to no recharge taking place from May to 

September. In contrast, the lowest streamflows (fig.10) were 

observed during October–December, which indicates a lag 

time between the recharge decline and the ensuing decline in 

streamflow. 

Table 3. Summary of the overall water budget in the test watershed with a well pumping in layer 1.

[gal/min, gallons per minute]

Flux into the test watershed

(gal/min)

Flux out of the test watershed

(gal/min)

Component Amount Component Amount

Recharge 108,253 Boundary of watershed 74,711 

Stream 3,650 Stream 27,740 

Storage 57,698 Well 7,710 

Storage 59,441 
Refined-Model Scenarios 17

Well pumps 71.4 gallons per minute from layer 1

Stream sourceStorage sourceLateral-flow sourceSum sources

Well pumps 71.4 gallons per minute from layer 2

-28-14

0

83

69

56422814

1991 1992 1993 1994 1995 1996 1997 1998 1999

1991 1992 1993 1994 1995 1996 1997 1998 1999

GALLONS PER MINUTE

1991 1992 1993 1994 1995 1996 1997 1998 1999

59.11%

53.26%

50.12%

7.65%

6.73%6.51%

33.23%

40.01%43.37%

Well pumps 71.4 gallons per minute from layer 5NET SOURCE OF WATER TO WELL 1991−99

EXPLANATION

ANNUAL SOURCE OF WATER TO WELL 1991−99

-28-14

0

83

69

56422814

-28-14

0

83

69

56422814

Figure 9. Distribution of the source of water to the pumping well for the stream-depletion scenarios.
18 Implementation of Local Grid Refinement (LGR) for the Lake Michigan Basin Regional Groundwater-Flow Model

Figure 10.

xxxxx

Recharge and streamflow distribution shown for three time periods for the test stream in the climate-change scenario.

xxx

oooo

x

xx

xxx

xxx xxxx x

xx

Jan

Feb

Mar

Apr

May

Jun

Jul

Aug

Sep

Oct

Nov

Dec

2.5 0 2. 5 1. 0 1. 5 0. 0 0.

50 50 50 50 51 50 50 50 50 49 50 50

1991−2000

DISCHARGE, IN CUBIC FEET PER SECOND

xxx

Jan

Feb

Mar

Apr

May

Jun

Jul

Aug

Sep

Oct

Nov

Dec

2.5 0 2. 5 1. 0 1. 5 0. 0 0.

50 50 50 50 50 50 50 50 50 50 50 50

2013−2022

DISCHARGE, IN CUBIC FEET PER SECOND

x

Jan

Feb

Mar

Apr

May

Jun

Jul

Aug

Sep

Oct

Nov

Dec

2.5 0 2. 5 1. 0 1. 5 0. 0 0.

50 50 50 50 50 50 50 50 50 50 50 50

2035−2044

DISCHARGE, IN CUBIC FEET PER SECOND EXPLANATION

o

x

Data value between 1.5 times the interquartile range and the last quartile

75th percentile

median

25th percentile

Outlier data value greater than 3 times the interquartile range 

Data value between 1.5 and 3 times the

 interquartile range

50 Number of samples

1991 2001 2011 2021 2031 2041

YEAR

50

1510

ANNUAL RECHARGE, IN INCHES

Annual recharge for inset model 1991−2044

Annual rechargeLOESS fit
Summary 19

simulated. The precipitation-runoff component of streamflow 

was not simulated in the groundwater-flow model, so that 

flow in streams was derived either from existing flow (at the 

boundaries of the intermediate model) or from discharge from 

the aquifer (base flow). Finally, the climate scenario used for 

this investigation should be viewed as one possible realization of the potential effects of climate change. There are many 

other climate scenarios available, which illustrate an array of 

potential climate-change outcomes that could have been used 

for this investigation. The focus of this investigation was to 

demonstrate a tool for assessing the potential effects of climate 

change on a hydrologic system, rather than to provide actual 

predictions of climate change. 

Summary

The U.S. Geological Survey is providing an assessment of the availability and use of water resources, and the 

factors that affect them, throughout the United States. The 

Great Lakes Basin was chosen as the pilot study area to test 

new techniques and methods for evaluating water resources. 

These methods and techniques will then be transferred to other 

water-availability assessment areas. As part the Great Lakes 

Basin pilot study, a regional groundwater-flow model was 

constructed for the Lake Michigan Basin. This regional model 

was used as the framework for developing refined groundwater-flow models using telescopic mesh refinement and local 

grid refinement techniques.

The refined models developed illustrate the capability 

to use regional models to assess local-scale problems. The 

refined-model properties were taken directly from the regional 

model; however, only the upper six layers of the regional 

model were used for the refined-model analysis and were 

not refined further. Two separate lateral refinements were 

executed, creating models in an area of central Michigan, to 

ensure adequate refinement for headwater-stream simulation 

and reduce the computational burden: an intermediate-scale 

model with 500 feet by 500 feet cells covering an area of 

454 square miles, and a local-scale model with 71.4 feet by 

71.4 feet cells covering an area of 21.7 square miles. 

Recharge for the refined models was estimated using the 

soil water balance (SWB) approach to provide initial recharge 

estimates for the regional model. Instead of using the same 

climate variables used for the regional model, climate inputs 

from an atmosphere-ocean coupled general circulation model 

(AOGCM) simulating high greenhouse-gas- emission conditions were used by the SWB code to estimate recharge. 

The surface-water network used for the refined models 

was developed based on the National Hydrography Dataset 

(NHDPlus) datasets. Surface-water bodies were simulated 

using the general-head boundary package, and streams were 

simulated using a streamflow-routing package. The NHDPlus data were processed to compile streambed elevation and 

geometry into one data set. This data set was processed further 

Assumptions and Limitations of the 

Models

There were several assumptions made for the purposes 

of this investigation. For instance, it was assumed that the 

streambed for all streams in the model area had constant thickness and constant hydraulic conductivity. This was done to 

facilitate the demonstration of the modeling methods used, but 

these characteristics could be refined using parameter estimation and field verification for actual case studies. It also was 

assumed that all recharge estimated using the SWB code for a 

dedicated stress period immediately infiltrated the water table 

(no unsaturated-zone flow or processes were simulated).

The assumption that the bottom of layer 6 acts as a 

no-flow boundary may not be adequate for simulations that 

take place over very long time intervals (hundreds to thousands of years). There is a low potential for flow through the 

bottom of layer 6, but over long time periods, the volume of 

flow into that layer may be appreciable. Likewise, the salinity conditions that also minimize vertical flow into and out 

of layer 6 may change over long periods of time. Therefore, 

using a no-flow boundary at the base of layer 6 may have 

to be reevaluated for simulations of very long time periods. 

Another assumption used in the model development was that 

the specified-head boundaries in the intermediate model were 

adequately representative of the current conditions and static 

during the model simulations. Although this may be more 

accurate during the 9-year stream-depletion scenarios, it poses 

a problem for the long-term climate scenarios. Similarly, the 

specified- head boundaries used in the intermediate model 

have the potential to affect heads and flows near the pumping 

well in the local model. However, the low pumping rate of the 

well in conjunction with the large distance (6.4 mi) of the well 

from the boundaries makes this unlikely but it was not tested 

in these scenario models. Additionally, for the stream depletion scenarios, the gage location used for recording streamflow 

hydrographs does not reflect the holistic effect of the pumping 

well on the stream reach. However, the output from the gage 

used for these scenarios does illustrate the effect of a well 

pumping near a stream.

One limitation of the refined models in this study was 

that there was no vertical refinement of any of the model 

layers. However, future studies that may require a finer vertical discretization, such as for simulated particle tracking, 

could utilize this feature of LGR. The SWB results used in 

the intermediate and local models did not take into account 

the effect of stomatal resistance of plants and the variation 

in stomatal resistance, which may affect evapotranspiration 

and ultimately groundwater recharge. The SWB code also 

did not take into account potential changes in future vegetation or land-cover types, so the current land use was assumed 

static throughout the period of interest analyzed by the SWB 

code. In addition, routing of overland flow from precipitation was not accounted for in the SWB code; therefore, any 

effects this may have on the distribution of recharge are not 
20 Implementation of Local Grid Refinement (LGR) for the Lake Michigan Basin Regional Groundwater-Flow Model

to limit the extent of the network to only the intermediate and 

local models. A program was used to assign starting flows to 

streams that originated from outside the intermediate- model 

boundary. Because stream-width data were not available for 

every stream segment, an option was used that uses a power 

function to relate streamflow to stream depth and width when 

calculating the geometry of the stream. The refined models 

then were used to demonstrate stream-depletion scenarios and 

a potential climate-change scenario. Demonstrating stream 

depletion was done using four scenarios: a natural-conditions 

scenario with no pumping, a scenario with pumping from 

layer 1, a scenario with pumping from layer 2, and a scenario 

with pumping from layer 5. Analyses of the scenarios indicate 

that pumping decreased streamflow in the test stream, with 

the wells pumping from the deeper layers having less of an 

effect than the well that was pumping from layer 1 and that 

was adjacent to the stream. Examining the scenarios through 

a mass-balance approach indicated that the net source of the 

water to the well ranged from 59 percent streamflow for the 

well in layer 1 to 50 percent streamflow for the well in layer 5. 

The difference largely is made up by more water taken from 

surrounding watersheds (lateral flow) as the well was placed 

deeper.

A climate scenario was developed to provide one realization of effects of potential climate change on streamflows. This 

scenario used the recharge developed from the AOGCM inputs 

to the SWB model for the period 1991–2044 and was adjusted 

so that the current (1991-2005) period was more representative 

of actual climate. Results of the scenario indicate that recharge 

and related streamflow during the model-simulation period 

trended slightly upward toward the middle of the simulation 

period and trended downward from the middle to the end. 

These scenarios illustrate the utility of taking a regional 

model and applying grid- refinement techniques to evaluate 

local-scale problems. The grid- refinement approach couples 

the local-scale to the larger-scale framework, thus helping 

ensure that boundary conditions are consistent and reasonable. 

Local-scale problems can be evaluated, such as determining 

whether enough streamflow can be maintained in a stream 

affected by pumping to support the aquatic ecosystem it supported prior to pumping. Likewise, competing uses among 

various water users can be better evaluated at a smaller scale. 

Using these various modeling tools that build upon an existing 

modeling framework means that many local-scale waterresource issues can be efficiently assessed by leveraging existing groundwater-flow models. 
References Cited 21

References Cited

Arihood, L.D., 2009, Processing, analysis, and general evaluation of well-driller records for estimating hydrogeologic 

parameters of the glacial sediments in a ground-water flow 

model of the Lake Michigan Basin: Scientific Investigations 

Report 2008–5184, 26 p.

Calver, Ann, 2001, Riverbed permeabilities—information 

from pooled data: Ground Water, v. 39, no. 4, p. 546–553.

Environmental Systems Research Institute, Inc., 2008, 

ArcMap, GIS software, Version 9.2: Redlands, Calif.

Feinstein, D.T., Hunt, R.J., and Reeves, H.W., in press, 

Regional ground-water-flow model of the Lake Michigan 

Basin in support of Great Lakes Basin water availability 

and use studies: U.S. Geological Survey Scientific 

Investigations Report 2010–5109.

Grannemann, N.G., and Reeves, H.W., 2005, Great Lakes 

Basin water availability and use–A study of the National 

Assessment of Water Availability and Use Program: U.S. 

Geological Survey Fact Sheet 2005–3113, 4 p.

Harbaugh, A.W., Banta, E.R., Hill, M.C., and McDonald, 

M.G., 2000, MODFLOW-2000, the U.S. Geological Survey 

modular ground-water model—User guide to modularization concepts and the Ground-Water Flow Process: U.S. 

Geological Survey Open-File Report 00–92, 121 p.

Harbaugh, A.W., 2005, MODFLOW-2005—The U.S. Geological Survey modular ground-water model—The groundwater flow process: U.S. Geological Survey Techniques and 

Methods 6–A16 [variously paged].

Hayhoe, K., Wake, C., Anderson, B., Liang, X.-L., Maurer, 

E., Zhu, J., Bradbury, L., DeGaetano, A., Stoner A., and 

Wuebbles, D., 2008, Regional climate change projections 

for the northeast USA: Mitigation and Adaptation Strategies 

for Global Change, v. 13, no. 5–6, p. 425–436.

Holtschlag, D.J., 2009, Application guide for AFINCH 

(analysis of flows in networks of channels) described by 

NHDPlus: U.S. Geological Survey Scientific Investigations 

Report 2009–5188, 106 p.

Horizon Systems Corporation, 2006, NHDPlus Home—

Horizon Systems Corporation NHDPlus database, accessed 

January 9, 2008, at http://www.horizon-systems.com/NHDPlus/.

Lampe, D.C., 2009, Hydrogeologic framework of bedrock 

units and initial salinity distribution for a simulation of 

ground-water flow for the Lake Michigan Basin: U.S. 

Geological Survey Scientific Investigations Report 

2009–5060, 49 p.

Langevin, C.D., Shoemaker, W.B., and Guo, Weixing, 2003, 

MODFLOW-2000, the U.S. Geological Survey modular 

ground-water model—Documentation of the SEAWAT-2000 

version with the variable-density flow process (VDF) and 

the integrated MT3DMS transport process (IMT): U.S. 

Geological Survey Open-File Report 2003–426, 43 p.

Leake, S.A., and Claar, D.V., 1999, Procedure and computer 

programs for telescopic mesh refinement using 

MODFLOW: U.S. Geological Survey Open-File Report 

99–238, 53 p.

Leopold, L.B., Wolman, M.G., and Miller, J.P., 1992, Fluvial 

processes in geomorphology: Mineola, New York, Dover 

Publications, Inc., 522 p.

Maurer, E.P., Wood, A.W., Adam, J.C., Lettenmaier, D.P., and 

Nijssen, B., 2002, A long-term hydrologically based dataset 

of land surface fluxes and states for the conterminous 

United States: Journal of Climate, v. 15, no. 11, 

p. 3237–3251.

Mehl, S.W., and Hill, M.C., 2006, MODFLOW-2005, the U.S. 

Geological Survey modular ground-water model—Documentation of shared node local grid refinement (LGR) and 

the boundary flow and head (BFH) package: U.S. Geological Survey Techniques and Methods 6–A12, 78 p.

Prudic, D.E., Konikow, L.F., and Banta, E.R., 2004, A new 

streamflow-routing (SFR1) package to simulate streamaquifer interaction with MODFLOW-2000: U.S. Geological 

Survey Open-File Report 2004–1042, 104 p.

Rumbaugh, J.O., and Rumbaugh, D.B., 2007, Groundwater 

Vistas, Version 5: Reinholds, Penn., Environmental Simulations, Inc., 375 p., accessed January 14, 2010, at http://www.

groundwatermodels.com/Groundwater_Vistas.php.

The Nature Conservancy, 2008, ELOHA fact sheet, accessed 

June 29, 2009, at http://conserveonline.org/workspaces/

eloha/documents/brochure-english.

U.S. Geological Survey, 2006, National Elevation Dataset (NED) 30-meter Digital Elevation Models (DEMs): 

National Mapping Division, scale 1:24,000, accessed January 14, 2010, at http://gisdata.usgs.net/NED.

Westenbroek, S.M., Kelson, V.A., Dripps, W.R., Hunt, R.J., 

and Bradbury, K.R., 2010, SWB—A modified Thornthwaite-Mather soil-water-balance code for estimating 

groundwater recharge: U.S. Geological Survey Techniques 

and Methods 6–A31, 67 p.

Appendix 1 23

Appendix 1. Surface-Water Network 

This appendix describes the processing steps used to convert stream information in the 

NHDPlus dataset to the SFR1 package format used to simulate streams in MODFLOW-LGR.

Streams in all model scenarios were simulated using the SFR1 package. Streamfeatures 

simulated with the SFR1 package are identified by reach and segment. A reach is defined as 

a portion of a stream that lies within a model cell. A segment is a group of connected reaches 

of a stream that can span multiple model cells (full description in Prudic and others, 2004). 

Segments and reaches must be ordered and numbered sequentially starting from the farthest 

upstream segment and terminating at the farthest downstream segment. The upstream segment 

must have a lower number than the downstream segment for flow to be correctly routed from 

the upstream segment to the downstream segment. The NHDPlus dataset was well suited for 

transitioning to the SFR1 package, because the geographic information system (GIS) stream 

coverage was organized as several different segments (stream arcs) correctly routed from 

upstream to downstream. 

Several geospatial data-processing steps in ArcMap were executed to prepare the NHDPlus stream data for conversion to an SFR1 file. The first step was to join the stream network 

to the value-added attribute and the flow-line attribute tables (Horizon Systems Corporation, 

2006) generated as part of the NHDPlus effort. These tables contain the stream-arc elevations 

and flow-routing information that ensure flow is correctly routed through the stream network. 

Once these data were appended to the stream segments, streams within the intermediatemodel area were extracted from the NHDPlus dataset within ArcMap (Environmental Systems 

Research Institute, Inc., 2008). 

Streams that originated outside the outer boundary of the intermediate model were split at 

the intersection of the model boundary and stream segment. In that case, the information associated with the affected stream segment was updated to reflect the changes. Specifically, a new 

stream-segment length was calculated for the affected segment, and an elevation was calculated 

at the newly created boundary node of the stream segment. The elevation of the stream arc at 

the model boundary was calculated using the slope of the stream arc and the elevation of the 

known end of the stream arc from the NHDPlus data along with the new stream-arc length. 

Flows into each segment that originated outside the external boundary of the intermediate 

model were estimated using the “Analysis of Flows In Networks of Channels” code (AFINCH) 

(Holtschlag, 2009). AFINCH estimates were needed because streamflow for those streams 

was not available from the regional model (and the SFR1 package was not used in the regional 

model). A flow of zero was assigned to all headwater segments. Flows were assigned to segments first by identifying the origin of each segment, which included creating a GIS coverage 

of points that represented either stream inflow or stream diversions (where a stream segment 

splits into two separate segments). All streams that originated outside the intermediate-model 

boundary and flowed into the model grid had to have an inflow point created. Once the flow 

calculated by AFINCH was assigned to the inflow points, a GIS coverage representing the 

model grid was created. The model-grid coverage was used by the ArcMap SFR1 creation program to locate the SFR reaches in the correct model cell and assign the stream elevation to each 

reach. Stream elevations in the NHDPlus data ensure proper flow routing, so these elevations 

were taken from the stream coverage and inserted into the model-grid coverage. This was done 

to make certain that the elevations were consistent between the model grid and the NHDPlus 

data. The ArcMap SFR1 creation tool then was executed, which (1) ordered and numbered all 

the stream segments and reaches according to SFR1 package requirements (Prudic and others, 

2004), (2) assigned a starting flow to each segment that started at the model boundary, and (3) 

translated all the related GIS information into an SFR1 package formatted text file. 

The SFR input file specifies if and by which technique stream stage is calculated using 

the ICALC variable (Prudic and others, 2004). In this study, ICALC was set equal to 3, which 

related streamflow to stream width and depth using a power function relation from Leopold 
24 Implementation of Local Grid Refinement (LGR) for the Lake Michigan Basin Regional Groundwater-Flow Model

and others (1992). The numerical constants CDPTH, FDPTH, AWDTH, and BWDTH (Prudic 

and others, 2004) were estimated using nonlinear regression of known widths and depths at 

five stream-gage locations within or near the intermediate-model area. Results of the regression yielded constants of 0.00083, 0.48, 0.63, and 0.31 for CDPTH, FDPTH, AWDTH, and 

BWDTH respectively.

 A beta version of MODFLOW-LGR was used in this study, which was able to route 

surface-water flow between the intermediate and local models (Steffen W. Mehl, U.S. Geological Survey, personal commun., 2008). To implement this version of MODFLOW-LGR, stream 

segments that cross the boundary of the intermediate and local models had to be split at the 

interface of the two models (fig. 1–1). This was done so that in a fully coupled LGR-simulation, 

flow from streams in the intermediate model would end at a segment at the boundary of the 

local model and be transferred to the analogous stream segment within the local model. Likewise, flow from the stream segments in the local model would end at the model boundary and 

be transmitted back across to the analogous stream segment within the intermediate model. 

The procedure for generating the SFR1 package for the local model was the same as that 

for the intermediate model with the exception that starting flows were not assigned to streams 

that originate outside the local-model boundary because that flow was routed from the intermediate model to the local model. To allow for the proper routing of flow between model grids, 

two variables were added to item 4b of the SFR1 package (Prudic and others, 2004; Steffen W. 

Mehl, U.S. Geological Survey, personal commun., 2008). The first variable, AUX LGRGRID, 

was an integer value indicating from which grid the flow would enter. The second variable, 

AUX LGRSEG, was an integer indicating which stream segment the flow would exit from the 

first grid. For instance, the grid number and the segment number of the stream segment feeding 

into the local-model stream segment were appended to the end of item 4b in the SFR1 file for 

local streams receiving flow from the parent grid (Prudic and others, 2004). For streams not 

receiving flow from the parent model, AUX LGRGRID and AUX LGRSEG were assigned the 

value 0. This process was the same for streams in the parent model receiving flow from the 

local model; however, the stream segments in the parent model that were covered by the localmodel grid were given values of -1 for both AUX LGRGRID and AUX LGRSEG. A -1 value 

for those variables indicates no simulation of flow because these segments are redundant in that 

they are simulated by the local grid in the fully-coupled model.
Appendix 1 25

Figure 1–1. Stream-segment construction across the intermediate and local-model boundary.

0 0.02 0.04 0.08 Miles

0 0.02 0.04 0.08 Kilometers

EXPLANATION

Local grid

Intermediate grid

Intermediate stream

Local stream




Hoard—Implementation of Local Grid Refinement (LGR) for the Lake Michigan Basin Regional Groundwater-Flow Model—Scientific Investigations Report 2010–5117

Printed on recycled paper