GIS-ENABLED MODELING OF MICHIGAN’S GROUNDWATER SYSTEMS

By
Mehmet Oztan

A DISSERTATION
Submitted to
Michigan State University
in partial fulfillment of the requirements
for the degree of
DOCTOR OF PHILOSOPHY
Civil Engineering
2011

ABSTRACT
GIS-ENABLED MODELING OF MICHIGAN’S GROUNDWATER SYSTEMS
By
Mehmet Oztan

In this paper we systematically evaluate a GIS-enabled and data-intensive modeling system to
assess the system’s capability to simulate the regional-scale ground water flow in 26 of 38 USGS
8-digit Hydrologic Unit Codes (HUCs) across Michigan’s Lower Peninsula. All models
simultaneously simulate the glacial and uppermost bedrock aquifer layer in the peninsula.
Specifically, the modeling system is used to simulate the long-term average static water levels
and base flow in the glacial aquifer, and long-term average piezometric heads in the bedrock
aquifer. The results are compared with the measured ground water levels from Michigan’s
statewide ground water database (MSGWD) and USGS-estimated base flow values. The
coupling between the models and data allows real-time and interactive analysis of the
assumptions and boundary conditions used in the conceptual models. Overall, this paper presents
the most extensive regional-scale modeling and evaluation effort in Michigan to date that
provides critical insight into the state’s ground water systems. We postulate that the results of
this research will have significant implications on the sustainable management of Michigan’s
ground water resources and ground water dependent ecosystems.

Keywords: Michigan’s ground water systems, Michigan’s statewide ground water database, GISenabled ground water modeling, regional-scale ground water flow, coupling between data and
model, Interactive Ground Water

MİŞİGIN’IN YERALTISUYU SİSTEMLERİNİN CBS TABANLI MODELLEMESİ

Bu araştırma, CBS tabanlı ve çok sayıda arazi verisine dayalı bir nümerik modelleme sisteminin
Mişigın Eyaleti’nin güney yarımadasında yer alan 26 büyük ölçekli yüzeysuyu havzasındaki
yeraltısuyu akımını temsil etmekteki yeterliliğini değerlendirmektedir. Araştırma için geliştirilen
bütün modeller eşzamanlı olarak havza dahilindeki buzul ve kaya akiferlerini
kavramsallaştırmaktadır. Modelleme sistemi, buzul akiferdeki uzun zamanlı yeraltısuyu
seviyeleri ve baz akım ile kaya akiferindeki piezometrik yükleri hesaplamakta kullanılmıştır.
Sonuçlar Mişigın’ın yeraltısuyu veritabanında yer alan yeraltısuyu seviye ölçümleri ve USGS
tarafından havza ölçeğinde gerçekleştirilmiş baz akım tahminleriyle karşılaştırılmıştır. Sayısal
modelleri oluşturmak için kullanılan kesintisiz veri akışı, kavramsal modelde kullanılan
varsayımların ve sınır koşullarının gerçek zamanlı ve etkileşimli analizini olası kılmaktadır.
Genel olarak, bu araştırma Mişigın Eyaleti’nde günümüze kadar gerçekleştirilmiş en kapsamlı
bölgesel ölçekli modelleme çalışması olup eyaletin yeraltısuyu sistemlerinin işleyişi hakkında
önemli bilgiler sunmaktadır. Araştırma sonuçlarının, gelecekte eyaletin yeraltısuyu kaynakları ile
devamlılığı bu kaynaklara bağlı olan ekosistemlerin sürdürülebilir yönetiminde önemli katkı
sağlayacağı öngörülmektedir.

Anahtar Kelimeler: Mişigın’ın yeraltısuyu sistemleri, Mişigın’ın yeraltısuyu veritabanı, CBS
tabanlı yeraltısuyu modellemesi, bölgesel ölçekli yeraltısuyu akımı, modelleme için kesintisiz
veri akışı, Interactive Ground Water

Copyright by

MEHMET OZTAN

2011

I dedicate this doctoral dissertation to my mother Ruveyde Leyla Oztan,
my father Prof. Sadi Yilmaz Oztan (1930-2003),
my uncle Prof. Yuksel Oztan (1933-2010), and to the graduate students who don’t give up
despite all difficulties they experience through their education
Mehmet Oztan, May 2011
Tampa, Florida

v

Hayatımın en önemli eseri olan bu doktora tezini
canım annem Rüveyde Leyla Öztan’a,
değerli bilim adamları babam Prof. Dr. Şadi Yılmaz Öztan (1930-2003) ile
amcam Prof. Dr. Yüksel Öztan (1933-2010)’ın ölümsüz hatıralarına ve
büyük zorluklarla mücadele ederek eğitimlerini sürdüren yüksek lisans öğrencilerine
ithaf ediyorum
Mehmet Öztan, Mayıs 2011
Tampa, Florida

vi

ACKNOWLEDGEMENTS

This dissertation would have not been possible without the help and support of
professors, funding agencies, colleagues, my family, and friends. First of all, I would like to give
my sincere thanks to Prof. Shu-Guang Li, my dissertation advisor, for his support and advice
throughout my graduate studies at MSU.

I am grateful to my Ph. D. committee members Prof. Amirpouyan Nejadhashemi of
Michigan State University and Dr. Howard W. Reeves of USGS for their detailed comments that
have significantly improved the research and manuscript. Other committee members, Prof. Roger
B. Wallace and Prof. Phanikumar Mantha also gave valuable feedback throughout the research.

I would like to give my sincere thanks to Prof. Ronald Harichandran, the chairperson of
Department of Civil and Environmental Engineering and Prof. Thomas Voice at MSU, for their
continuous support towards my degree completion.

I would like to thank Prof. Huasheng Liao, who is a visiting professor at MSU, and a
professor at Sichuan University (China), for his invaluable comments and help throughout my
research. I consider myself very lucky to know such a special person and admiring scientist.

This research was funded and made possible by the Michigan Department of
Environmental Quality (MDEQ) and the National Science Foundation. Other funders include the
Office for International Students and Scholars at MSU, and the American Water Works
Association. I would like to thank Mr. Richard Mandle from the MDEQ for his help and
guidance regarding the research tasks I completed.

vii

The groundwater modeling research group at MSU is a big family. I am very happy to be
a part of this family and especially thankful to fellow researchers Hassan Abbas, Dr. Dipa Dey,
and Adrian Zhou with who I closely worked in the last four years.

Lori Larner, administrative assistant in the Department of Civil and Environmental
Engineering at MSU, was the person who helped me with every single administrative issue that
alone I wouldn’t be able to resolve.

I am grateful to Prof. M. Yavuz Corapcioglu who encouraged and supported me to come
to MSU for my Ph. D. research.

I am very grateful to Mrs. K. Setenay Arikan, Dr. M. Firat Arikan, and Mr. Ilhan Kubilay
Geckil for their valuable friendship and moral support throughout my graduate studies at MSU.

Finally, I would like to thank my family, specifically my mother Ms. Ruveyde Leyla
Oztan and my cousin Ms. Aysegul Kucukgocmen for their existence and endless support, my
wife Prof. Amy S. Thompson, Dr. Lowery Thompson, Mrs. Barbara Thompson, Ms. Laura
Thompson, and all family friends who have been my encouragement not only in this endeavor,
but all in others.

There is someone else who keeps living deep in my heart, my father Prof. Sadi Yilmaz
Oztan, who I am very thankful for being my role model for hard work and for teaching me that
honesty and human values are above any scientific research.

Mehmet Oztan
May 2011
viii

TABLE OF CONTENTS

LIST OF TABLES.............................................................................................................. x
LIST OF FIGURES .......................................................................................................... xii
LIST OF ABBREVIATIONS........................................................................................ xviii
GIS-BASED REGIONAL-SCALE MODELING OF GROUND WATER FLOW IN
MICHIGAN’S LOWER PENINSULA .............................................................................. 1
INTRODUCTION ................................................................................................... 1
MODEL DEVELOPMENT..................................................................................... 4
Conceptual Model....................................................................................... 4
Governing Equations and Boundary Conditions. ....................................... 8
Numerical Model. ..................................................................................... 10
MODEL EVALUATION ...................................................................................... 14
Statistical analysis..................................................................................... 22
Conceptual messages. ............................................................................... 26
Base flow comparison............................................................................... 32
CONCLUSIONS ................................................................................................... 34
APPENDICES .................................................................................................................. 37
Appendix A: Michigan’s Statewide Groundwater Database ................................. 38
Appendix B: Conceptual Difference between HCOND1 and HCOND2 .............. 68
Appendix C: Data Analysis and Management....................................................... 71
Appendix D: Analysis of the Hydraulic Connection between the Great Lakes and
Regional-Scale Heads .......................................................................................... 101
Appendix E: Analysis of the Simulated Heads.................................................... 107
Appendix F: Analysis of Pumping in the Southern Lower Peninsula ................. 151
Appendix G: Definition of the Additional Model Statistics ................................ 155
Appendix H: Detailed List of the Model Statistics.............................................. 158
Appendix I: MATLAB Script for the Calculation of the Model Statistics.......... 162
LITERATURE CITED ................................................................................................... 168
ix

LIST OF TABLES

Table 1: Model settings and basic information for model analysis. ............................................. 13
Table 2: Explanation of the abbreviations given in Figure 3........................................................ 16
Table 3: Summary of the significant messages illustrated in Figure 4. ........................................ 18
Table 4: Summary of the significant messages illustrated in Figure 5. ........................................ 20
Table 5: Summary of the significant messages illustrated in Figure 6. ........................................ 22
Table 6: Summary of the model statistics for the glacial aquifer simulations.............................. 25
Table 7: Adjusted L ' w values...................................................................................................... 33
Table 8: Hydraulic conductivity classification based on the land systems of Michigan (GWIM,
2006). ............................................................................................................................................ 57
Table 9: Essential modeling parameters recorded in the NHD..................................................... 60
Table 10: Basic specifications of IGW. ........................................................................................ 72
Table 11: Basic features of the main computer used in the research............................................ 73
Table 12: Basic features of the advanced computer used for some of the simulations. ............... 73
Table 13: Criteria used for eliminating the bad data..................................................................... 75
Table 14: Ranges for the color transmissivity map. ..................................................................... 96
Table 15: Prescribed head values assigned for the shorelines of the Great Lakes. .................... 100
Table 16: Temporal analysis of the SWL based on the Great Lakes' water level periods.......... 103
Table 17: Temporal analysis of the SWL based on the annual quarters..................................... 103

x

Table 18: Statistics that are used to quantify the head comparisons........................................... 155
Table 19: Statistics for the model residuals, glacial aquifers...................................................... 158
Table 20: Statistics for the model performances, glacial aquifers. ............................................. 159
Table 21: Statistics for the model residuals, bedrock aquifers. .................................................. 160
Table 22: Statistics for the model performances, bedrock aquifers............................................ 161

xi

LIST OF FIGURES

Figure 1: Illustration of the boundary conditions used in the conceptual model. For interpretation
of the references to color in this and all other figures, the reader is referred to the electronic
version of this dissertation. ............................................................................................................. 7
Figure 2: Conceptual use of the statewide data, and the general cross-section of the conceptual
model............................................................................................................................................... 8
Figure 3: Locations of the modeled watersheds. .......................................................................... 15
Figure 4: Simulated (y axis) and measured (x axis) water levels in the watersheds across the
regional-scale recharge zones. GL: Glacial aquifer, BR: Bedrock aquifer. Scale of the color map
is non-linear. ................................................................................................................................. 17
Figure 5: Simulated (y axis) and measured (x axis) water levels in the watersheds across the
transition zones between regional-scale recharge and discharge zones. GL: Glacial aquifer, BR:
Bedrock aquifer. Scale of the color map is non-linear.................................................................. 19
Figure 6: Simulated (y axis) and measured (x axis) water levels in the watersheds with shale
formations, karst formations, and confining layers. GL: Glacial aquifer, BR: Bedrock aquifer.. 21
Figure 7: a) Distribution of the simulated piezometric heads, Shiawassee Watershed, and b)
Simulated and measured piezometric heads. ................................................................................ 27
Figure 8: a) Simulated piezometric heads in the shale formations, St. Joseph Watershed, and b)
Simulated and measured piezometric heads. ................................................................................ 28
Figure 9: a) Model boundaries before and after the boundary expansion, and b) Comparison
between the simulated and measured water levels before and after the boundary expansion...... 29
Figure 10: Simulated and measured water levels in the Manistee Watershed, glacial aquifer..... 30
Figure 11: Simulated and measured water levels in the bedrock aquifer along the northern tip of
the Lower Peninsula: a) Black Watershed, and b) Cheboygan Watershed................................... 31
Figure 12: Correlation between the simulated and USGS-estimated base flow values................ 33

xii

Figure 13: Comparison between the seepage discharge values when streams are modeled as oneway and two-way head dependent flux model features. ............................................................... 34
Figure 14: Glacial thickness map for Michigan's Lower Peninsula. Values are subject to change
due to the resolution of the raster-based data................................................................................ 40
Figure 15: Major bedrock aquifers in Michigan. .......................................................................... 43
Figure 16: Time-based spatial distribution of the water wells in the glacial aquifers in Michigan.
....................................................................................................................................................... 45
Figure 17: Time-based spatial distribution of the water wells in the bedrock aquifers in
Michigan. ...................................................................................................................................... 46
Figure 18: Time-based change in the representation of the hydraulic conductivity (K) field in the
Kalamazoo County........................................................................................................................ 47
Figure 19: Distribution of the oil and gas wells in Michigan. ...................................................... 49
Figure 20: Raster-based bedrock surface map of Michigan. ........................................................ 50
Figure 21: Watersheds of Michigan.............................................................................................. 52
Figure 22: Sub-watersheds of the Muskegon Watershed. ............................................................ 53
Figure 23: Conceptualization of the calculation interval for the HCOND1 parameter. ............... 54
Figure 24: Conceptualization of the calculation interval for the HCOND2 parameter. ............... 55
Figure 25: Stream segments with missing stream stage data in the NHD of Michigan. .............. 61
Figure 26: Raster-based recharge map of Michigan. .................................................................... 64
Figure 27: Locations of the stream gauging stations in Michigan................................................ 67
Figure 28: Hydraulic conductivity field of the Tittabawassee Watershed. Above HCOND1, and
Below: HCOND2.......................................................................................................................... 69
Figure 29: Comparison between the model results with HCOND1 and HCOND2. Black Dots:
Results for HCOND1, and Blue Circles: Results for HCOND2. ................................................. 70

xiii

Figure 30: Error map for the DEM information in Wellogic........................................................ 77
Figure 31: Location of an area in the western Lower Peninsula where DEM error is extremely
high. .............................................................................................................................................. 79
Figure 32: Simulated and measured heads for the model area shown in Figure 31. Above:
comparison using vector-based DEM, and Below: comparison using raster-based DEM with 10
m (sampled from 30 m) resolution................................................................................................ 80
Figure 33: Identifying bad SWL data using model-based analysis, Thornapple Watershed,
bedrock aquifer. ............................................................................................................................ 81
Figure 34: Identifying the outliers for the Clinton Watershed, bedrock aquifer. ......................... 83
Figure 35: Analysis of the outliers based on the model results in the Clinton Watershed, bedrock
aquifer. .......................................................................................................................................... 84
Figure 36: Error locations for the bedrock surface map of the Upper Peninsula. ........................ 86
Figure 37: Magnifying the error in the bedrock surface map and available data points in the
Upper Peninsula. ........................................................................................................................... 87
Figure 38: Conceptual picture of the possible reason for the error in the bedrock surface map. . 88
Figure 39: USGS 90 m raster-based DEM. .................................................................................. 91
Figure 40: Raster-based hydraulic conductivity maps interpolated from the data. ...................... 93
Figure 41: Color map for the transmissivity ranges generated based on the aquifer test data. .... 95
Figure 42: High order streams with artificial lines along the shores of the Great Lakes. ............ 97
Figure 43: Conceptualization of the stream cells in the models. .................................................. 99
Figure 44: Location of the test area for temporal SWL analysis. ............................................... 102
Figure 45: Temporal SWL analysis for annual quarters in the Kalamazoo and Van Buren County.
..................................................................................................................................................... 104
Figure 46: Temporal analysis for the Great Lakes water level periods in the Kalamazoo and Van
Buren County. ............................................................................................................................. 105
xiv

Figure 47: Analysis of the simulated heads in the Au Gres-Rifle Watershed, glacial aquifer. .. 107
Figure 48: Analysis of the simulated heads in the Au Gres – Rifle Watershed, bedrock aquifer.
..................................................................................................................................................... 108
Figure 49: Analysis of the simulated heads in the Au Sable Watershed, glacial aquifer. .......... 109
Figure 50: Analysis of the simulated heads in the Birch/St. Clair watersheds, glacial aquifer. . 110
Figure 51: Analysis of the simulated heads in the Birch/St. Clair watersheds, bedrock aquifer.111
Figure 52: Analysis of the simulated heads in the Black Watershed, glacial aquifer................. 112
Figure 53: Analysis of the simulated heads in the Black Watershed, bedrock aquifer. ............. 113
Figure 54: Analysis of the simulated heads in the Cass Watershed, glacial aquifer. ................. 114
Figure 55: Analysis of the simulated heads in the Cass Watershed, bedrock aquifer. ............... 115
Figure 56: Analysis of the simulated heads in the Cheboygan Watershed, glacial aquifer........ 116
Figure 57: Analysis of the simulated heads in the Cheboygan Watershed, bedrock aquifer...... 117
Figure 58: Analysis of the simulated heads in the Clinton Watershed, glacial aquifer. ............. 118
Figure 59: Analysis of the simulated heads in the Clinton Watershed, bedrock aquifer............ 119
Figure 60: Analysis of the simulated heads in the Flint Watershed, glacial aquifer. ................. 120
Figure 61: Analysis of the simulated heads in the Flint Watershed, bedrock aquifer. ............... 121
Figure 62: Analysis of the simulated heads in the Huron Watershed, glacial aquifer................ 122
Figure 63: Analysis of the simulated heads in the Huron Watershed, bedrock aquifer.............. 123
Figure 64: Analysis of the simulated heads in the Kalamazoo Watershed, glacial aquifer........ 124
Figure 65: Analysis of the simulated heads in the Kalamazoo Watershed, bedrock aquifer...... 125
Figure 66: Analysis of the simulated heads in the Kawkawlin Watershed, glacial aquifer........ 126
xv

Figure 67: Analysis of the simulated heads in the Kawkawlin Watershed, bedrock aquifer. .... 127
Figure 68: Analysis of the simulated heads in the Lower Grand River Watershed, glacial aquifer.
..................................................................................................................................................... 128
Figure 69: Analysis of the simulated heads in the Lower Grand River Watershed, bedrock
aquifer. ........................................................................................................................................ 129
Figure 70: Analysis of the simulated heads in the Manistee Watershed, glacial aquifer. .......... 130
Figure 71: Analysis of the simulated heads in the Maple Watershed, glacial aquifer................ 131
Figure 72: Analysis of the simulated heads in the Maple Watershed, bedrock aquifer.............. 132
Figure 73: Analysis of the simulated heads in the Muskegon Watershed, glacial aquifer. ........ 133
Figure 74: Analysis of the simulated heads in the Pere-Marquette Watershed, glacial aquifer. 134
Figure 75: Analysis of the simulated heads in the Pigeon Watershed, bedrock aquifer............. 135
Figure 76: Analysis of the simulated heads in the Pine Watershed, glacial aquifer................... 136
Figure 77: Analysis of the simulated heads in the Raisin Watershed, glacial aquifer................ 137
Figure 78: Analysis of the simulated heads in the Raisin Watershed, bedrock aquifer.............. 138
Figure 79: Analysis of the simulated heads in the Shiawassee Watershed, glacial aquifer........ 139
Figure 80: Analysis of the simulated heads in the Shiawassee Watershed, bedrock aquifer. .... 140
Figure 81: Analysis of the simulated heads in the St. Joseph Watershed, glacial aquifer.......... 141
Figure 82: Analysis of the simulated heads in the St. Joseph Watershed, bedrock aquifer........ 142
Figure 83: Analysis of the simulated heads in the Thornapple Watershed, glacial aquifer........ 143
Figure 84: Analysis of the simulated heads in the Thornapple Watershed, bedrock aquifer. .... 144
Figure 85: Analysis of the simulated heads in the Thunder Bay Watershed, glacial aquifer. .... 145
Figure 86: Analysis of the simulated heads in the Thunder Bay Watershed, bedrock aquifer... 146
xvi

Figure 87: Analysis of the simulated heads in the Tittabawassee Watershed, glacial aquifer. .. 147
Figure 88: Analysis of the simulated heads in the Tittabawassee Watershed, bedrock aquifer. 148
Figure 89: Analysis of the simulated heads in the Upper Grand River Watershed, glacial aquifer.
..................................................................................................................................................... 149
Figure 90: Analysis of the simulated heads in the Upper Grand River Watershed, bedrock
aquifer. ........................................................................................................................................ 150
Figure 91: Interpolated SWL in the city of Lansing. .................................................................. 152
Figure 92: Revisiting the model results for the Upper Grand River. Circled area indicate the
pumping effect. ........................................................................................................................... 153
Figure 93: Identifying biases similar to that of the Upper Grand River case. Circled areas might
potentially indicate the pumping effect. Above: Flint Watershed, Below: Shiawassee Watershed.
..................................................................................................................................................... 154

xvii

LIST OF ABBREVIATIONS

1/d: One per day
ac: acre
DEM: Digital Elevation Model
ft: feet
ft/d: feet per day
2

ft /d: square feet per day
3

ft /s: cubic feet per second
GIS: Geographic Information System
GWIM: Groundwater Inventory and Mapping Project
HUC: Hydrological Unit Code
in/yr: inch per year
K: Hydraulic conductivity
LERCMM: Laboratory of Excellence for Real-time Computing and Multi-scale Modeling
m: meter
m/d: meter per day
m/s: meter per second
3

m /s: cubic meter per second
2

mi : square mile
MDEQ: Michigan Department of Environmental Quality

xviii

MDNRE: Michigan Department of Natural Resources and Environment
MSGWD: Michigan’s Statewide Groundwater Database
MSU: Michigan State University
NED: National Elevation Dataset
NHD: National Hydrography Dataset
NOAA: National Oceanic and Atmospheric Administration
RASA: Regional Aquifer Systems Analysis
SRTM: Shuttle Radar Topography Mission
SWL: Static water level
T: Transmissivity
USGS: United States Geological Survey

xix

GIS-BASED REGIONAL-SCALE MODELING OF GROUND WATER FLOW IN
MICHIGAN’S LOWER PENINSULA

INTRODUCTION

In the last two decades, advances in information and computational technologies and
increasing demands in ground water use have shifted the focus of ground water investigations
from stand-alone models towards GIS-based modeling (Barazzuoli et al., 1999; Batelaan et al.,
1993; El-Kadi et al., 1994; Gogu et al., 2001; Henriksen et al., 2003; Hoaglund et al., 2002;
Hojberg et al., 2007; Li et al., 2009; Liao et al., 2010; Lubczynski and Gurwin, 2005; Mende et
al., 2007; Michael et al., 2005; Refsgaard et al., 2009; Tim et al., 1996). Recently, Li et al. at
Michigan State University (MSU) developed a data-integrated, interactive and graphically
oriented modeling system, Interactive Ground water (IGW), to investigate Michigan’s ground
water systems in collaboration with Michigan Department of Natural Resources and
Environment (MDNRE) (Li and Liu, 2006). IGW enables real-time two-dimensional and threedimensional modeling, visualization, and analysis for many ground water problems (Li and Liu,
2006; Zheng, 2010). Most importantly, IGW provides a coupling between process-based models
and Michigan’s existing statewide ground water database (MSGWD) in a way that allows the
definition of the model parameters almost everywhere in the state and at the time of numerical
model creation, a unique feature that makes the regional-scale statewide ground water modeling
possible. In this paper, we take full benefit of this feature to evaluate the modeling system’s
capability to simulate regional-scale ground water flow in Michigan’s Lower Peninsula. This
evaluation helps us to improve the understanding of the regional-scale ground water flow in the
glacial aquifer and its interaction with the bedrock aquifer across the peninsula.
1

MSGWD is a pre-processed GIS-based database and has been assembled by the
Michigan Department of Environmental Quality (now MDNRE) to contain most of the data
needed for regional-scale ground water flow simulations across the state (GWIM, 2006). The
database has been created based on the requirement made by the Public Act 148, Section 32802
(2003) that aims to improve the statewide environmental and water resources management. All
data are available in geo-referenced digital files throughout our application domain, and spatial
variations of the aquifer parameters in the models are defined with these data. Specifically,
MSGWD includes the data for defining numerical model domains (e.g. watershed boundaries),
aquifer elevations (e.g. digital elevation models (DEMs), aquifer lithologies), aquifer stresses
(e.g. estimated ground water recharge, streams), aquifer properties (e.g. hydraulic conductivities,
transmissivities), and calibration data (e.g. static water levels (SWL), base flow).

We focus on the ground water flow in two aquifer layers in the Lower Peninsula: the
glacial and uppermost bedrock aquifer. The most critical conceptual questions we ask are the
following: Can we simulate the dominant variations in regional-scale ground water flow in these
two aquifers and explain the variations in measured ground water levels by assigning only one
computational layer for each aquifer in the numerical model? What is the controlling factor of
the bedrock flow along the regional-scale discharge zones and in the watersheds with shale units,
relatively small glacial thickness and local confining layers? Can we conceptualize the bedrock
aquifer as a collection of geologic sub-crops that are in contact with the glacial aquifer and that
represent the effective transmissivity field? What is the impact of modeling glacial aquifer with
one computational layer in the regional-scale recharge zones? Can we use surface-watershed
boundaries to represent the regional-scale no-flow ground water boundaries everywhere in the
peninsula? For which watersheds is the modeling system insufficient to simulate the ground
2

water flow, and why? We seek the answers to these conceptual questions by simulating ground
water flow for 26 of 38 USGS 8-digit Hydrologic Unit Codes (hereafter referred as watersheds)
across the peninsula. We tested the modeling system on a large number of watersheds in order to
ensure that the capability of this system is generalizable. Extending the analysis over such a large
area enables the coverage of the aquifers that are apart and that -at a regional-scale- have
different hydrogeologic and physical properties (e.g. different recharge/discharge patterns,
different aquifer thicknesses, different lithologies, different watershed sizes, and variations in the
magnitude of vertical flow). Furthermore, maximizing the number of the modeled watersheds
provides the opportunity to systematically test the capability of the simple two-layer conceptual
models to represent the regional-scale ground water flow in the glacial and bedrock aquifers in
the peninsula. Additionally, we assume that the watershed boundaries coincide with the regionalscale ground water basin boundaries to evaluate the accuracy of the boundary conditions used in
the models.

We use IGW v. 5.20 to develop all of the models. Along with IGW, Arc/Info and
MATLAB software applications are used for the model evaluation. Goodness of fit of the
simulated long-term static water levels in the glacial aquifers and long-term piezometric heads in
the bedrock aquifers to measured water level data from MSGWD is evaluated with model
statistics. Similarly, simulated base flow values are systematically compared to USGS-estimated
base flow values to evaluate the goodness of model calculations. Additionally, the comparisons
between the simulated and measured water levels in each watershed are used to point out the
most critical hydrogeologic aspects of the regional-scale ground water flow in the peninsula, and
to evaluate the assumptions used in the conceptual models.

3

MODEL DEVELOPMENT

Conceptual Model. The conceptual model focuses on two regional aquifer layers: the
glacial aquifer and uppermost bedrock aquifer. The watershed boundaries (USGS, 2011) are the
model boundaries, and are initially assumed to coincide with the regional-scale ground water
basin boundaries (Figure 1). Main purpose of this assumption is to test whether the watershed
boundaries can be used for regional-scale ground water modeling across the Lower Peninsula of
Michigan. Another assumption is that there is no lateral ground water flux coming to or from the
glacial and bedrock aquifer through the watershed boundaries. Lateral flow only occurs into the
Great Lakes through the watershed boundaries along the lake shorelines. The base of the bedrock
aquifer is assumed to be impermeable.

The upper elevation of the glacial aquifer (i.e. land surface elevation) is defined by the
USGS 90 m National Elevation Dataset (NED) (USGS, 2006). Where applicable, land surface is
conceptualized as the seepage surface. That is, discharge to the land surface is zero when water
levels in the model cells adjacent to the land surface are equal to or less than the assigned land
surface elevation (Anderson and Woessner, 2002). All polygon-based surface water bodies (e.g.
inland lakes, ponds, wetlands) are also assumed to be the components of the seepage surface.
Steady-state water levels in these water bodies are assigned from the USGS 90 m NED. The
bottom elevation of the glacial aquifer (bedrock surface) across the state is defined based on
almost 150,000 well logs of water, and oil and gas wells. The locations and information in the
water wells are provided with the state’s electronic water well database, Wellogic (GWIM,
2006). The database for the oil and gas wells is compiled by the Michigan Department of Natural
Resources and Environment (GWIM, 2006). The statewide effective horizontal hydraulic
4

conductivity field of the glacial aquifer is defined on the basis of lithologic descriptions recorded
in more than 200,000 water well records from the statewide Wellogic database (GWIM, 2006).
The vertical hydraulic conductivity of the glacial aquifer is defined based on the vertical
anisotropy given with

Kx K y

 K z where K x is the effective hydraulic conductivity in x
10
10

direction, K y is the effective hydraulic conductivity in y direction, and K z is the effective
hydraulic conductivity in z direction. The statewide effective transmissivity (Tx for the effective
transmissivity in x direction and Ty for the effective transmissivity in y direction) field of the
bedrock aquifer is pre-defined based on the spatial interpolation of more than 1,400 aquifer test
data provided by the USGS. The bedrock aquifer thickness is assumed to be constant and
approximately 76 m. Since the Tx and Ty terms are pre-defined, changing the thickness only has
an impact on the vertical flux [Equation (4)] between the glacial and bedrock aquifer. However,
based on the available well depths that are provided with the USGS aquifer test database, a
bedrock aquifer thickness of 76 m was found to be reasonable, even though the depths of the
wells in the database are not necessarily the same.

It is assumed that the source of ground water to the glacial aquifer is from precipitation
recharge. The statewide long-term average ground water recharge to the glacial aquifer is
estimated by the USGS with hydrograph separation for more than 400 stream flow hydrographs,
and landscape attributes (Neff et al., 2005). The only source of water to the bedrock aquifer is
maintained through the vertical leakage from the overlying glacial aquifer. The locations of the
streams are defined with the USGS National Hydrography Dataset (NHD) for all 1st to 6th order
streams (USGS, 2011). Steady-state water levels in the stream segments are assigned to all

5

 

streams from the USGS 30 m NED (USGS, 2006). The leakance L' parameter is used to define
the hydraulic connection between the seepage surface/streams and ground water [Equation (1)].

L' 

K'
b'

1

where K ' is the hydraulic conductivity of the stream/seepage surface bed and b' is the thickness
of stream/seepage surface bed. All streams are assumed to be effluent streams since the stream
gauging stations and miscellaneous measurements in Michigan show that most streams are
gaining streams in the state (Holtschlag and Nicholas, 1998). Where applicable, the water
elevations along the shorelines of the Great Lakes are assumed to be equal to the long-term
average of the observed lake levels reported by the National Oceanic and Atmospheric
Administration (NOAA) (NOAA, 2011). These levels are represented in the watershed models
for both the glacial and bedrock aquifers as prescribed-head boundaries (Figure 1).

Conceptual use of the data and a general cross-section of the conceptual model are
illustrated in Figure 2.

6

BC2
BC1

State boundary
BC1 Watershed boundary/No-flow ground water boundary
BC2 Great Lakes shoreline/prescribed head boundary, where applicable
Great Lakes
Figure 1: Illustration of the boundary conditions used in the conceptual model. For interpretation
of the references to color in this and all other figures, the reader is referred to the electronic
version of this dissertation.

7

SWL, K, aquifer elevations
(Wellogic)
Stream stage
(USGS NHD)
Great Lakes’ levels
(NOAA)

Bedrock surface
(oil/gas/water well)

Recharge
(USGS)

Lake/Wetland level
(USGS NED)

Glacial
Base flow aquifer
(USGS)

Land surface
(USGS NED)

Seepage
LH

LM

Saginaw Aquifer
Coldwater
Shale

Coldwater
Shale
Other
formations

Other
Marshall formations
Aquifer

Marshall
Aquifer

Collection of the sub-crops of the bedrock lithologies/Conceptual and computational
bedrock layer
Impermeable bedrock base
LM Lake Michigan

LH Lake Huron

Figure 2: Conceptual use of the statewide data, and the general cross-section of the conceptual
model.

Governing Equations and Boundary Conditions. Ground water flow in the glacial and
bedrock aquifer is represented with [Equation (2)] and [Equation (3)], respectively:

hg   
hg 
 
K x hg  e
  K y hg  e


  qg  qv  0
x 
x  y 
y 









 2

where hg is the ground water level in the glacial aquifer, e is the bottom elevation of the glacial
aquifer, q g is the source (i.e. recharge) and sinks (i.e. seepage, pumping where applicable), and

8

qv is the vertical flux between the glacial and bedrock aquifer and expressed with [Equation
(4)].

  hb    hb 




 Tx x   y  Ty y   qb  qv  0
x 




 3

where hb is the ground water level in the bedrock aquifer, qb is the sinks (i.e. pumping where
applicable).

'
Kz
qv 
hg  hb
d





 4

'
where K z is the conductance between the glacial and bedrock aquifer, d is the arithmetic average

of the aquifer thicknesses, and

'
Kz 

d



K zg K zb bg  bb



 5

bg K zb  bb K zg

bg  bb

 6

2

where K zg is the vertical hydraulic conductivity of the glacial aquifer, K zb is the vertical
hydraulic conductivity of the bedrock aquifer, bg is the glacial aquifer thickness, and bb is the
bedrock aquifer thickness.

9

Based on the conceptual model, two different boundary conditions given in [Equation
(7)] and [Equation (8)] are used to solve [Equation (2)] and [Equation (3)]:

hg
n 
1

hg

2



hb

n 
1

 hb

2

7

0

 H  x, y 

 8

where 1 is watershed boundary,  2 is the Great Lakes’ shoreline, and H  x, y  is the Great
Lakes’ long-term average water level.

Numerical Model.

Quasi three-dimensional finite-difference numerical models are

used to simulate the ground water flow in 26 of the 38 watersheds across Michigan’s Lower
Peninsula. The models are represented by rectangular grids discretized into rows and columns
which construct the model cells. The values for the model input parameters in the model cells are
assigned from the raster and vector-based GIS files (MSGWD) that are coupled with the
numerical models to allow the representation of the aquifer heterogeneities. In the models, each
aquifer layer (i.e. glacial and bedrock aquifer layer) is represented by one computational layer
that simulates one average water level for the entire thickness of the aquifer.

At a regional scale, configuration of the water table is commonly a subdued replica of the
land-surface topography (USGS, 1999). One computational layer is used to simulate the glacial
aquifer based on the general practical rule that dictates that the Dupuit-Forchheimer
approximation is acceptable when the distance between hydrologic boundaries, L, (here tens of

10

kilometers in most watersheds) is large compared to the saturated aquifer thickness, H,
(Haitjema,2005) [Equation (9)]:

L

Kh
Kv

9 

H

where Kh is the horizontal hydraulic conductivity, and Kv is the vertical hydraulic conductivity.

All models represent the steady-state flow conditions for which the simulated ground
water levels are the long-term average values and water levels don’t change with time. However,
measured water levels in the wells span many years, in some cases over fifty years, and models
do not incorporate the temporal variations in measured water levels due to the impacts of
seasonal fluctuations and ground water withdrawals.

Cell sizes used in the models are not fine enough to capture the head variations induced
by small-scale spatial heterogeneities (e.g. perched aquifers, local confining layers). Even though
these small-scale details are not the focus of this paper, it should be noted that they contribute to
the error in the simulations. Additionally, streams become strong boundary conditions as the
stream order and the hydraulic connection (i.e. leakance) between the stream and ground water
increases, and the ground water levels simulated in the models converge to the stream stages in
these stream cells. Cell sizes are not fine enough to capture the small-scale variations in stream
levels along the stream courses.

Most of the wells recorded in Wellogic do not have a pumping rate. More than 90% of
the wells are household wells, and if it is assumed that each well in the watersheds is active and
11

-6

pump ground water with an average discharge rate of 4x10 m3/d, then in none of the
watersheds, magnitude of the ratio between total discharge through the wells and total ground
-2

water recharge is more than 10 . Furthermore, at a regional-scale, it is not possible to capture
the actual drawdown created with small-scale pumping. Because of these reasons, pumping wells
are not simulated in the models.

The information regarding the number of the cells, and the grid sizes used in the models,
the model area and the number of the water-well records (i.e. number of measured water level
data) used for the analysis of the simulations are given in Table 1.

12

Table 1: Model settings and basic information for model analysis.
Watershed
Name

NX

NY

∆X
(m)

∆Y
(m)

Model
Area
(km )

No. of
Glacial
Wells

No. of
Bedrock
Wells

2

Au Gres – Rifle

130

87

654

653

2,593

1,834

562

Au Sable

150

103

837

835

5,323

3,184

7

Birch – St. Clair

150

363

437

437

4,552

2,770

1,234

Black

50

91

768

768

1,560

369

158

Cass

150

108

636

635

2,335

109

2,159

Cheboygan

71

103

768

768

2,429

1,446

77

Clinton

150

138

404

405

1,903

12,057

181

Flint

150

116

569

566

3,499

3,611

9,750

Huron

150

139

574

573

2,432

19,938

2,358

Kalamazoo

150

97

982

978

5,220

12,730

2,884

Kawkawlin

150

245

238

238

1,294

623

1,735

L. Grand River

150

124

732

732

5,185

14,162

1,298

Manistee

150

137

845

847

4,925

5,222

2

Maple

150

106

571

572

2,413

1,173

1,106

Muskegon

150

153

1,022

1,025

7,066

15,054

17

Pere – Marquette

150

212

521

522

5,192

10,227

6

Pigeon

150

147

503

502

2,239

22

1,308

Pine

150

121

581

580

2,677

4,305

15

Raisin

150

96

636

633

2,641

4,403

1,733

Shiawassee

150

185

630

629

3,198

10,307

2,829

St. Joseph

250

101

676

673

6,174

15,448

972

Thornapple

150

109

525

525

2,203

5,034

1,508

Thunder Bay

100

91

768

768

3,257

1,035

230

Tittabawassee

150

186

544

543

3,685

5,744

223

U. Grand River

150

173

612

610

4,500

2,032

17,579

13

MODEL EVALUATION

At a regional-scale, resolution of the water-well data (i.e. Wellogic) is several magnitudes
greater than that of the traditionally collected data, which is a unique feature that makes the GISenabled regional-scale ground water modeling possible almost everywhere in the State of
Michigan. The GIS interface of IGW maximally benefits from the data-rich structure of
Wellogic. The measured water levels are readily available in Wellogic, and can be brought into
IGW through the GIS interface in real-time and interactively. The coupling between the data and
models is used to evaluate the simulations by systematically comparing the simulated water
levels to the measured water levels (Figure 4, Figure 5, and Figure 6). IGW calculates the
simulated water levels at the nodes of each cell, and then using the information at the nodes,
applies bi-linear interpolation to predict the water levels at the location of the water wells.
Measured water levels that are used for model comparison are the water levels recorded in the
water wells after the wells are developed, prior to pumping (GWIM, 2006). Please refer to Figure
3 and Table 2 for the locations and identifications of the simulated watersheds. Summary of the
significant messages illustrated in Figure 4, Figure 5, and Figure 6 is provided in Table 3, Table
4, and Table 5.

14

Figure 3: Locations of the modeled watersheds.

15

Table 2: Explanation of the abbreviations given in Figure 3.
Watershed Name
Abbreviation
Au Gres – Rifle

GR

Au Sable

SB

Birch – St. Clair

BI

Black

BL

Cass

CA

Cheboygan

CH

Clinton

CL

Flint

FL

Huron

HU

Kalamazoo

KL

Kawkawlin

KA

Lower Grand River

LG

Manistee

MN

Maple

MA

Muskegon

MU

Pere-Marquette

MR

Pigeon

PI

Pine

PE

Raisin

RA

Shiawassee

SH

St. Joseph

JO

Thornapple

TO

Thunder Bay

TU

Tittabawassee

TI

Upper Grand River

UG

16

400

400

GL

400

GL

400

GL

300

200

300

300

200

300

200

200

200 300 400
Manistee (MN)
400

400

300

300

150

300

GL
250

250

400

GL

GL

300
200

200
175 250 325
Pine (PE)

200 275 350
Thunder Bay (TU)

200 300 400
Tittabawassee (TI)

250

175

300

200

150 250 350
Marquette (MR)

325

GL

275

200

250

200

200 300 400
Black (BL)
350

GL

GL

200 300 400
Muskegon (MU)

200 300 400
Cheboygan (CH)

200 300 400
Au Sable (SB)
350

GL

GL

200 250 300
Huron (HU)

200 300 400
Au Gres (GR)
300
GL

GL
225

200
200
300

150

250 300

BR

150

PI

280

250

300

BR

230

200
200 250 300
U. Grand R. (UG)

225

180
DEM (m)
160
505
Uplands (regional-scale recharge zones)

180 230 280
Thornapple (TO)

Figure 4: Simulated (y axis) and measured (x axis) water levels in the watersheds across the
regional-scale recharge zones. GL: Glacial aquifer, BR: Bedrock aquifer. Scale of the color map
is non-linear.
17

Table 3: Summary of the significant messages illustrated in Figure 4.
Northern upland (NU): Largest glacial thickness, significant vertical flow, deep bedrock
aquifers (i.e. very few wells for model validation)
- On NU: Simulated glacial heads are averaged over large aquifer depths and systematically
under predict measured heads (e.g. MN, MR, MU, SB)
- Around NU: Simple no-flow boundaries cannot be used on watershed boundaries to model
glacial flow (e.g. BL, CH, GR, MR, TI, TU)
Southern upland (SU): Glacial thickness is moderate with respect to NU, significant vertical
flow
- On SU: Systematic model under prediction still exists but not as strong as in NU (e.g. HU
(GL), UG (GL), PE (GL), TO (GL))
- Around SU: Simple no-flow boundaries cannot be used on watershed boundaries for
ground water modeling, bias at high heads are less due to smaller gradient/lateral flow at the
boundaries (e.g. HU (GL), TO (BR), UG (GL))
- In major pumping centers, simulated heads locally over predict the measured heads since
pumping wells are not modeled (e.g. circled area in UG (BR))

18

280

Glacial Thickness (m)
1
377

GL
230

325

GL

250
175

180
180
250

230 280

175

250 325

300

BR

225

250

200

BR

200
200 225 250
Maple (MA)

200 250 350
Shiawassee (SH)

350

260

GL

325

250

GL
220

GL

250

150
150

250

350

180

175

350

180 220 260
260
175

BR
275

325

260

325

200

250

250

220

BR

200 275 350
Kalamazoo (KL)

350
GL

180
180 220 260
Cass (CA)

175
175 250 325
Flint (FL)
300

BR

GL

325

GL

300

250

225

250

250

150

150

175

GL

200

150 250 350
St. Joseph (JO)

175 250 325
Clinton (CL)

150 225 300
L. Grand R. (LG)

200 250 300
Raisin (RA)

Figure 5: Simulated (y axis) and measured (x axis) water levels in the watersheds across the
transition zones between regional-scale recharge and discharge zones. GL: Glacial aquifer, BR:
Bedrock aquifer. Scale of the color map is non-linear.
19

Table 4: Summary of the significant messages illustrated in Figure 5.
Transition zones:
- Glacial thickness gradually decreases from uplands towards discharge zones along
lowlands (regional-scale drainage along the central peninsula), and the Great Lakes shores
- At low elevations, major streams function as strong boundary conditions and control glacial
flow (e.g. CL, FL, KL, LG, SH)
- Through vertical hydraulic connection, major streams also control bedrock flow at lower
elevations and in regions where bedrock is very close to land surface (e.g. CA, FL, KL, SH)
- Even in watersheds where glacial aquifer is intercalated with confining layers at low
elevations, major streams control bedrock flow (e.g. FL, CA, SH)
- Underneath recharge mounds close to regional-scale discharge zones, simulated glacial
heads are averaged over moderate glacial thickness (e.g. CL, KL, LG, MA, RA)
- At high elevations, simulated heads under predict due to boundary conditions (e.g. CA (GL
and BR), FL (GL and BR), JO (GL), KL (GL and BR), LG (GL), MA (GL and BR), SH
(GL))

20

280

260

250

350

BR
230

220

200

300

150

250

BR

BR

180

180
180 220 260
Cheboygan (CH)

150 200 250
L. Grand R. (LG)

180 230 280
Black (BL)
280

240

St. Joseph (JO)
325

325

BR

BR

BR

230
BR

180 210 240
Thunder Bay (TU)

250

225

180

210
180

175

125

180 230 280
Au Gres (GR)

175 250 325
Raisin (RA)

260

300

BR

BR

BR

250

250

BR

220

210

200

190

160
175 250 325
Huron (HU)

125 225 325
Clinton (CL)

160 210 260
Pigeon (PI)

250

190 220 250
Tittabawassee (TI)

300

GL

GL

250

210

200

170
170

210

250
300

200 250 300
BR

220
250
190
160

200

BR

160 190 220
Kawkawlin (KA)

200 250 300
Confining layer
Birch-St. Clair (BI)
Shale
Limestone and dolomite
Figure 6: Simulated (y axis) and measured (x axis) water levels in the watersheds with shale
formations, karst formations, and confining layers. GL: Glacial aquifer, BR: Bedrock aquifer.
21

Table 5: Summary of the significant messages illustrated in Figure 6.
Shale formations:
- Glacial thickness above shale units gradually decreases from southern uplands towards the
Great Lakes shorelines
- Hydraulic connection between glacial and bedrock aquifer is maintained through vertical
flow
- Major streams in glacial aquifers connect with bedrock aquifers through cracks in shale
formations and control bedrock flow (e.g. CL, HU, JO, RA)
Confining layers along the shores:
- Glacial thickness significantly decreases in the watersheds close to the Great Lakes
- At high elevations, ground water in bedrock aquifers regionally discharges into the Great
Lakes through deep bedrock flow system underneath the confining layers along the shores
- Artesian wells exist in the watersheds next to the Great Lakes (e.g. GR, KA)
- Systematic over prediction in glacial aquifers above the confining layers is partially related
to pumping near the shores since wells are not simulated (e.g. KA, LG)
- Systematic over prediction in bedrock aquifers underneath the confining layers may
potentially be related to the under estimation of transmissivity values due to low data
resolution (e.g. BI, GR, KA, PI, RA, TI)
Karstic formations: Significant bias in the simulated heads, Darcy’s Law doesn’t apply to the
karstic bedrock aquifers along the northern tip of the peninsula (e.g. CH, BL, TU)

Statistical analysis. The goodness-of-fit of the simulated water levels (Figure 4, Figure
5, and Figure 6) to the measured data is evaluated with the model statistics. These statistics are
used to demonstrate the overall performance of each model through evaluating the average error
in the models (Anderson and Woessner, 2002). Following statistics are used for the evaluation:

1) The arithmetic mean error (AME). Even though the positive and negative residuals (i.e.
difference between the simulated and measured ground water levels) may cancel each other out,
the AME allows quantifying the average model bias through a simple calculation (Anderson and
Woessner, 1991; El-Kadi (Ed.), 1995).

22

n

AME 

  hmi  hsi 
i 1

10 

n

where n is the total number of measurements in a watershed, hmi is the measured ground water
level, and hsi is the simulated ground water level.

2) The normalized standard deviation of the residuals (NSDR). Generally, if the range of the
measured water-level data is large, there is a larger standard deviation in residual errors. In
ground water investigations, a good fit to the data would be reflected if this ratio is
approximately 0.1 or less indicating that the residuals are generally less than 10 percent of the
range in altitude of the observations (Kuniansky et al. 2004; Hunt et al. 2003).

n

  ri   
i 1
NSDR 

2

n 1
DR

11

where ri is the particular absolute residual, µ is the residual mean, and DR is the data range
(Walton, 1992).
2

2

3) The coefficient of determination (R ). R indicates the strength of the fit between the
2

simulated and measured ground water levels through the dispersion in the comparison. R can be
expressed as:

23

n

  hmi  hm  hsi  hs 
i 1

2

R 

n

n

  hmi  hm    hsi  hs 
i 1
i 1
2

12

2

where hm is the average of the measured ground water levels, and hs is the average of the
2

simulated ground water levels. The range of R changes between 0 and 1where a value of 0
indicates no correlation and, ideally, a value of 1 means that the dispersion of the prediction is
equal to that of the observation (Krause et al. 2005).

A summary of the model statistics is provided in Table 2. The AME, and NSDR are
2

generally within 10% of the total head difference (DR) within the watersheds, and R values
indicate that for most watersheds more than 85% of the total variation in the graphical
comparisons (Figure 4, Figure 5, and Figure 6) can be explained with the linear relation between
the simulated and measured water levels.

24

Table 6: Summary of the model statistics for the glacial aquifer simulations.
Watershed Name

Glacial Aquifer

Bedrock Aquifer
2

2

DR
(m)

AME
(m)

NSDR
(%)

R

DR
(m)

AME
(m)

NSDR
(%)

R

Au Gres – Rifle

232

5.7

4.8

0.89

105

-0.1

4.4

0.96

Au Sable

268

6.1

4.4

0.96

-

-

-

-

Birch – St. Clair

139

-0.7

5.0

0.92

125

-1.5

4.8

0.94

Black

228

6.9

7.4

0.91

100

-13.2

17.4

0.26

Cass

83

3.9

9.5

0.78

102

-0.1

5.7

0.87

Cheboygan

253

16.6

8.5

0.96

105

-1.9

14.8

0.67

Clinton

185

3.9

3.8

0.98

212

-0.9

6.2

0.88

Flint

170

7.1

6.0

0.95

166

-2.4

3.7

0.93

Huron

139

3.8

4.4

0.83

140

1.2

3.6

0.95

Kalamazoo

205

3.4

3.9

0.90

171

3.3

4.1

0.95

Kawkawlin

77

-3.8

8.5

0.80

63

-4.7

6.8

0.80

L. Grand River

164

8.2

6.0

0.86

102

-3.2

8.6

0.73

Manistee

294

20.4

6.5

0.91

-

-

-

-

Maple

102

3.3

6.0

0.89

64

-0.2

7.9

0.84

Muskegon

284

11.5

4.4

0.96

-

-

-

-

Pere – Marquette

230

12.0

5.4

0.89

-

-

-

-

Pigeon

-

-

-

-

90

-5.2

5.2

0.91

Pine

168

5.3

4.3

0.97

-

-

-

-

Raisin

156

2.4

4.1

0.94

172

-5.7

5.0

0.99

Shiawassee

174

3.2

3.7

0.97

144

-1.5

3.4

0.97

St. Joseph

205

2.7

3.1

0.95

141

2.7

5.1

0.93

Thornapple

140

8.0

5.8

0.86

107

2.1

5.2

0.90

Thunder Bay

206

12.7

8.3

0.82

75

-3.1

13.4

0.71

Tittabawassee

217

17.6

10.1

0.89

74

-1.9

8.7

0.89

U. Grand River

139

2.8

4.0

0.91

133

-0.1

4.4

0.89

25

Conceptual messages.

When combined with the model statistics listed in Table 6,

evaluation of the comparisons shown in Figure 4, Figure 5, and Figure 6 provides the following
conceptual messages regarding the ground water systems in the Lower Peninsula of Michigan:

-

The conceptual model assumes that the watershed boundaries are no-flow ground water
boundaries. In other words, the hydraulic connection between the glacial and bedrock
aquifer is only maintained through the vertical flow between two aquifer layers. Based on
this assumption and the evaluation of the simulated heads, it is concluded that contrary to
common belief, in the southern Lower Peninsula, bedrock aquifers are well-connected
with the surface water network, especially where the glacial thickness is relatively small
and even in areas where the aquifer lithology is intercalated with confining geologic
th

th

formations. Specifically, highly permeable beds of the major (i.e. 4 to 6 order) streams
embedded in the glacial aquifer maintain a strong hydraulic connection with bedrock
aquifer, and control the ground water flow in the aquifer (Figure 5). This hydraulic
connection can be observed from Figure 7b where the simulated water levels have a very
good fit to data, and bend around and follow the streamlines (Figure 7a).

26

Simulated head (m)

Measured head (m)
DR=144 m
AME=-1.5 m
NSDR (%)=3.4
2

R =0.97
Stream
220

Pieozometric head (m)

a)
b)
Figure 7: a) Distribution of the simulated piezometric heads, Shiawassee Watershed, and b)
Simulated and measured piezometric heads.

In some watersheds across the southern Lower Peninsula, ground water is located in the
low-yielding shale formations. Water wells are drilled in the cracks of these formations where
the glacial thickness is relatively smaller and the bedrock aquifer is closer to the land surface.
Based on the evaluation of the simulations, it is concluded that the strong hydraulic connection
between the surface water network and shale formations also exists through the vertical flow
between two aquifer layers (Figure 8).

27

Simulated head (m)

Measured head (m)
Stream

220

Shale

Simulated head (m)
Bedrock well

DR=141 m
AME=2.7 m
NSDR (%)=5.1
2

R =0.93
a)

b)

Figure 8: a) Simulated piezometric heads in the shale formations, St. Joseph Watershed, and b)
Simulated and measured piezometric heads.

-

No-flow (i.e. Neumann) ground water boundary conditions can be used on the watershed
boundaries to model the ground water flow in both aquifer layers, except for several
watersheds that are located on the edge of Michigan’s two major uplands in northern and
southern Lower Peninsula (Figure 4). These two uplands represent the regional-scale
recharge areas and there is a significant amount of flow across watershed boundaries
from topographic highs to the discharge zones located at lower elevations. To reduce the
model deviations that occur due to the watershed boundaries, a careful examination is
needed. Specifically, numerical model boundaries need to be expanded to the extent for
which the real no-flow ground water basin boundaries are maintained (Figure 9).

28

Simulated head (m)

Measured head (m)

Original model (watershed) boundary

Comparison for the original model

Extended model boundary

* Comparison for the expanded model

Stream
a)

b)

Figure 9: a) Model boundaries before and after the boundary expansion, and b) Comparison
between the simulated and measured water levels before and after the boundary expansion.

-

In general, simulating glacial aquifer flow with one computational layer causes variable
magnitudes of model under predictions. Specifically, ground water flow in the
watersheds located on the uplands must be modeled with caution because of the presence
of the significant vertical flow. In these watersheds, simulated water levels represent the
average values over a large aquifer depth and are not comparable to the measured water
levels in the water wells (Figure 10). To capture more details vertically, the
computational layer needs to be divided into multiple hydrostratigraphic layers.
Similarly, since the depth of the wells in this region range between 30-45 m, hydraulic
conductivity data used in the models only represent a small portion of the glacial
deposits. As a result, any non-aquifer materials that may exist at deeper locations are not

29

included in the conductivity estimation which may potentially cause model under

Simulated head (m)

predictions in the region.

DR=294 m
AME=20.4 m
NSDR (%)=6.5 m
2

R =0.91

Measured head (m)
Figure 10: Simulated and measured water levels in the Manistee Watershed, glacial aquifer.

-

The models can be used to identify the watersheds for which the modeling system is
insufficient to simulate the ground water flow. For example, along the northern tip of the
Lower Peninsula, the bedrock aquifers in the lowlands of the watersheds are composed of
karstic limestone and dolomite geologic formations (Rheaume, 1991). Through a realtime analysis, it was found out that the two-layer conceptual models that adopt Darcy’s
Law for flow simulation do not apply to these watersheds. Because of this, there is a
significant bias between the simulated and measured water levels (Figure 11).

30

Simulated head (m)

Simulated head (m)

Measured head (m)
Measured head (m)
a)
b)
Figure 11: Simulated and measured water levels in the bedrock aquifer along the northern tip of
the Lower Peninsula: a) Black Watershed, and b) Cheboygan Watershed.

-

In some watersheds (i.e. Au Sable, Manistee, Muskegon, Pere-Marquette, Pine) located
in the northern Lower Peninsula, because of the high glacial thickness, bedrock aquifer is
rarely used for ground water supply (Rheaume, 1991). Since the number of data is very
few in these watersheds, the performance of the modeling system in the deep bedrock
aquifer is not validated.

-

Relatively speaking, in the Lower Peninsula, glacial aquifers have more spatial
heterogeneity than the bedrock aquifers. Because of this, data definitions of the
parameters that have an impact on flow variations are smoother in the bedrock aquifers
than that of the glacial aquifers in a way that helps the models predict the ground water
levels better in the bedrock aquifers. Head comparisons and model statistics also indicate
that defining effective transmissivity fields for bedrock aquifers through a collection of
different geologic sub-crops is sufficient for modeling regional-scale bedrock flow in
most watersheds.

-

Data in MSGWD can be relatively crude and involve different sources of uncertainties
such as driller variability, inaccuracies in well locations, and processing and recording
31

errors. Therefore before evaluating the simulations, certain queries (e.g. excluding the
water-wells for which the measured water levels are lower than the well’s bottom
elevation, excluding negative static water levels) are applied to eliminate particular data
errors. On the other hand, it is not possible to detect all “bad data” in the database, and
these data adversely affect the goodness-of-fit of the models.

Base flow comparison.

The purpose of the base flow comparison is to show that

the models presented in this paper can reproduce the USGS-estimated base flow values within
each watershed, after testing model parameters for a selected set of watersheds.

In the models, the source of water that enters the glacial aquifer layer is ground water
recharge from precipitation. This recharge is eventually discharged from the glacial aquifer as
base flow into streams. To define the connection between ground water and the streams, and to
match the USGS-estimated base flow values, different sets of leakance values were tested in
selected watersheds. In nature, stream width and streambed conductivity vary frequently along
the course of a stream so that it becomes practically very difficult to measure these parameters.
Because of this, the leakance values were introduced through a lumped parameter, L ' w , where

L ' represents the streambed leakance, and w represents the stream width. Different sets of values
were tested for the L ' w parameter based on the goodness of the comparison between the modelsimulated and USGS-estimated base flow values. Once an optimum set of values was
determined, those values were applied to all watershed models. In other words, the leakance
parameter was not thoroughly calibrated for each watershed but adjusted to provide an
acceptable prediction of the USGS base flow estimates in the watersheds. The resulting set of
L ' w varies depending on the stream orders (Table 7).
32

Table 7: Adjusted L ' w values.
Stream order

L'w (m/d)

1

3

2

5

3

10

4

20

5

50

6

100

To check the accuracy of the base flow calculations for all watersheds, the modelsimulated base flow values are compared to the USGS-estimated values (Figure 12). The
2

coefficient of determination (R ) for this comparison is found to be 0.98. This value indicates

Predicted Base flow (m3/s)

that the same set of L ' w values can be applied to all watershed models.

60
R² = 0.98

50
40
30
20
10
0
0

20

40

60

USGS Base flow (m3/s)
Figure 12: Correlation between the simulated and USGS-estimated base flow values.

33

st

nd

Additionally, it was found out that modeling the small (i.e. 1 and 2

order) streams as

two-way head-dependent flux features lead to gross water budget errors (Figure 13). Instead,
conceptualizing all streams in the watersheds as effluent streams is a justifiable assumption since
it provides accurate base flow and seepage discharge calculations.

Streams modeled as one-way head-dependent flux
Streams modeled as two-way head-dependent flux

Seepage discharge (m3/s)

60.0
50.0
40.0
30.0
20.0
10.0
0.0
Au Sable

Flint

Kalamazoo L. Grand Muskegon U. Grand
River
River

Figure 13: Comparison between the seepage discharge values when streams are modeled as oneway and two-way head dependent flux model features.

CONCLUSIONS

Li et al. at MSU recently designed a modeling system that maximally benefits from the
data-rich structure of MSGWD and that maintains the coupling between the data and ground
water models. In this paper, we systematically test the capability of this system to simulate the
regional-scale ground water flow in most watersheds across Michigan’s Lower Peninsula.
34

Overall, our investigation shows that the modeling system represents the regional-scale aquifer
dynamics and significantly enhances our ability to understand the critical hydrogeologic aspects
of the ground water systems in the peninsula. Model statistics show that the AME and NSDR are
generally less than 10 percent. The values of R2 are usually above 85 percent. For the watersheds
where the model performances are relatively low, the comparisons between the simulated and
measured ground water levels gave us the opportunity to visually and interactively examine the
accuracy of the assumptions used in the conceptual model.

Specifically, the comprehensive modeling and evaluation effort presented in this research
provides the following messages:

-

The relatively simple two-layer models are able to capture the dominant variations in
regional-scale ground water flow, including the interaction between the upper glacial and
lower bedrock aquifers.

-

Contrary to common belief, the bedrock aquifers are well connected with the surface
water network, even in the watersheds where confining layers exist, and generally the
th

th

high (i.e. 4 to 6 ) order streams that are linked to the glacial aquifer control the bedrock
flow.
-

No-flow boundary conditions can be imposed on the watershed boundaries to simulate
ground water flow in the glacial and bedrock aquifers, except for the watersheds around
the northern and southern uplands (or the regional-scale recharge zones) of the Lower
Peninsula, due to the substantial amount of ground water flow from high to low
elevations.

35

-

Simulated heads in the glacial aquifers on the uplands represent the average heads over
large aquifer depths, and are not comparable to the measured static water levels in the
region, because of the significant amount of vertical flow.

-

Performance of the modeling system in the deep bedrock aquifer underneath the uplands
is not validated because of lack of measured data.

-

Modeling the small streams as two-way head-dependent flux features lead to gross water
budget errors, and these streams should be represented as effluent streams.

-

Darcy’s Law cannot be applied to simulate the ground water flow in the bedrock aquifers
along the northern tip of the peninsula, because of the karstic nature of the bedrock in this
region.

The GIS-enabled ground water models presented in this research significantly reduced
the time and financial costs needed to collect field data, and shifted the focus to conceptual
modeling and critical thinking. This provides a very useful tool that can be utilized for costeffective ground water management in Michigan. In the future, results of this research can be
used to calibrate the simulations to the static water levels from the well records and the stream
flow data from the USGS gauging stations for site-specific ground water problems.

36

APPENDICES

37

Appendix A: Michigan’s Statewide Groundwater Database

Michigan’s statewide groundwater database (MSGWD) consists of four main data
categories:



Geological/hydrological framework which provides the spatial distributions of the aquifer
parameters



Physical framework which helps to define the geometry of the aquifers



Groundwater systems’ stress factors which are imposed to the models through the source
and sink terms



Calibration parameters which are used for the model validation

MSGWD uses the Michigan Georeference projected coordinate system and the GCS
North American geographic coordinate system. Unit systems used in the database are the US
Customary System and the International System of Units.

Data categories are available in two different geographical model formats:



Raster data models which appear in finite difference models and represent twodimensional spatial variations by assigning a single value to each cell in a fixed
rectangular array. For example, a layer might represent the spatial variation of depth to
groundwater in a layer constructed by conducting interpolation between data located at
different measurement locations (Goodchild, 1996).



Vector data models which represent spatial variation through irregularly distributed
points, lines, and polygons. The meaning of layer for vector systems is more complex
38

than that of raster data models. While mainly a vector layer provides a parameter (e.g.
hydraulic conductivity, SWL, stream elevation) to a groundwater model; it can also
facilitate display (e.g. county boundary, road) by grouping one or more classes of points,
polylines, or polygons (Goodchild, 1996).

Current trend in water resources management in the United States is to develop a
hydrological data network (Consortium of Universities for the Advancement of Hydrologic
Science, CUAHSI) that could be accessed by anyone through web services using a stand-alone
computer. MSGWD is readily available to be connected to this system since the state’s
groundwater data can be downloaded and mapped through the state’s online geographic
information system (http://www.michigan.gov/cgi) and the GWIM’s webpage
(http://gwmap.rsgis.msu.edu/) which is the result of a collaborative effort between the MDEQ,
the USGS, and Michigan State University. Most important advantages of these online services
are that all data provided through the services have consistent data structures and that these data
can be directly visualized and processed through GIS tools. Data that were produced for the
GWIM are briefly discussed in the following sections.

Hydrogeological Setting of Michigan

The regional aquifer system of Michigan’s Lower Peninsula consists of the
unconsolidated glacial aquifer, and the two major bedrock aquifers (Saginaw and Marshall
Aquifer) separated by confining layers. The Coldwater Shale forms the base of the regional
aquifer system (Barton et al., 1996; Westjohn and Weaver, 1998).

39

The glacial deposits cover nearly all of the Lower Peninsula. The glacial thickness map
of the Lower Peninsula is shown in Figure 14. This map is generated by subtracting the rasterbased bedrock surface map (with a resolution of 500 m) provided with the GWIM from the
raster-based USGS 90 m SRTM land surface map.

Glacial Thickness (m)
376

1
Scale is non-linear

Figure 14: Glacial thickness map for Michigan's Lower Peninsula. Values are subject to change
due to the resolution of the raster-based data.
40

The glacial deposits are the largest reservoir of fresh groundwater in the Lower Peninsula
and primarily consist of thick sequences of sand and gravel, but in some areas also consist of
sand and gravel layers within till or other fine-grained deposits (Barton et al., 1996; Westjohn
and Weaver, 1998).

The Saginaw Aquifer consists of sandstone lenses, and typically ranges in thickness from
100 to 350 ft in areas where this unit is used for water supply. In the western-central part of the
aquifer system, the Saginaw Aquifer is separated from the glacial deposits by 100 to 150 ft of red
beds. The red beds consist of mudstone and poorly consolidated red shale, gypsum, and minor
amounts of sandstone. Together with fine-grained glacial and lacustrine deposits, the red beds
form the sub-regional confining units that separate the glacial aquifer from the underlying
Saginaw Aquifer (Barton et al., 1996; Westjohn and Weaver, 1998).

The Marshall Aquifer, the basal aquifer in the basin, includes sandstones that overlie the
Coldwater confining unit. The Marshall Aquifer consists of two or more permeable sandstones in
the central part of the basin; but in and near areas where the aquifer is a sub-crop of the glacial
deposits, it consists of one permeable sandstone unit. In areas where more than one sandstone
unit exists, alternating carbonate, shale, siltstone, and evaporites separate permeable sandstones.
The thickness of permeable sandstones is typically from 75 to 200 ft. The Marshall Aquifer is
regionally confined except where it is a sub-crop of the glacial aquifer (Barton et al., 1996).

The bedrock aquifers (Traverse Group) along the northern tip of the Lower Peninsula are
covered with relatively thin glacial lacustrine deposits. Traverse Group consists of limestone and
dolomite formations which are the primary source of most groundwater supplies in the area. The

41

transmissivities of these karstic aquifers depend on the interconnection of fractures and solution
channels in the aquifers (Rheaume, 1991).

Bedrock layer is exposed in most parts of the Upper Peninsula. Here, major aquifers are
the crystalline-rock aquifer and the Jacobsville Aquifer. The cyrstalline-rock aquifer is located in
the western part of the peninsula. Groundwater generally moves through secondary openings,
such as joints, fractures, or faults, in the cyristalline rocks. The Jacobsville Aquifer is found
running under large portions of the Upper Peninsula along the Lake Superior shore. It is
composed of feldspathic and quartzitic sandstone and shale formations. Although the aquifer
consists of sandstone, it is well cemented and has a low permeability. Groundwater moves
through joints and fractures which decrease to insignificance at depths of 100 and 150 ft. The
aquifer is considered a confined aquifer (Marsicek, 2002).

Approximate boundaries of the major bedrock aquifers beneath the surficial glacial
aquifer are shown in Figure 15.

42

Red Beds
Saginaw Aquifer
Marshall Aquifer
Traverse Group
Coldwater Shale
Jacobsville Aquifer
Cyristalline-Rock Aquifer
Non-aquifer material

Figure 15: Major bedrock aquifers in Michigan.

Wellogic

The paper logs of the water wells that had been drilled in Michigan since 1960’s were
scanned and processed by the MDEQ to make almost 300,000 (totally in the glacial and the
bedrock aquifers) data points (as of 2006) usable in Wellogic. MSU and the MDEQ
systematically work together to update the database with the most recent well data. The data is
stored in the shapefiles where the wells are recorded with a unique Well ID along with the
43

attributes that summarize the location of the wells and the spatial data to provide the information
for the hydrogeological and the physical framework of Michigan’s aquifers.

Number of the water wells in Wellogic increases by approximately an amount of 1,500
every month (GWIM, 2006). Even though data resolution is not uniform throughout the state,
since the amount of the data in the database is steadily increasing, in the near future, in the areas
where the data amount is poor, resolution will be much better than that of today. Change of the
number of the statewide data with respect to time is illustrated in Figure 16 and Figure 17.

Wellogic comes with the following advantages:



Data is available almost everywhere in the state, free, and open to public



If made full use of it, massive data structure reduces both the monetary costs and the time
required to conduct field investigations



Data resolution is significantly higher than that of the traditionally collected data



Even though it is not uniform and because of the uncertainties involved, it is not as good
as of the traditionally collected data, data quality improves with time (Figure 18) shows
the increase in the number of the data points from 1970 to 2006 and how this
improvement translates into the map of the hydraulic conductivity field for the
Kalamazoo County, MI.

44

120,167 wells

197,548 wells
1970

Glacial Well

1980

14,001 wells

1990

52,762 wells

2000

93,421 wells

2006

Figure 16: Time-based spatial distribution of the water wells in the glacial aquifers in Michigan.

45

47,125 wells

70,312 wells

6,367 wells
1970

1980

Bedrock Well

1990

23,367 wells

2000

36,289 wells

2006

Figure 17: Time-based spatial distribution of the water wells in the bedrock aquifers in
Michigan.
46

K (m/d)
38.4

9.6
Scale is non-linear
a) 264 data points (1970)

K (m/d)
43.2

3.5
Scale is non-linear
b) 11,298 data points (2006)
Glacial well
Figure 18: Time-based change in the representation of the hydraulic conductivity (K) field in the
Kalamazoo County.

47

Physical Framework

Land surface elevation.

All wells recorded in Wellogic have an attribute that

provides the land surface elevation at the well’s location. This attribute was generated from the
raster-based 30 m (7.5 minute) Digital Elevation Model (DEM) provided in the USGS National
Elevation Dataset, NED. For the GWIM, following steps were applied to extract the elevation
data for each well by using ArcGIS 8.3 (GWIM, 2006).



NED coverage of the state was downloaded in multiple geographical tiles from the
USGS’s related website



The raster data in the tiles were converted into integer grids



Tiles were merged into a statewide DEM map



The statewide map was projected into Michigan’s coordinate system and the map unit
was converted into decimal meters



The unit of the data was converted into ft and these new values were added to the
database



The statewide map was divided into the counties (totally 83) with an exterior buffer of 1
mile to the outside of the county boundaries



The county-based NED data grid was used to sample the land surface elevation at the
location of the water wells

More detailed information for the NED and a discussion of the method that is used to generate
the raster-based DEM database for this research can be found in Appendix C (p. 89).

48

Surface elevation of the bedrock layers.

Statewide oil and gas wells (i.e. deep wells)

compiled by the MDEQ, and the water wells which intersect the bedrock layers where the glacial
thickness is small (the Upper Peninsula or the southern Lower Peninsula), were used to map the
surface of the bedrock aquifer (GWIM, 2006). Total number of the wells that were used for
mapping was approximately 53,000 (Figure 19).

Oil/Gas Well

Figure 19: Distribution of the oil and gas wells in Michigan.

49

In the point-based well shapefile, the maximum recorded value of the surface elevation of
the bedrock layer is 1,430 ft where the minimum is -982 ft. A raster map with a 500 m of
resolution was generated by interpolating the information provided in the database (Figure 20).

Bedrock Surface Elevation (m)
585

-10
Scale is non-linear

Figure 20: Raster-based bedrock surface map of Michigan.

A discussion on the existing bedrock surface map can be found in Appendix B (p. 85).

50

Watershed boundaries.

The USGS delineated the watersheds in the United States

using a nationwide hydrological system that divides and subdivides the country into hierarchical
levels. These levels of sub-divisions, which are used for organization of the hydrological data,
are called hydrological units. Hydrological units are used to collect and organize the data and
these units can be defined as the land areas that catch precipitation and drain it to rivers, streams,
wetlands, lakes, or groundwater (Seaber et al., 1987).

The largest unit in the hydrological unit system is called region and represents either the
drainage area of a major river or the combined drainage areas of a series of rivers. Watersheds
are the smallest subdivisions of the regions. Each watershed includes the area drained by a river
system or a section of a river (i.e. stream) and its tributaries in that reach (Seaber et al., 1987).

The watersheds are identified with an 8-digit identification number which is called
hydrological unit code (HUC). A hydrological unit code consists of 2 digits for each
hydrological level that describes the relation of the hydrological units to each other to represent
the way the smaller watersheds get together to form the larger watersheds. The hydrological
system has a scale of 1:2,000,000 (1 cm = 20 km) and divides the country into 2,264 watersheds,
60 of which are located in Michigan (Figure 21). The USGS provides not only the boundaries
but also the surface areas (in acres) of the watersheds in the corresponding shapefiles (Seaber et
al., 1987).

51

Watershed Boundary
State Boundary
Major River

Figure 21: Watersheds of Michigan.

The USGS also provides a statewide GIS database for sub-watersheds that divides the
existing watersheds into smaller drainage areas. This GIS database was generated from the
1:24,000 (1 cm = 240 m) scale topographic maps which are also known as 7.5-minute
quadrangles. Totally there are 2,287 sub-watersheds in Michigan’s sub-watershed system

52

(Seaber et al., 1987). A set of sub-watersheds located in the Muskegon Watershed is shown in
Figure 22.

Watershed Boundary
State Boundary
Stream/River

Figure 22: Sub-watersheds of the Muskegon Watershed.

Hydrogeological Framework

Hydraulic conductivity of the glacial aquifer.

In GWIM, equivalent horizontal

hydraulic conductivity for each water well drilled in the glacial aquifer was estimated (in ft/d)
from the lithology information reported in Wellogic (GWIM, 2006). Two hydraulic conductivity

53

values were calculated for each well: Horizontal Hydraulic Conductivity 1 (HCOND1) and
Horizontal Hydraulic Conductivity 2 (HCOND2).

HCOND1 is defined as the equivalent hydraulic conductivity of the lithologies from the
bottom of the well to the first confining layer (Figure 23). However, HCOND1 recorded in the
system is most of the times equal to the lithology in the well screen. Since the well screen is built
within the most conductive geological material, magnitude of HCOND1 is most of the times has
2

ranges between 10 and 10 .

Water table

HCOND1

K1, B1

calculation (screen) interval

Bottom of the well
Figure 23: Conceptualization of the calculation interval for the HCOND1 parameter.

Calculation of HCOND2 is different than that of HCOND1. HCOND2 value was
calculated based on the lithology from the water table to the bottom of the well or the top of
bedrock for the wells completed in the bedrock aquifers underlying the glacial aquifers. In other
words, HCOND2 represents the hydraulic conductivity of the saturated thickness (Figure 24).

54

Water table

K1, B1
HCOND2

K2, B2

calculation (screen) interval
K3, B3
Bottom of the well
Figure 24: Conceptualization of the calculation interval for the HCOND2 parameter.

Assuming that all lithologies within the wells are layered horizontally, equivalent
horizontal hydraulic conductivity in the wells was calculated with the thickness-weighted
average formula [Equation (13)] (Freeze and Cherry, 1979):

n

 Ki Bi

K h  i 1

13

B

where
L
K h : equivalent horizontal hydraulic conductivity in the water well ,  
T 
n : number of lithologies within the calculation interval ,
L
Ki : hydraulic conductivity of the lithological unit ,  
T 
Bi : thickness of the lithological unit ,  L 
B : saturated thickness,  L 

55

In general, the equivalent horizontal hydraulic conductivity is controlled by the lithology
with the highest hydraulic conductivity. The hydraulic conductivity values that were assigned to
the lithologies had been adopted from the Source Water Assessment Program conducted by the
MDEQ in 2004 by using a textbook range of values. For each lithology, three different
(minimum, intermediate, maximum) hydraulic conductivity values were defined. Depending on
the land system that the water wells were drilled in, one of the three values was assigned for the
corresponding lithology. Full list of these values are available in the GWIM Technical Report
(GWIM, 2006). Some of the selected values along with the land system classification are given
in Table 8.

For some wells, equivalent vertical hydraulic conductivity value was also calculated
depending on the availability of data, and using [Equation (14)] (Freeze and Cherry, 1979):

Kv 

m

B
Bk

14 



k 1 K k

where
L
K v : equivalent vertical hydraulic conductivity in the water well ,  
T 
m : number of lithologies within the calculation interval ,
L
K k : hydraulic conductivity of the lithological unit ,  
T 
Bk : thickness of the lithological unit ,  L 
B : thickness of the calculation interval ,  L 

56

Table 8: Hydraulic conductivity classification based on the land systems of Michigan (GWIM,
2006).
Land system

Lithology

Bedrock
Lacustrine fine
Lakes

Ice-marginal till
Lacustrine coarse
Lodge till/fine
supraglacial drift

Coastal dunes
Ice-contact
outwash
Proglacial
outwash

Hydraulic conductivity (ft/d)
Minimum

Intermediate

-4

Clay
Clay & Gravel

3x10

Clay & Sand

3x10

Clay & Silt

10

-4

10

Gravel

10

-4

-3

10

-4

-2

10

-4

-4

Maximum
-4

10

-2

10

-2

10

-4

10

100

50
-2

10

300

-2

-2

Gravel & Clay

10

Gravel & Sand

1

50

100

Gravel & Silt

1

1

1

Sand

1

50

100

10

0.1

1

Sand & Gravel

1

50

100

Sand & Silt

1

1

10

1

1

Sand & Clay

10

-4

-1

Silt

10

Silt & Clay

10

Silt & Gravel

10

Silt & Sand

10

-4

-3

10

-2

10

10

-1

1

1

-1

1

1

57

The conceptual difference between the HCOND1 and the HCOND2 parameters is
demonstrated in Appendix B. Final products that are generated to implement HCOND1 and
HCOND2 into the groundwater models developed in this research are discussed in Appendix C.

Transmissivity of the bedrock aquifer.

Transmissivity map of the bedrock

layer was generated by interpolating a total number of 233 point-based transmissivity data
(GWIM, 2006). The sources of the transmissivity data used in the interpolation are given below:



The MDEQ collected the paper records of the aquifer tests conducted by the consulting
companies in various municipal wells across the state



The USGS compiled the aquifer test records throughout the state as a part of the
Michigan Basin Regional Aquifer System Analysis (RASA) Program

Existing transmissivity map of the bedrock aquifers in the state is improved for this
research. Further discussion on the improvement can be found in Appendix C (p. 85).

Stress Factors

Surface water network.

MSGWD uses the medium resolution (1:100,000) USGS

NHD for the identification of Michigan’s surface water network. The NHD provides
comprehensive digital spatial data that represent the surface water network of the United States
through common features such as streams, rivers, lakes, and ponds. These data are designed to be
used in general mapping and the analysis of surface water systems and surface watergroundwater interactions through GIS. For mapping and analysis purposes, the NHD is used with

58

other data themes such as the surface elevation and feature boundaries to produce general
reference maps (USGS, 2000).

The features in the NHD were organized into polylines and polygons. The stream lines
were broken up into shorter segments stretching from confluence to confluence and these
segments were linked together to trace the flow of water. Artificial lines that go through inside
lakes represent the flow of water through a lake. The polygons typically represent the water
bodies such as lakes and ponds where the lines were used to digitize the streams (USGS, 2000).

Streams were labeled with a unique and permanent identifier (reach code) which gives
the streams an identity for inventory and analysis. This makes it easier to perform calculations in
GIS such as determining the distance to the stream flow gauging stations and calculating base
flow. Linear referencing also makes it easy to link additional data to the NHD without having to
customize it (USGS, 2000).

The NHD is available for different resolutions. In this research, a medium resolution
NHD which was generated with 1:100,000 (1 cm = 1 km) scale topographic mapping is used for
groundwater modeling (USGS, 2000).

The NHD GIS files were updated by appending quite a few different parameters for every
stream segment and lake recorded in Michigan (GWIM, 2006). However, number of the essential
parameters for modeling is limited to a few. Table 9 shows the attribute names and the
definitions of these parameters.

59

Table 9: Essential modeling parameters recorded in the NHD.
Feature

Parameter

Definition

Unit

Stream

Order (Strahler
number)

Stream size based on the hierarchy of the
tributaries (available 1 through 6)

-

F_Elev

Upstream elevation of the stream segment

m

T_Elev

Downstream elevation of the stream segment

m

Base flow

Groundwater discharge to streams

ft /s

Name

Stream name

-

Length

Length of the stream segment

m

Lake_Elev

Lake elevation

m

Area

Surface area of the lake

ac

Name

Lake name

-

Lake

3

In the NHD GIS files, for some stream segments a generic value of -9,999 for the
attributes that define the upstream and the downstream elevations of those segments is detected
(Figure 25). This value indicates the stream segment locations where there is no data available
for the stream stages. Stream stages in the vector-based stream segments were generated by
extracting the land surface elevations from the USGS 30 m NED (GWIM, 2006). For this
research, statewide raster-based stream stage maps are used to directly extract the stage values.
Final product for the stream network used in the groundwater models is discussed in Appendix C
(p. 96).

Lately, the USGS, in cooperation with the U.S. Environmental Protection Agency
(USEPA) evaluated different methods to delineate catchments of the medium resolution

60

(1:100,000-scale) NHD with a goal to estimate stream flow and velocity to support water quality
modeling. The results of this study are the first national dataset that links the mapped stream
network to the landscape. The data, called NHDPlus, are distributed for all states, except Alaska,
through the website http://www.epa.gov/waters (Johnston et al., 2009).

Stream segments with missing data

Figure 25: Stream segments with missing stream stage data in the NHD of Michigan.

Estimated shallow groundwater recharge.

In this research, groundwater

recharge is defined as the precipitation that infiltrates the land surface, moves downward through

61

the unsaturated zone, and enters the water table. The long-term average base flow values were
assumed by the USGS to be equal to the long-term average groundwater recharge (i.e. total water
entering the soil minus losses to evapotranspiration before it reaches the water table) in the
surrounding watershed. Steps for calculating the estimates of recharge are given below (Neff et
al., 2006):



Long term stream flow records from the 162 stream flow gauging stations in the Lower
Peninsula and the 46 stream flow gauging stations in the Upper Peninsula out of totally
416 stations (Figure 27) were used to determine the stream flow component that results
from groundwater discharge (base flow) into the streams



Base flow was supplemented by direct runoff during and immediately after precipitation
or snowmelt events, resulting in peaks on a hydrograph showing stream flow through
time



Peaks on the hydrographs were divided into base flow and stream flow components by
using the hydrograph separation technique



Watersheds were delineated for each of the 208 gauging stations, and the state was
divided into three regions (i.e. the eastern Lower Peninsula, the western Lower Peninsula,
and the Upper Peninsula) based on the distribution of base flow yield (base flow per unit
area)



For the ungauged areas, multivariate linear regression models were used to relate the base
flow information in the gauged areas to landscape attributes (e.g. agricultural land use,
urban land use, annual growing degree days, annual precipitation, percentage of the
watershed underlain by the lacustrine deposits) in the surrounding watersheds for each of
the three regions
62



A base recharge (mean base flow yield) value was estimated within each of the three
regions by dividing the base flow to the total area of each of the three regions



As an initial estimate, mean base flow yield value was multiplied by the watershed area
for each station



A residual value was then calculated by subtracting the measured base flow at each
station from the area-based prediction



A second set of linear regression models was created to predict these residuals based on
the same watershed attributes used in the previous regression model



The resulting regression models were then applied to estimate the residual in each section
2

of 1 mi of area


Finally, the residuals were added to the initial estimated yields to come up with an
estimate of the total recharge for each section

The spatial resolution of the estimates of recharge in the Upper Peninsula is not as
detailed as of the Lower Peninsula because the number of stream flow gauging stations in the
Upper Peninsula is not good enough to provide a sufficient number of observations to support
the incorporation of land cover or surface geology into the models.
2

The estimated recharge values for each recharge polygon with an area of 1 mi were
tabulated into a shapefile (GWIM, 2006). For this research, this shapefile is converted to a rasterbased map to be used in the groundwater models (Figure 26). It is observed that, in the original
vector-based GIS file, there are local zones for which recharge data is missing. Locations and
extensions of these polygons are also shown in Figure 26.

63

Recharge (in/yr)
22

1
Scale is non-linear
No data

Figure 26: Raster-based recharge map of Michigan.

Calibration Parameters

Static water level.

Wellogic contains more than 45 years of water level data which

makes this data set ideal to demonstrate the long-term regional-scale behavior of the
groundwater level distributions across the state. The static water level (SWL, in ft) recorded in
the database is defined as the depth from the land surface to the water table in the corresponding

64

well. The SWL reported here is the water level in the well after the well is developed, prior to
pumping.

In this research, SWL information is used to conduct the comparison between the
groundwater levels simulated with the watershed-scale models and measured groundwater levels
in the water wells.

Base flow.

Base flow can be defined as groundwater discharge into streams.

Precipitation that is not returned to atmosphere by evapotranspiration either flows across the land
surface, discharging directly into wetlands, lakes, and rivers, or infiltrates into land surface,
recharging groundwater system and ultimately resulting in groundwater discharge that
contributes to base flow. The partitioning of precipitation between runoff and infiltration is
controlled by the capacity of groundwater system to receive and transmit water. This partitioning
is a result of landscape attributes (p. 61), and the characteristics of the shallow subsurface
materials which are the principal constraints in the estimation. Coarse-textured unconsolidated
sediments and permeable and fractured bedrock are typically associated with above average
infiltration, whereas fine-textured unconsolidated sediments and low-permeability bedrock are
typically associated with above average surface runoff (Neff et al., 2005; Neff et al., 2006).

The USGS used two fundamental models to calculate a base flow index (BFI, i.e. ratio of
base flow to stream flow) for the estimation of base flow in each watershed in Michigan (Neff et
al., 2005): (1) G-Model which develops an empirical relation (from stepwise regression
analysis), between each surface material class (bedrock, tills, organic sediments, coarse-textured
sediment, fine-textured sediment) and the BFI in each gauged watershed. Here, it was assumed
that groundwater discharge is most closely related to surface geology, and that the averaged
65

values of BFI approximate the portion of stream flow that is the result of groundwater discharge.
(2) G-SW model, which along with the geological interpretation used in G-model, takes the
surface water features into account by incorporating the proportion of surface water network
within a gauged watershed, and develops an empirical relation to calculate the BFI.

The USGS used the stream flow records from the stream flow gauging stations (Figure
27) to estimate the long term average stream flow in each HUC (here watershed) across the state.
Base flow was then determined by multiplying model-simulated BFI with the estimated total
stream flow for each HUC. Further details for the model equations and the parameter weights
used in the models can be found in (Neff et al., 2005).

66

Õ

Stream Gauging Station
Major Stream

Figure 27: Locations of the stream gauging stations in Michigan.

In gauged areas (i.e. areas that are tributaries to a given gauge or gauges), base flow
estimates were calculated based on the hydrograph separation using the stream flow data
extracted from the USGS Automated Data Processing System (ADAPS). To develop regression
models to describe base flow in ungauged areas (watersheds that are represented by 8-digit
HUCs which may also contain gauged areas), only the USGS stream flow gauging stations with
a length of record of at least 36 months were used (Neff et al., 2005).

67

Appendix B: Conceptual Difference between HCOND1 and HCOND2

As it was discussed in Appendix A, the HCOND1 values in the water wells represent the
lithology of the screen intervals where the HCOND2 values represent the saturated thickness in
the water wells (or the saturated thickness of the glacial aquifer where the aquifer thickness is
small). Difference between the values of the HCOND1 and the HCOND2 parameters is shown in
Figure 28.

68

K (m/d)
82.9

2.4
Scale is non-linear

K (m/d)
63.5

0.4
Scale is non-linear

Figure 28: Hydraulic conductivity field of the Tittabawassee Watershed. Above HCOND1, and
Below: HCOND2.

69

All models presented in this research utilize the HCOND2 parameter to generate the flow
fields. For demonstration purposes, the model of the Tittabawassee Watershed is redone using
the HCOND1 parameter, and the head comparisons from both conductivity fields are shown on
the same plot (Figure 29). In this figure, it can be clearly seen that the model generated with the
HCOND1 parameter predicts lower water levels (black dots) at the same locations than those of
the HCOND2 parameter (blue circles). Even this comparison alone indicates that the HCOND1
parameter cannot be used for aquifer characterization and modeling purposes. The HCOND1 can
be used for the applications where the screen interval of the water wells is of great importance in
terms of determining the flow direction when the wells are pumped (e.g. well-head delineation).

Figure 29: Comparison between the model results with HCOND1 and HCOND2. Black Dots:
Results for HCOND1, and Blue Circles: Results for HCOND2.

70

Appendix C: Data Analysis and Management

Introduction

Model prediction is always uncertain to some degree, and it is of critical importance to
have a reliable dataset for which data uncertainties are reduced so that the model’s ability to
represent the reality can be improved.

In MSGWD, while some parameters were presented as directly measured field data, some
of them were calculated from those directly measured data, or estimated through interpolation of
the existing data by using various geospatial modeling techniques (e.g. Kriging interpolation).
While directly measured Wellogic data has uncertainties such as temporal variability, location
inaccuracy, and driller variability, calculation and processing errors could also be introduced for
the indirectly calculated parameters.

Interpolation can be justified if the data is abundant (Lubczynski and Gurwin, 2005).
However, even in cases where data is abundant, generated maps come with problems that
originate from the data errors and uncertainties, and the methods that are used for mapping.
Because of all these potential problems, before developing the conceptual models of the state’s
watersheds, the data provided with the GWIM is revisited to check possible data problems, a step
that is important to improve the quality of the models. Also, where it is projected that data
smoothing won’t significantly change the regional-scale data definitions while the computational
efficiency will be improved, some vector data are converted into raster data by benefiting from
the Kriging and the inverse distance weighting interpolation techniques.

71

In this research, for data analysis and groundwater modeling, mainly two programs are
used: ArcGIS (version 9.2) and Interactive Groundwater (IGW, version 5.20P). IGW is a realtime, interactive, and hierarchical simulation environment with a GIS interface that allows data
integration for deterministic, stochastic, and multi-scale finite difference groundwater flow and
contaminant transport modeling. IGW has been developed and continuously improved in the
Laboratory of Excellence for Real-time Computing and Multi-scale Modeling (LERCMM),
Department of Civil and Environmental Engineering, Michigan State University. Other basic
specifications of the program are summarized in Table 10.

Table 10: Basic specifications of IGW.
Developers: Shu-Guang Li, Huasheng Liao, Qun Liu
Contact: A133 Engineering Research Complex, Michigan State University, East
Lansing, MI 48824. E-mail: lishug@egr.msu.edu
Year first available: 1997
Hardware required: IBM compatible PC
Operating system: PC with Windows (NT/95/98/ME/2000/XP)
Programming language: Visual Basic 6, Visual FORTRAN 5
Space required for installation: 50 MB
Current version: IGW 5.2P
Availability: Downloadable version (version 4.7), with user’s manual and
supporting material, available at http://www.egr.msu.edu/igw/

Specifications of the computer that was used for most of the research are given in Table
11.

72

Table 11: Basic features of the main computer used in the research.
Processor: Pentium D CPU 3.20 GHz
Memory: 3.50 GB
Operating system: Microsoft Windows XP Professional (32-bit)

For more challenging analysis and modeling issues that needed better computational
resources, the computer of which specifications are given in Table 12 was used:

Table 12: Basic features of the advanced computer used for some of the simulations.
Processor: Intel Xeon CPU X5482 3.20 GHZ
Memory: 8 GB
Operating system: Microsoft Windows XP Professional (64-bit)

Identifying and Eliminating the Data Uncertainties

Data uncertainties in Wellogic.

Wellogic provides valuable information regarding

the hydrogeological and the physical framework of Michigan’s groundwater systems. On the
other hand, efficiency of groundwater flow and contaminant transport models is significantly
influenced by the “bad data” and uncertainties involved in the data (Simard, 2007; Tiedeman et
al., 2003). This fact comes along with the requirement of a systematic analysis for the
identification of the data uncertainties which helps to improve the groundwater management.

Some of the uncertainties that are introduced with Wellogic are given below:

73



Temporal variability: The wells were drilled in different times and measured timedependent variables (e.g. SWL) only reflect the conditions at the drilling time. However,
at a regional scale with the help of large amount of data, temporal variability is averaged
and gives long-term trend that may be used in steady-state groundwater models.



Vertical uncertainty: The data recorded (e.g. SWL, hydraulic conductivity) in the water
wells reflect the conditions in the well screen. For the groundwater models developed in
this research, the spatial data in the water wells are considered as they represent the
vertical average values for the corresponding aquifer (i.e. glacial or bedrock) layer.



Well location inaccuracy: The quality and location of Wellogic data is generally worse
than that of the traditionally collected data especially because of the well locations that
were recorded based on address matching, before the global positioning system
technologies became popular.



Driller variability: The water wells in the database are drilled by different drillers.
Because of this, data quality from different drillers is different.

Vertical uncertainty discussed above can be explained by an example regarding the
definition of the HCOND2 parameter provided in MSGWD. While HCOND2 is a good
estimation for the saturated thickness in the water wells, where glacial deposits extend below the
water wells and glacial thickness is extremely high (e.g. northern part of the Lower Peninsula),
deeper lithologies are not incorporated into HCOND2.

In addition to the uncertainties discussed above, there are systematic inconsistencies, bad
and missing values, and recording, processing, and typographic errors, along with wells with no
data for specific parameters. Even though the number of these data points is not too many, it is
74

necessary to systematically eliminate them to improve the quality of the data-based analysis and
the model results. Data elimination criteria and the number of the data eliminated using those
criteria are shown in Table 13. The K data with zero and negative values don’t exist in the
current (2006) version of Wellogic, and the elimination criteria of K ≤ 0 is added as a precaution.

Table 13: Criteria used for eliminating the bad data.
Data Elimination Criteria

No. of Wells
Glacial

Bedrock

SWL ≤ 0

5196

2630

SWL ≥ 999

-

1680

DEM = 0

120

95

(DEM – SWL) ≤ 0

-

1640

(SWL – Well Depth) ≥ 0

97

77

K≤0

-

-

Another observation regarding the data uncertainties is that for some water wells, well
depth data with zero values exist, where the SWL data for those wells are recorded as positive
values i.e. (SWL can’t be bigger than well depth). This doesn’t necessarily mean that the
recorded SWL data for those wells are wrong, however, data inconsistencies must be recognized
for this kind of cases with different comparable data combinations.

At this point, it is important to revisit the land surface elevation data provided with the
water wells. The accuracy and reliability of the land surface data are very important since this
parameter represents the upper boundary of the glacial aquifer, and it is needed to calculate the

75

measured heads above sea level. Furthermore, the stream stages in the 2006 version of Wellogic
were derived from the land surface elevation data stored in the water wells.

To test the accuracy of the land surface elevation data in the water wells provided with
the GWIM, following procedure is applied:



First, the most recent raster-based USGS DEM (p. 89) with 10 m resolution (sampled
from 30 m DEM) is used to extract the DEM data at each well location in Wellogic



Then, a new shapefile, with the attributes of the DEM data from the GWIM and the DEM
data extracted from the new datasets at the same well locations (almost 300,000 data
points), is produced



Finally, in the resulting shapefile, the data from the GWIM is subtracted from the new
data to calculate the residuals



For the areas where there are huge differences between the two datasets, the new DEM
data is cross-validated by comparing the surface elevation data stored in the oil and gas
wells

Residual map produced through the procedure discussed above is shown in Figure 30.

76

DEM Error (m)
-187

123
Scale is non-linear

Figure 30: Error map for the DEM information in Wellogic.

Figure 30 shows that the residual between the two datasets can produce extremely big
values at some locations. Final product provided with the GWIM doesn’t directly use the raster
data but instead converts it to point-based vector data that is defined at every water well location,
after multiple procedures that go back and forth. Creating such a middle product, where the
locations of the wells have also uncertainties, could introduce a new source of error. In order to
demonstrate how the error in land surface elevation can adversely influence the simulated heads
77

(i.e. land surface elevation minus depth to water table), a test site is chosen where the DEM error
becomes the highest (Figure 31).

Here, while the water table calculated from the point-based land surface data that was
generated from the USGS raster-based 30 m DEM dataset produce extremely bad values, water
level elevations calculated directly from the raster-based DEM data with the same resolution
provides a much better final product (Figure 32).

78

KALAMAZOO

BARRY
DEM Error (m)
-187

123
Scale is non-linear

VAN BUREN

ALLEGAN

KENT
OTTAWA

Figure 31: Location of an area in the western Lower Peninsula where DEM error is extremely
high.

79

Figure 32: Simulated and measured heads for the model area shown in Figure 31. Above:
comparison using vector-based DEM, and Below: comparison using raster-based DEM with 10
m (sampled from 30 m) resolution.

80

Another example that helped to identify the data uncertainties through a reverse (modelbased) analysis is shown in Figure 33.

Figure 33: Identifying bad SWL data using model-based analysis, Thornapple Watershed,
bedrock aquifer.

Figure 33 demonstrates the head comparison between the simulated and measured heads
in the bedrock aquifer in the Thornapple Watershed (southern Lower Peninsula). The strong bias
detected on the left side of the forty-five degree line indicates that there might be an uncertainty
that needs to be examined by going back to the GIS files and look for a possible source of error.
At the end, the bias in the comparison is found to be a result of the erroneous water well data for
which the SWL information was recorded as water levels above sea level instead of depth to
water table (i.e. possibly a processing and/or a recording error).

81

Outlier analysis.

Elimination of the bad data in Wellogic was discussed before (p.

73). Errors in the calculation of the water level elevations are reduced by introducing the rasterbased USGS DEM. However, it is a difficult task to detect the water wells with incorrect and
anomalously high or low static water levels and these wells can sometimes be a major source of
error. In this research, for the model evaluations, a systematic outlier analysis isn’t conducted so
that the effectiveness of the data can be judged in a better sense. On the other hand, the modelbased analysis significantly helps to identify the data points with extreme values. Some of these
points are can be identified with the NOSD parameter defined in Appendix G.

Hill-Rowley et al. (2003) states that using Wellogic data at local scales without
processing might present “confused” groundwater flow patterns which are not consistent with the
hydrogeological framework of the aquifer. Even though most of the issues related to these local
data are lost at watershed scale, it is still possible to identify some of those points by comparing
the data and models. In Figure 34, two extreme values in the bedrock wells can be easily detected
using the comparative analysis.

82

Figure 34: Identifying the outliers for the Clinton Watershed, bedrock aquifer.

The SWL values that were recorded in these wells (red box) are compared to those of the
surrounding wells (Figure 35). Since within a very short distance, SWL is not expected to
fluctuate anomalously, these two data points can be considered as the “bad data”.

83

SWL5=30 ft
Bedrock Well

Stream

SWL3=20 ft

SWL2=280 ft

SWL1=380 ft

SWL4=30 ft

SWLmean=36 ft

Figure 35: Analysis of the outliers based on the model results in the Clinton Watershed, bedrock
aquifer.

84

Raster-based bedrock surface map.

The accuracy of the bedrock surface map is

important because this parameter defines the saturated thickness (water table minus bedrock
surface) of the glacial aquifers. Raster-based bedrock elevation map was generated by
interpolating the bedrock surface elevation data recorded in the oil, gas, and also deep water
wells (GWIM, 2006). However, the data resolution for these wells (with bedrock surface
elevation data) is not very good and especially for the Upper Peninsula, data is sparse.

In order to test the accuracy of the raster-based bedrock surface map (with 500 m of
resolution), following steps are applied:



Using IGW, 90 m SRTM land surface data for a portion of the Upper Peninsula, is
remapped onto a grid with a cell size of 500 m resolution



The resulting DEM grid is extracted into ArcGIS



The raster-based bedrock surface data is subtracted from the raster-based DEM data



The final product gives the glacial thickness

The reason that a statewide glacial thickness map with 500 m of resolution isn’t
generated is because it is computationally impossible to map the DEM data at state scale by
using such small grids. Glacial thickness map created for the accuracy test, along with the wells
that are used to produce the bedrock surface map, is shown in Figure 36.

85

Error (m)
-480

Data Point

-4
Scale is non-linear

Figure 36: Error locations for the bedrock surface map of the Upper Peninsula.

Figure 37 zooms into the red box shown in Figure 36 to have a closer look at the data
resolution used in the interpolation. It is observed that, at certain locations where there is no data
available, error produced with the interpolation increases, as expected.

86

Data point

Figure 37: Magnifying the error in the bedrock surface map and available data points in the
Upper Peninsula.

The cause of the interpolation artifacts shown in Figure 37 is also illustrated in Figure 38.
While the interpolation uses the data at the well locations, if the data density in between those
data points is not good enough, the interpolation might end up with very bad representations of
the reality.

87

Well with
bedrock surface data

Interpolated bedrock surface

Bedrock surface
in the field

Figure 38: Conceptual picture of the possible reason for the error in the bedrock surface map.

Since in the Lower Peninsula, data resolution is much better than that of the Upper
Peninsula, it is assumed that the raster bedrock surface map for the Lower Peninsula is a good
representation of the reality.

Results from the analysis on the accuracy of the bedrock surface map indicate that there
is a need for the revision of this map. However, for the groundwater models developed for this
research, since there isn’t a dataset available with a better resolution to generate a new bedrock
surface map, existing map is used to define the bottom elevation of the glacial aquifer.

Reproduced GIS Files

To improve the computational efficiency, to develop seamless raster-based data
definitions, and where possible, to improve the conceptual representation of the model

88

parameters, some of the statewide data files are reproduced. Final database is coupled with
models to provide a flexible and interactive data environment that allows the simultaneous use of
the vector-based and the raster-based geospatial data.

Raster-based land surface elevation.

The USGS and Earth Resources Observation

and Science (EROS) Data Center (EDC) provides access to The National Map Seamless Server
(available online at http://seamless.usgs.gov/) which makes it possible to explore and retrieve
countrywide land surface elevation data from the National Elevation Dataset (NED). NED is the
primary elevation product of the USGS and it is a seamless dataset with the best available raster
elevation data of the United States. The NED data is distributed in geographic coordinates in
units of decimal degrees, and in conformance with the North American Datum of 1983 (NAD
83). All elevation values are in meters and the NED of Michigan is available for 30 m (1 arcsecond) of resolution everywhere in the state. NED provides data with a resolution of 10 meters
(1/3 arc-second) or higher for only 85% of the United States, and most of Michigan is among the
15% where the 10 m DEM coverage was generated from oversampling of 30-meter DEM source
data (USGS, 2006; 2009).

Other than the NED data, the USGS also provides Shuttle Radar Topography Mission
(SRTM) data for the land surface elevation. SRTM is a joint project between the National
Aeronautics and Space Administration (NASA) and National Geospatial-Intelligence Agency
(NGA) to map the world in three dimensions. SRTM generates a digital topographic map of the
Earth’s land surface with data points spaced every 3 arc-second (about 90 meters) for global
coverage of latitude and longitude (USGS, 2008).

89

Even though the DEM data provided by the USGS is pre-processed, it is still necessary to
arrange the data in a way that makes it possible to use it for the groundwater models developed
in this research. Following steps are applied to generate the GIS files for the model integration:



The USGS raster-based NED and SRTM data for 10 m (sampled from 30 m DEM), 30 m,
and 90 m resolutions are downloaded (from the National Map Seamless Server) in 485,
42, and 10 geographic tiles, respectively



The numbers of the geographic tiles are optimized (e.g. each tile for 10 m dataset stores
50 MB of raster data) based on the limitations of the computational resources that are
used to extract the data for the models



A statewide polygon shapefile, that allows to locate the raster data for data transfer
depending on the geographic extensions of the model area, is created

One of the final products, statewide 90 m DEM is shown in Figure 39.

90

Land Surface Elevation (m)
584

174
Scale is non-linear

Figure 39: USGS 90 m raster-based DEM.

Raster-based hydraulic conductivity field.

Along with the SWL data, hydraulic

conductivity data is the most critical component needed for groundwater investigations. The
distribution of water levels and hydraulic conductivity dictates the groundwater speed and
direction. Both of these data are needed for most hydrogeological analyses such as contaminant
transport predictions, remediation design, and well-head delineation (Simard, 2007).
91

Statewide raster-based hydraulic conductivity maps for HCOND1 and HCOND2
parameters are produced through the following steps:



The state map of Michigan is divided into 8 seamless rectangular geographic tiles in
ArcGIS



The water wells (totally almost 200,000) for each tile are extracted in ArcGIS



The HCOND1 and HCOND2 data in the water wells are interpolated (with a number of
300 nearest data points and a cell resolution of 540 m) within the tiles by using IGW’s
Inverse Distance Weighting interpolation interface



After the conductivity fields for each tile are interpolated, they are extracted as the rasterbased maps in IGW



The resulting maps are merged in ArcGIS

Final products for the raster-based maps of the HCOND1 and HCOND2 parameters are
shown in Figure 40.

92

Figure 40: Raster-based hydraulic conductivity maps interpolated from the data.

93

0.39
Scale is non-linear

K (m/d)
88.5

a) HCOND1

0.24
Scale is non-linear

K (m/d)
80.1

b) HCOND2

Raster-based transmissivity.

The statewide transmissivity map is improved through the

following steps:



The USGS point-based aquifer test analysis database is updated with the current available
data which makes the use of 1,423 statewide data points where the previous map
(GWIM) was generated based on 233 data points



The polygons (USGS shapefiles), that were used to map the extent of the bedrock units
based on the locations of the aquifer test wells and the interpretations of the county-based
hydrogeological reports, are maintained from the MDEQ [personal communication with
Richard Mandle from the MDEQ]

Even though the number of the data points used in the production of the transmissivity
map is still not too many, it is a lot more than that of the ones provided with the GWIM, and the
new map is considered as the best that can be created with the current available information.

The locations of the data points along with the bedrock geology polygons are shown in
Figure 41. Here, the color map gives a general qualitative idea about the transmissivity
distribution of the bedrock layers (Table 14).

94

Transmissivity
Very low
Low
Intermediate
High
No data / Coldwater Shale

Figure 41: Color map for the transmissivity ranges generated based on the aquifer test data.

95

Table 14: Ranges for the color transmissivity map.
2

Class

Average T (ft /d)

Very low

< 500

Low

715

Intermediate

2,344

High

13,802

2

The minimum transmissivity value recorded in the database is 0.32 ft /d and located in
2

the western Upper Peninsula, where the maximum value is 142,062 ft /d and located in the
2

Lower Peninsula (Marshall Aquifer). For the Coldwater Shale unit, a constant value of 1 ft /d
was assigned to define the transmissivity used in the models.

After the data points are arranged based on the geological layers they represent, these
data points are interpolated through the Kriging interpolation interface of IGW within the
polygon coverage they are located in. The resulting transmissivity contour maps are then
converted to the raster maps to produce the final statewide transmissivity database.
rd

Raster-based stream data. In the NHD GIS files, for the streams of 3 and greater
orders, some polyline segments that follow the shorelines of the Great Lakes are detected. These
segments are called the coastline reaches and they are artificial paths that were added to
delineate the coasts (here the Great Lakes) of the country (USGS, 2000). Even though these
artificial segments don’t have a conceptual meaning in the groundwater models developed in this
research, in the database,stream stages are also stored for these polyline features (Figure 42). In
96

th

th

the figure, only the 5 and 6 order streams with artificial lines on the boundaries of the Great
Lakes are shown.

th

5

th

and 6 order streams

Figure 42: High order streams with artificial lines along the shores of the Great Lakes.

IGW still recognizes these artificial segments as streams in a way that might introduce
artificial fluxes in the models. To resolve this issue, artificial paths along the shorelines of the

97

Great Lakes (Michigan’s state boundaries) are removed, and a new stream data set for the stream
network is generated.

In the raster-based DEM maps, land surface elevation in a stream segment of the NHD
can be considered as the stream stage at that location. Since the raster-based USGS DEM data
provided with this research improves the accuracy of the land surface elevation data, it can also
quantitatively improve the stream stage information across the state. Within this context, Prof.
Shu-Guang Li et al. at MSU developed a code to directly extract the stage information based on
the NHD’s polyline stream network coverage projected onto the raster-based DEM. According to
this code, if the polyline passes through a raster cell, the code assigns that cell as a stream cell
and records the DEM value stored in the raster cell’s center as the stream stage for that cell
(Figure 43).

98

Non-stream cell

Stream cell

Stream cell

Stream cell

Stream cell
Stream cell

Non-stream cell

Non-stream cell

Stream cell

Figure 43: Conceptualization of the stream cells in the models.

DEM-based stream stage data is extracted, for all stream orders, from 30 m and 90 m
USGS DEM. Currently, a stream stage data set based on the 10 m USGS DEM (sampled from 30
m DEM) is also available for only the streams with the orders of equal to and greater than 3.

Great Lakes boundary.

A polyline coverage is created to delineate the shorelines of

the Great Lakes. Long-term average lake levels that were provided by NOAA-Great Lakes
Environmental Research Laboratory (website available at http://www.glerl.noaa.gov/) are
assigned to the poylines along the shorelines (Table 15).

99

Table 15: Prescribed head values assigned for the shorelines of the Great Lakes.
Lake Name

Water Level (ft)

Erie

571

Huron-Michigan

578

St. Clair

574

Superior

600

100

Appendix D: Analysis of the Hydraulic Connection between the Great Lakes and RegionalScale Heads

The primary driving forces that change the Great Lakes’ levels are precipitation and
evaporation. Water stored in the Great Lakes can change depending on the amount of the
precipitation falling on the lakes and the runoff contributing from their surrounding watersheds.
Based on the change in the water stored, the Great Lakes have had several major lake level (high
or low) periods in the past. These temporary periods can be divided into six: 1960 – 1970, 1970 –
1976, 1976 – 1983, 1983 – 1989, 1989 – 1998, and 1998 – 2006 (GWIM, 2006; NOAA, 2010).
To observe the temporal variations in the static water levels in connection with the lake level
periods, glacial wells in the Van Buren County next to the Lake Michigan and the Kalamazoo
County (Figure 44) next to the Van Buren County were separated into those periods, using IGW.
Then the average SWL values for the corresponding periods were calculated for these wells and
the results were plotted in Figure 45. Same water wells were also divided into annual quarters
(January – March, April – June, July – September, October – December) for all times to observe
the seasonal changes in the static water levels (Figure 46). The numbers of the wells for the
periods used in the both analyses and the average static water levels from the analyses are shown
in Table 16 and Table 17, respectively. Temporal comparisons show that while between the
annual quarters, average SWL changes within a few feet of range (Figure 45), fluctuations within
the Great Lakes level periods is a lot more (especially for the Van Buren County which is closer
to the Lake Michigan) than that of the seasonal changes (Figure 46).

101

KALAMAZOO
VAN BUREN
Figure 44: Location of the test area for temporal SWL analysis.

102

Table 16: Temporal analysis of the SWL based on the Great Lakes' water level periods.
Lake Level Period

Number of the Glacial Wells

Average SWL (ft)

Kalamazoo

Van Buren

Kalamazoo

Van Buren

1960 – 1970

266

283

843.50

723.48

1970 – 1976

531

914

844.77

719.04

1976 – 1983

1,320

808

846.97

710.74

1983 – 1989

2,094

901

849.29

729.57

1989 – 1998

4,088

83

847.41

753.72

1998 – 2006

3,036

2,206

842.26

725.49

Table 17: Temporal analysis of the SWL based on the annual quarters.
Annual Quarter

Number of the Glacial Wells

Average SWL (ft)

Kalamazoo

Van Buren

January – March

1,973

744

847.26

724.96

April – June

3,055

1,565

846.31

726.51

July – September

3,458

1,611

845.83

722.11

October – December

2,477

1,255

845.52

718.99

103

Kalamazoo Van Buren

Kalamazoo County, Average SWL (ft)

847.25

846.75

846.25

845.75

845.25
JAN-MAR

APR-JUN

JUL-SEP

OCT-DEC

Annual Quarters

Van Buren County, Average SWL (ft)

727
726
725
724
723
722
721
720
719
718
JAN-MAR

APR-JUN

JUL-SEP

OCT-DEC

Annual Quarters
Figure 45: Temporal SWL analysis for annual quarters in the Kalamazoo and Van Buren County.

104

Kalamazoo County, Average SWL (ft)

850
849
848
847
846
845
844
843
842
1960-1970 1970-1976 1976-1983 1983-1989 1989-1998 1998-2006
The Great Lakes Level Periods

Van Buren County, Average SWL (ft)

759
754
749
744
739
734
729
724
719
714
709
1960-1970 1970-1976 1976-1983 1983-1989 1989-1998 1998-2006
The Great Lakes Level Periods
Figure 46: Temporal analysis for the Great Lakes water level periods in the Kalamazoo and Van
Buren County.
105

The fluctuations in the groundwater levels are consistent with the high and low lake level
trends in the corresponding periods. The results shown in Figure 46 are not surprising since
groundwater and the Great Lakes are interconnected components of a regional hydrological
system. Effects of the periodic climate changes in this system are observed in the Great Lakes
and groundwater as water level changes. The connection between the Great Lakes’ levels and the
groundwater levels in both the glacial and bedrock aquifers along the shorelines of the Great
Lakes can also be observed through the SWL values recorded in the water wells along the
shores. The strong connection between the regional-scale groundwater flow dynamics and the
Great Lakes maintains the main argument for the assumption of modeling the Great Lakes’
shorelines as prescribed-head boundary condition.

106

Appendix E: Analysis of the Simulated Heads

2

Data density=0.71/km

Mean of residuals
1 standard deviation
2 standard deviations
3 standard deviations
Figure 47: Analysis of the simulated heads in the Au Gres-Rifle Watershed, glacial aquifer.
107

2

Data density=0.22/km

Mean of residuals
1 standard deviation
2 standard deviations
3 standard deviations
Figure 48: Analysis of the simulated heads in the Au Gres – Rifle Watershed, bedrock aquifer.

108

2

Data density=0.59/km

Mean of residuals
1 standard deviation
2 standard deviations
3 standard deviations
Figure 49: Analysis of the simulated heads in the Au Sable Watershed, glacial aquifer.

109

2

Data density=0.61/km

Mean of residuals
1 standard deviation
2 standard deviations
3 standard deviations

Figure 50: Analysis of the simulated heads in the Birch/St. Clair watersheds, glacial aquifer.

110

2

Data density=0.27/km

Mean of residuals
1 standard deviation
2 standard deviations
3 standard deviations

Figure 51: Analysis of the simulated heads in the Birch/St. Clair watersheds, bedrock aquifer.

111

2

Data density=0.24/km

Mean of residuals
1 standard deviation
2 standard deviations
3 standard deviations

Figure 52: Analysis of the simulated heads in the Black Watershed, glacial aquifer.

112

2

Data density=0.10/km

Mean of residuals
1 standard deviation
2 standard deviations
3 standard deviations

Figure 53: Analysis of the simulated heads in the Black Watershed, bedrock aquifer.

113

2

Data density=0.05/km

Mean of residuals
1 standard deviation
2 standard deviations
3 standard deviations

Figure 54: Analysis of the simulated heads in the Cass Watershed, glacial aquifer.

114

2

Data density=0.93/km

Mean of residuals
1 standard deviation
2 standard deviations
3 standard deviations

Figure 55: Analysis of the simulated heads in the Cass Watershed, bedrock aquifer.

115

2

Data density=0.59/km

Mean of residuals
1 standard deviation
2 standard deviations
3 standard deviations

Figure 56: Analysis of the simulated heads in the Cheboygan Watershed, glacial aquifer.

116

2

Data density=0.03/km

Mean of residuals
1 standard deviation
2 standard deviations
3 standard deviations

Figure 57: Analysis of the simulated heads in the Cheboygan Watershed, bedrock aquifer.

117

2

Data density=6.33/km

Mean of residuals
1 standard deviation
2 standard deviations
3 standard deviations

Figure 58: Analysis of the simulated heads in the Clinton Watershed, glacial aquifer.

118

2

Data density=0.09/km

Mean of residuals
1 standard deviation
2 standard deviations
3 standard deviations

Figure 59: Analysis of the simulated heads in the Clinton Watershed, bedrock aquifer.

119

2

Data density=1.03/km

Mean of residuals
1 standard deviation
2 standard deviations
3 standard deviations

Figure 60: Analysis of the simulated heads in the Flint Watershed, glacial aquifer.

120

2

Data density=2.79/km

Mean of residuals
1 standard deviation
2 standard deviations
3 standard deviations

Figure 61: Analysis of the simulated heads in the Flint Watershed, bedrock aquifer.

121

2

Data density=8.19/km

Mean of residuals
1 standard deviation
2 standard deviations
3 standard deviations

Figure 62: Analysis of the simulated heads in the Huron Watershed, glacial aquifer.

122

2

Data density=1.04/km

Mean of residuals
1 standard deviation
2 standard deviations
3 standard deviations

Figure 63: Analysis of the simulated heads in the Huron Watershed, bedrock aquifer.

123

2

Data density=2.44/km

Mean of residuals
1 standard deviation
2 standard deviations
3 standard deviations

Figure 64: Analysis of the simulated heads in the Kalamazoo Watershed, glacial aquifer.

124

2

Data density=0.55/km

Mean of residuals
1 standard deviation
2 standard deviations
3 standard deviations

Figure 65: Analysis of the simulated heads in the Kalamazoo Watershed, bedrock aquifer.

125

2

Data density=0.48/km

Mean of residuals
1 standard deviation
2 standard deviations
3 standard deviations

Figure 66: Analysis of the simulated heads in the Kawkawlin Watershed, glacial aquifer.

126

2

Data density=1.34/km

Mean of residuals
1 standard deviation
2 standard deviations
3 standard deviations

Figure 67: Analysis of the simulated heads in the Kawkawlin Watershed, bedrock aquifer.

127

2

Data density=2.73/km

Mean of residuals
1 standard deviation
2 standard deviations
3 standard deviations

Figure 68: Analysis of the simulated heads in the Lower Grand River Watershed, glacial aquifer.

128

2

Data density=0.25/km

Mean of residuals
1 standard deviation
2 standard deviations
3 standard deviations

Figure 69: Analysis of the simulated heads in the Lower Grand River Watershed, bedrock
aquifer.
129

2

Data density=1.06/km

Mean of residuals
1 standard deviation
2 standard deviations
3 standard deviations

Figure 70: Analysis of the simulated heads in the Manistee Watershed, glacial aquifer.

130

2

Data density=0.49/km

Mean of residuals
1 standard deviation
2 standard deviations
3 standard deviations

Figure 71: Analysis of the simulated heads in the Maple Watershed, glacial aquifer.

131

2

Data density=0.46/km

Mean of residuals
1 standard deviation
2 standard deviations
3 standard deviations

Figure 72: Analysis of the simulated heads in the Maple Watershed, bedrock aquifer.

132

2

Data density=2.13/km

Mean of residuals
1 standard deviation
2 standard deviations
3 standard deviations

Figure 73: Analysis of the simulated heads in the Muskegon Watershed, glacial aquifer.

133

2

Data density=1.97/km

Mean of residuals
1 standard deviation
2 standard deviations
3 standard deviations

Figure 74: Analysis of the simulated heads in the Pere-Marquette Watershed, glacial aquifer.

134

2

Data density=0.58/km

Mean of residuals
1 standard deviation
2 standard deviations
3 standard deviations

Figure 75: Analysis of the simulated heads in the Pigeon Watershed, bedrock aquifer.

135

2

Data density=1.61/km

Mean of residuals
1 standard deviation
2 standard deviations
3 standard deviations

Figure 76: Analysis of the simulated heads in the Pine Watershed, glacial aquifer.

136

2

Data density=1.67/km

Mean of residuals
1 standard deviation
2 standard deviations
3 standard deviations

Figure 77: Analysis of the simulated heads in the Raisin Watershed, glacial aquifer.

137

2

Data density=0.66/km

Mean of residuals
1 standard deviation
2 standard deviations
3 standard deviations

Figure 78: Analysis of the simulated heads in the Raisin Watershed, bedrock aquifer.

138

2

Data density=3.23/km

Mean of residuals
1 standard deviation
2 standard deviations
3 standard deviations

Figure 79: Analysis of the simulated heads in the Shiawassee Watershed, glacial aquifer.
139

2

Data density=0.88/km

Mean of residuals
1 standard deviation
2 standard deviations
3 standard deviations

Figure 80: Analysis of the simulated heads in the Shiawassee Watershed, bedrock aquifer.

140

2

Data density=2.50/km

Mean of residuals
1 standard deviation
2 standard deviations
3 standard deviations

Figure 81: Analysis of the simulated heads in the St. Joseph Watershed, glacial aquifer.

141

2

Data density=0.16/km

Mean of residuals
1 standard deviation
2 standard deviations
3 standard deviations

Figure 82: Analysis of the simulated heads in the St. Joseph Watershed, bedrock aquifer.

142

2

Data density=2.28/km

Mean of residuals
1 standard deviation
2 standard deviations
3 standard deviations

Figure 83: Analysis of the simulated heads in the Thornapple Watershed, glacial aquifer.

143

2

Data density=2.28/km

Mean of residuals
1 standard deviation
2 standard deviations
3 standard deviations

Figure 84: Analysis of the simulated heads in the Thornapple Watershed, bedrock aquifer.

144

2

Data density=0.32/km

Mean of residuals
1 standard deviation
2 standard deviations
3 standard deviations

Figure 85: Analysis of the simulated heads in the Thunder Bay Watershed, glacial aquifer.
145

2

Data density=0.07/km

Mean of residuals
1 standard deviation
2 standard deviations
3 standard deviations

Figure 86: Analysis of the simulated heads in the Thunder Bay Watershed, bedrock aquifer.

146

2

Data density=1.56/km

Mean of residuals
1 standard deviation
2 standard deviations
3 standard deviations

Figure 87: Analysis of the simulated heads in the Tittabawassee Watershed, glacial aquifer.

147

2

Data density=0.06/km

Mean of residuals
1 standard deviation
2 standard deviations
3 standard deviations

Figure 88: Analysis of the simulated heads in the Tittabawassee Watershed, bedrock aquifer.

148

2

Data density=0.45/km

Mean of residuals
1 standard deviation
2 standard deviations
3 standard deviations

Figure 89: Analysis of the simulated heads in the Upper Grand River Watershed, glacial aquifer.

149

2

Data density=3.90/km

Mean of residuals
1 standard deviation
2 standard deviations
3 standard deviations

Figure 90: Analysis of the simulated heads in the Upper Grand River Watershed, bedrock
aquifer.

150

Appendix F: Analysis of Pumping in the Southern Lower Peninsula

The SWL data that are used for the head comparisons conducted in this research
represent the long-term averages of the groundwater levels in the watersheds. These water levels
are recorded after the water wells are drilled and before the pumping in those wells starts.
However, all of the water wells in the watersheds are not drilled at the same time and the
construction dates change within a wide range between 1960s (sometimes even at earlier dates)
and 2006.

Groundwater levels in the aquifers vary naturally, both seasonally and from year to year,
over a range of several feet in most places (GWIM, 2006). Also, the water level that had been
recorded in a well in 1970 doesn’t reflect the changes due to pumping in a neighbor well that was
drilled in 1980. In other words, transient impacts of the human induced stresses are not
superimposed onto the data.

The cone of depression that has evolved over time within the Lansing area is a good
example to demonstrate the temporal variability in the water levels. Bedrock aquifer in the
Lansing area is the primary source of water for domestic, commercial, and industrial use
(Holtschlag et al., 1996). While in the 1970s, drawdown in the area reached its maximum levels
because of the significant pumping, after late 1980s, locations of the pumping wells were
changed to distribute the stress to a larger area and recover the water levels (Holtschlag et al.,
1996). To demonstrate the significant water level changes in the area, SWL data recorded in the
bedrock wells that were drilled in two different periods (between 1970 and 1980, and between
1990 and 2000) is interpolated and shown in (Figure 91).

151

SWL (m)
269

215
Scale is non-linear
a) 1970-1980 period

SWL (m)
267

233
Scale is non-linear
b) 1990-2000 period

Bedrock well
Stream
Lansing area
Figure 91: Interpolated SWL in the city of Lansing.

152

Revisiting the model of the Upper Grand River Watershed, one can clearly see the effect
(strong local bias) of the pumping in the bedrock aquifer in the Lansing area (Figure 92).

Figure 92: Revisiting the model results for the Upper Grand River. Circled area indicate the
pumping effect.

Since the water wells are not simulated in any of the models, simulated heads cannot
capture the impacts of pumping. However, pumping effect (drawdown) can be identified from
the comparisons, especially when the total pumping discharge is significant over a long-term
period. Furthermore, the bias shown in Figure 92 can be recovered if the pumping wells are
enabled in the conceptual model. Even though there is not enough evidence to say that it is also a
result of pumping, a similar bias for which the model systematically and at a local scale
overestimates the data is also detected in the head comparisons for the bedrock aquifers in the
Flint and Shiawassee Watersheds (Figure 93).

153

Figure 93: Identifying biases similar to that of the Upper Grand River case. Circled areas might
potentially indicate the pumping effect. Above: Flint Watershed, Below: Shiawassee Watershed.

154

Appendix G: Definition of the Additional Model Statistics

Additional statistics that help to quantify the comparison between the simulated and
measured heads are summarized in Table 18. A MATLAB script (Appendix I) is provided to
generate the statistics in this table and to plot the graphics for the head comparisons illustrated in
this dissertation.

Table 18: Statistics that are used to quantify the head comparisons.
Parameter (unit)
Definition
Residual Min (m)

Minimum value of residuals

Residual Max (m)

Maximum value of residuals

µ (m)

Mean value of residuals

σ (m)

One standard deviation of residuals around µ

<0

Number of negative residuals

>0

Number of positive residuals

Data Min (m)

Minimum measured head

Data Max (m)

Maximum measured head

Model Min (m)

Minimum simulated head

Model Max (m)

Maximum simulated head

RMSE (m)

Root mean square error

NOSD

Number of outliers

NRMSE (m)

Normalized root mean square error

R

2

NS

Coefficient of determination
Nash-Sutcliffe model efficiency coefficient

Residuals used in the analysis represent the true error between the measured and
simulated heads (i.e. measurement minus prediction). The NOSD in Table 18 represents the

155

number of the outliers for the condition of

absolute  residuals   absolute  mean  residuals   3  residuals   .

Root mean square error (RMSE) is frequently used in the groundwater modeling
applications to measure the difference between the simulated and measured heads at the
locations of the measurements. IGW’s GIS interface allows the automatic extraction of the
measured head values (above sea level) by subtracting the depth to water level recorded in
Wellogic from the land surface elevation value extracted from the USGS-10 m DEM (sampled
from the USGS-30 m DEM) at the well’s location. RMSE is defined by the following formula
[Equation (15)]:

n

  hs  hm 

RMSE  i 1

2

15 

n

where n is the total number of measurements in a watershed, hs is the simulated groundwater
level, and hm is the measured groundwater level.

NRMSE is the normalized value of RMSE which is calculated by dividing RMSE to the
data range [Equation (16)]:

NRMSE 

RMSE
DR

16 

where DR is the data range.

156

The smaller the NRMSE for a model is the better the model predicts the reality.

Another measure of the ability of the model to predict the head field is Nash-Sutcliffe model
efficiency coefficient (NS). NS is calculated with [Equation (17)]:

n

  hs  hm 

2

17 

NS  1  i 1
2
hm  hm





where n is the total number measurements in a watershed, hs is the simulated groundwater level,
hm is the measured groundwater level, and hm is the average of the measured groundwater
levels.

The values of NS can range from   to 1. An efficiency value of 1 corresponds to a
perfect match between the simulated and measured heads. An efficiency value of zero (NS=0)
indicates that the simulated heads are as accurate as the mean of the measured heads. An
efficiency value less than zero (NS<0) indicates that the measurements give a better
representation of the flow field than the model predictions.

157

Appendix H: Detailed List of the Model Statistics

Table 19: Statistics for the model residuals, glacial aquifers.
Watershed
Glacial Aquifer
Name
No. of
Residual
Data
µ
σ
<0
Points Min Max

>0

NOSD

Au Gres – Rifle

1,834

-35.4

101.5

5.7

11.3

493

1,341

56

Au Sable

3,184

-46.5

67.1

6.1

11.8

780

2,404

81

Birch – St. Clair

2,770

-25.3

32.1

-0.7

6.9

1,631

1,139

41

Black

369

-28.9

91.9

6.9

16.9

122

247

5

Cass

109

-22.7

27.4

3.9

19.6

33

76

0

Cheboygan

1,446

-33.5

102.8

16.6

21.6

258

1,188

3

Clinton

12,057

-43.9

46.9

3.9

7.1

3568

8,489

116

Flint

3,611

-25.6

51.7

7.1

10.2

886

2,725

12

Huron

19,938

-41.3

56.4

3.8

6.2

4,694

15,244

212

Kalamazoo

12,730

-50.1

55.0

3.4

8.1

4563

8,167

138

Kawkawlin

623

-23.0

20.1

-3.8

6.5

511

112

12

L. Grand River

14,162

-40.2

57.4

8.2

9.9

2,634

11,528

127

Manistee

5,222

-40.4

1,283

20.4

19.0

484

4,738

39

Maple

1,173

-25.8

31.6

3.3

6.1

338

835

18

Muskegon

15,054

-35.6

117.4

11.5

12.6

1,844

13,210

118

Pere – Marquette

10,227

-39.3

91.9

12.0

12.4

1,098

9,129

146

Pine

4,305

-47.6

59.1

5.3

7.2

795

3,510

41

Raisin

4,403

-28.5

41.0

2.4

6.4

1,495

2,908

38

Shiawassee

10,307

52.3

55.1

3.2

6.5

3,434

6,873

80

St. Joseph

15,448

-34.3

46.4

2.7

6.3

5,703

9,745

263

Thornapple

5,034

-34.3

55.5

8.0

8.1

774

4,260

17

Thunder Bay

1,035

-48.5

86.7

12.7

17.0

223

812

7

Tittabawassee

5,744

-45.7

107.6

17.6

21.9

1,038

4,706

3

U. Grand River

2,032

-15.2

34.6

2.8

5.5

648

1,384

16

158

Table 20: Statistics for the model performances, glacial aquifers.
Watershed
Glacial Aquifer
Name
Head
Data

Model
2

Min

Max

Min

Max

RMSE

NRMSE

R

NS

Au Gres – Rifle

174

406

175

308

12.66

0.05

0.89

0.85

Au Sable

146

414

175

391

13.29

0.05

0.96

0.95

Birch – St. Clair

165

304

175

275

6.98

0.05

0.92

0.92

Black

170

398

180

320

18.20

0.08

0.91

0.85

Cass

178

261

180

253

8.72

0.11

0.78

0.72

Cheboygan

159

412

176

334

27.25

0.11

0.96

0.85

Clinton

165

350

173

313

8.09

0.04

0.98

0.97

Flint

173

343

177

300

12.38

0.07

0.95

0.90

Huron

196

335

204

306

7.24

0.05

0.83

0.76

Kalamazoo

145

350

175

310

8.76

0.04

0.90

0.88

Kawkawlin

167

244

176

225

7.57

0.10

0.80

0.71

L. Grand River

147

311

175

292

12.90

0.08

0.86

0.76

Manistee

150

444

176

378

27.85

0.09

0.91

0.80

Maple

186

288

193

264

6.96

0.07

0.89

0.85

Muskegon

161

445

175

399

17.04

0.06

0.96

0.92

Pere – Marquette

144

374

175

287

17.25

0.07

0.89

0.78

Pine

172

340

183

312

8.95

0.05

0.97

0.94

Raisin

179

335

200

316

6.85

0.04

0.94

0.93

Shiawassee

174

350

177

297

7.21

0.04

0.97

0.96

St. Joseph

155

360

176

349

6.90

0.03

0.95

0.94

Thornapple

157

297

187

272

11.42

0.08

0.86

0.73

Thunder Bay

169

375

199

316

21.26

0.10

0.82

0.67

Tittabawassee

173

390

175

286

28.11

0.13

0.89

0.64

U. Grand River

197

336

194

30.4

6.21

0.04

0.91

0.88

159

Table 21: Statistics for the model residuals, bedrock aquifers.
Watershed
Bedrock Aquifer
Name
No. of
Residual
Data
Min
Max
µ
σ
<0
Points

>0

NOSD

Au Gres – Rifle

562

-18.4

19.4

-0.1

4.6

305

257

10

Birch – St. Clair

1,234

-27.8

28.8

-1.5

6.1

860

374

29

Black

158

-64.2

15.8

-13.2

17.4

121

37

17

Cass

2,159

-37.3

24.2

-0.1

5.7

1,129

1,030

44

Cheboygan

77

-56.3

47.9

-1.9

15.5

38

39

3

Clinton

181

-115.2

33.1

-0.9

13.3

85

96

3

Flint

9,750

-48.9

38.7

-2.4

6.1

6,584

3,166

263

Huron

2,358

-26.4

25.8

1.2

5.1

926

1,432

23

Kalamazoo

2,884

-28.7

43.7

3.3

7.0

929

1,955

44

Kawkawlin

1,735

-29.0

10.1

-4.7

4.3

1,520

215

335

L. Grand River

1,298

-39.9

44.1

-3.2

8.7

908

390

26

Maple

1,106

-19.3

25.4

-0.2

5.1

564

542

15

Pigeon

1,308

-26.4

17.2

-5.2

4.7

1,179

129

277

Raisin

1,733

-36.5

24.6

-5.7

8.5

1,236

497

72

Shiawassee

2,829

-68.9

25.5

-1.5

4.9

1,818

1,011

56

St. Joseph

972

-51.3

42.2

2.7

7.1

347

625

18

Thornapple

1,508

-28.6

22.9

2.1

5.5

533

975

11

Thunder Bay

230

-54.5

22.3

-3.1

10.0

142

88

10

Tittabawassee

223

-25.9

13.0

-1.9

6.4

131

92

4

U. Grand River

17,579

-49.4

27.4

-0.1

5.8

8,201

9,378

254

160

Table 22: Statistics for the model performances, bedrock aquifers.
Watershed
Bedrock Aquifer
Name
Head
Data

Model
2

Min

Max

Min

Max

RMSE

NRMSE

R

NS

Au Gres – Rifle

172

277

175

271

4.57

0.04

0.96

0.96

Birch – St. Clair

170

295

176

272

6.25

0.05

0.93

0.94

Black

174

274

187

264

21.84

0.22

0.26

-0.30

Cass

173

275

178

253

5.71

0.05

0.87

0.87

Cheboygan

165

270

176

273

15.54

0.15

0.67

0.55

Clinton

110

322

173

312

13.26

0.06

0.88

0.87

Flint

163

329

177

300

6.55

0.04

0.93

0.91

Huron

167

308

173

290

5.25

0.04

0.95

0.94

Kalamazoo

179

351

183

310

7.74

0.04

0.95

0.94

Kawkawlin

166

229

176

222

6.38

0.10

0.88

0.61

L. Grand River

150

251

176

251

9.27

0.09

0.73

0.69

Maple

189

253

193

244

5.08

0.08

0.84

0.84

Pigeon

163

253

176

239

6.97

0.08

0.91

0.79

Raisin

168

340

187

316

10.26

0.06

0.99

0.96

Shiawassee

163

308

177

296

5.17

0.03

0.97

0.97

St. Joseph

224

364

235

348

7.62

0.05

0.93

0.91

Thornapple

179

286

188

272

5.93

0.05

0.91

0.89

Thunder Bay

175

250

181

248

10.48

0.14

0.71

0.67

Tittabawassee

178

252

180

253

6.70

0.09

0.89

0.86

U. Grand River

194

328

194

304

5.83

0.04

0.89

0.88

161

Appendix I: MATLAB Script for the Calculation of the Model Statistics

% Statistics and head comparison for the watershed models
data=xlsread(‘WATERSHED.NAME.DATA.xlsx');
model=xlsread('WATERSHED.NAME.MODEL.xlsx');
% Eliminating zeros for data
data1=data(:);
b=isnan(data1); %Detect NaN cells
c=find(b==0); %Find non-zero cells
datanew=data1(c); %Assign a new vector for non-zero cells
datanew1=datanew(:); %New data column vector
ndv=[1:1:size(datanew1(:))]; %New data vector
% Eliminating zeros for model
model1=model(:);
b1=isnan(model1); %Detect NaN cells
c1=find(b1==0); %Find non-zero cells
modelnew=model1(c1); %Assign a new vector for non-zero cells
modelnew1=modelnew(:); %New model column vector
% Residuals
res=datanew1-modelnew1; %Residuals
plot1=plot(ndv,res)
plot1=xlabel('No. of data points','FontSize',24,'FontName','Times','Fontweight','b')
plot1=ylabel('Residual (m)','FontSize',24,'FontName','Times','Fontweight','b')
plot1=title('Residuals for the Au Sable Watershed (Drift Aquifer)',...
'FontSize',24,'FontName','Times','Fontweight','b')

162

hold on
axis([1,max(ndv),min(res)-5,max(res)+5]);
hold on
mr=mean(res); %Mean of residuals
mrvec=ones(size(res)); %Mean residual vector
mrvec1=mr*mrvec;
plot2=plot(ndv,mrvec1)
set(plot2,'Color','Red','LineWidth',2.5);
hold on
% Standard deviation calculations
% 1 std
sdr1=std(res); %Standard deviation of residuals
msdr11=mr+sdr1; %Mean + 1 std
sdr1vec=ones(size(res)); %Standard deviation vector
sdr1vec1=msdr11*sdr1vec;
msdr12=mr-sdr1; %Mean - 1 std
sdr1vec2=msdr12*sdr1vec;
plot3=plot(ndv,sdr1vec1)
set(plot3,'Color','Black','LineWidth',2.5);
hold on
plot4=plot(ndv,sdr1vec2)
set(plot4,'Color','Black','LineWidth',2.5);
hold on
% 2 stds
msdr21=mr+2*sdr1; %Mean + 2 std
163

sdr2vec1=msdr21*sdr1vec;
msdr22=mr-2*sdr1 %Mean - 2 std
sdr2vec2=msdr22*sdr1vec;
plot5=plot(ndv,sdr2vec1,':')
set(plot5,'Color','Black','LineWidth',2.5);
hold on
plot6=plot(ndv,sdr2vec2,':')
set(plot6,'Color','Black','LineWidth',2.5);
hold on
% 3 stds
msdr31=mr+3*sdr1; %Mean + 3 std
sdr3vec1=msdr31*sdr1vec;
msdr32=mr-3*sdr1; %Mean - 3 std
sdr3vec2=msdr32*sdr1vec;
plot7=plot(ndv,sdr3vec1,'-.')
set(plot7,'Color','Black','LineWidth',2.5);
hold on
plot8=plot(ndv,sdr3vec2','-.')
set(plot8,'Color','Black','LineWidth',2.5);
hold off
legend([plot2,plot3,plot5,plot7],'Mean of residuals',...
'1 standard deviation','2 standard deviations','3 standard deviations')
set(legend,'FontName','Times','FontSize',24,'Fontweight','b','Location','SouthWest')
set(gca,'FontSize',20, 'FontName', 'Times')

164

% RMSE calculation
ressq=res.^2;
par=sum(ressq);
RMS=sqrt(par/length(datanew1)); %RMS
NRMS=RMS./(max(datanew1)-min(datanew1)) %Normalized RMS
% Head comparison graph
plot9=plot(datanew1,modelnew1,'*','MarkerSize',2)
plot9=xlabel('Observation (m)','FontSize',24,'FontName','Times','Fontweight','b')
plot9=ylabel('Model Prediction (m)','FontSize',24,'FontName','Times','Fontweight','b')
tchar=num2str(max(ndv));
title(({'Head Comparison for the Au Sable Watershed (Drift Aquifer)';[tchar,' data points']}),...
'FontSize',24,'FontName','Times','Fontweight','b')
set(gca,'FontSize',20, 'FontName', 'Times')
hold on
%Coefficient of determination (R-square) calculation
mo=mean(datanew1); %Mean of observations
mp=mean(modelnew1); %Mean of predictions
odiff=datanew1-mo;
pdiff=modelnew1-mp;
term1=sum(odiff.*pdiff)
term2=sqrt(sum(odiff.^2));
term3=sqrt(sum(pdiff.^2));
rsq=(term1./(term2.*term3)).^2;
% 45 degree line
if min(datanew1) < min(modelnew1) && max(datanew1) > max(modelnew1)
165

x1=[min(datanew1)-5:1:max(datanew1)+5];
y1=[min(datanew1)-5:1:max(datanew1)+5];
else
x1=[min(modelnew1)-5:1:max(modelnew1)+5];
y1=[min(modelnew1)-5:1:max(modelnew1)+5];
end
plot10=plot(x1,y1,'r','LineWidth',2.5)
if min(datanew1) < min(modelnew1) && max(datanew1) > max(modelnew1)
axis([min(datanew1)-10 max(datanew1)+10 min(datanew1)-10 max(datanew1)+10])
else
axis([min(modelnew1)-10 max(modelnew1)+10 min(modelnew1)-10 max(modelnew1)+10])
end
axis square
grid on
hold off
% Outliers
out=find(abs(res)>=abs(msdr31));
outs=size(out);
% Nash Sutcliffe model efficiency coefficient
NS=1-(sum((datanew1-modelnew1).^2)/sum((datanew1-mean(datanew1)).^2));
% Statistics summarized
RMS
NRMS
nofdata=size(datanew1)
mindata=min(datanew1)
166

minmodel=min(modelnew1)
maxdata=max(datanew1)
maxmodel=max(modelnew1)
nneg=find(res<0); %Number of negative residuals
sizeneg=size(nneg)
npos=find(res>0); %Number of positive residuals
sizepos=size(npos)
minres=min(res)
maxres=max(res)
meanres=mr
stdres=sdr1
std2res=msdr21
std3res=msdr31
outs
NS
rsq

167

LITERATURE CITED

168

Alley, W., T. E. Reilly, and O. L. Franke, 1999. Sustainability of ground water resources. U.S.
Geological Survey Circular 1186, 86 p. http://pubs.usgs.gov/circ/circ1186/.
Anderson, M. P. and W. W. Woessner, 2002. Applied ground water modeling: Simulation of
flow and advective transport. Academic Press, San Diego, USA, ISBN: 0-12-059485-4.
Barazzuoli, P., M. Bouzelboudjen, S. Cucini, L. Kiraly, O. Menicori, and N. Salleolini, 1999.
Olocenic alluvial aquifer of the River Cornia coastal plain (southern Tuscany, Italy): Database
design for ground water management. Environmental Geology 39(2): 123-143.
Barton, G., R. J. Mandle, and J. M. A. Baltusis, 1996. Predevelopment fresh heads in the
glaciofluvial, Saginaw, Marshall aquifers in the Michigan Basin. U.S. Geological Survey OpenFile Report 95-319, 21 p.
Batelaan, O, F. de Smedt, M. N. O. Valle, and W. Huybrechts, 1993. Development and
application of a ground water model integrated in the GIS GRASS. In HydroGIS 93: Application
of Geographic Information Systems in Hydrology and Water Resources. Vienna, Austria.
El-Kadi, A. I., A. A. Oloufa, A. A. Eltahan, and H. U. Malik, 1994. Use of geographic
information systems in site-specific ground water modeling. Ground Water 32(4): 617-625.
El-Kadi, Aly (Ed.), 1995. Ground water models for resources analysis and management. CRC
Press, Boca Raton, USA, ISBN: 0-56670-100-7.
Freeze, A. R. and J. A. Cherry, 1979. Groundwater. Prentice Hall, New Jersey, ISBN: 0-13365312-9, 604 p.
Gogu, R. C., G. Carabin, V. Hallet, V. Peters, and A. Dassargues, 2001. GIS-based
hydrogeological databases and ground water modeling. Hydrogeology Journal 9: 555-569.
Goodchild, M. F., 1996. The application of advanced information technology assessing
environmental impacts, SSSA Special Publications (48).
Ground water Mapping and Inventory Project (GWIM), 2006. Technical Report.
http://gwmap.rsgis.msu.edu/. [Accessed March 2010].
Haitjema, H. M. and S. M. Bruker, 2005. Are water tables a subdued replica of the topography?
Ground Water 43(6): 781-786.
Henriksen, H. J., L. Troldborg, P. Nyegaard, T. O. Sonnenborg, J. C. Refsgaard, and B. Madsen,
2003. Methodology for construction, calibration, and validation of a national hydrological model
of Denmark. Journal of Hydrology 280: 52-71.
169

Hill-Rowley, R., T. McClain, and M. Malone, 2003. Static water level mapping in East Central
Michigan. Journal of American Water Resources Association 39(1): 99-111.
Hoaglund, J. R., G. C. Huffman, and N. G. Grannemann, 2002. Michigan Basin regional ground
water flow discharge to three Great Lakes. Ground Water 40(4): 390-405.
Hojberg, A. L., J. C. Refsgaard, F. van Geer, L. F. Jorgensen, and I. Zsuffa, 2007. Use of models
to support the monitoring requirements in the water framework directive. Water Resources
Management 21: 1649-1672.
Holtschlag, D. J. and J. R. Nicholas, 1998. Indirect ground water discharge to the Great Lakes.
U.S. Geological Survey Open Report: 98-579, 25p.
Hunt, R. J., D. A. Saad, and D. M. Chapel, 2003. Numerical simulation of ground water flow in
La Crosse County, Wisconsin, and into Nearby pools of the Mississippi River. U.S. Geological
Survey Water Resources Investigations Report 03-4154, 44 p.
http://pubs.usgs.gov/wri/wri034154/.
Johnston, C. M., T. G. Dewald, T. R. Bondelid, B. B. Worstell, L. D. McKay, A. Rea, R. B.
Moore, and J. L. Goodall, 2009. Evaluation of catchment delineation methods for the mediumresolution national hydrography dataset. U.S. Geological Survey Scientific Investigations Report
2009-2533, 88 p. http://pubs.usgs.gov/sir/2009/5233/.
Krause, P., D. P. Boyle, and F. Base, 2005. Comparison of different efficiency criteria for
hydrological model assessment. Advances in Geosciences 5:89-97.
Kuniansky, E. L., F. G. Gomez, and S. T. Gonzalez, 2004. Effect of aquifer development in
irrigation practices on ground water availability in the Santa Isabel area, Puerto Rico. U.S.
Geological Survey Water Resources Investigations Report 03-4303, 65 p.
http://pubs.usgs.gov/wri/wri034303/.
Li, S. G., H. S. Liao, S. Afshari, M. Oztan, H. Abbas, and R. Mandle, 2009. A GIS-enabled
hierarchical patch dynamics paradigm for modeling complex ground water systems across
multiple spatial scales. In Modelling of Pollutants in Complex Environmental Systems (Advanced
Topics in Environmental Science). ILM Publications. Hertfordshire, UK.
Li, S. G. and Q. Liu, 2006. Interactive Ground Water (IGW). Environmental Modelling &
Software 21:417-418.
Liao, H. S., S. G. Li, and R. Mandle, 2010. GIS-enabled multiscale modeling of a complex
remediation site. In: American Water Resources Association Spring Specialty Conference:
Geographic Information Systems and Water Resources VI.

170

Lubczynski, M. W. and J. Gurwin, 2005. Integration of various data sources for transient ground
water modeling with spatio-temporally variable fluxes: Sardon study case, Spain. Journal of
Hydrology 306: 71-96.
Marsicek, M., 2002. Geology of Michigan’s Upper Peninsula.
http://www.msstate.edu/dept/geosciences/CT/TIG/WEBSITES/LOCAL/Spring2002/Michael_M
arsicek/index.html. [Accessed March 2010].
Mende, A., A. Astorga, and D. Neumann, 2007. Strategy for ground water management in
developing countries: A case study in northern Costa Rica. Journal of Hydrology 339: 109-124.
Michael, W. J., B. S. Minsker, D. Tcheng, A. J. Valocchi, and J. J. Quinn, 2005. Water
Resources Research 41: 1-14.
National Oceanic and Atmospheric Administration (NOAA), 2011. Great Lakes water levels.
http://www.glerl.noaa.gov/data/now/wlevels/levels.html. [Accessed November 2010].
Neff, B. P., A. R. Piggott, and R. A. Sheets, 2006. Estimation of shallow groundwater recharge
in the Great Lakes Basin. U.S. Geological Survey Scientific Investigations Report 2005-5284, 31
p. http://pubs.usgs.gov/sir/2005/5284/.
Refsgaard, J. C., A. L. Hojberg, I. Moller, M. Hansen, and V. Sondergaard, 2010. Ground water
modeling in integrated water resources management: Visions for 2020. Ground Water 48(5):
633-648.
Neff, B. P., S. M. Day, A. R. Piggott, and L. M. Fuller, 2005. Base flow in the Great Lakes
Basin. U.S. Geological Survey Scientific Investigations Report 2005-5217, 32p.
http://pubs.usgs.gov/sir/2005/5217/.
Rheaume, S. J., 1991.Hydrologic provinces of Michigan. U.S. Geological Survey Water
Resources Investigations Report 91-4120, 73 p.
Seaber, P. R., P. F. Kapinos, and G. L. Knapp, 1987. Hydrologic unit maps. U.S. Geological
Survey Water Supply Paper 2294, 66 p. http://pubs.usgs.gov/wsp/wsp2294/.
Simard, A., 2007. Predicting groundwater flow and transport using Michigan’s statewide
Wellogic database, Ph.D. dissertation, Michigan State University, 125 p.
Tiedeman, C. R., M. C. Hill, F. A. D’Agnese, and C. C. Faunt, 2003. Methods for using
groundwater model predictions to guide hydrogeologic data collection, with application to the
Death Valley regional groundwater flow system. Water Resources Research (39): 1-17.

171

Tim, U. S., D. Jain, and H. H. Liao, 1996. Interactive modeling of ground water vulnerability
within a geographic information system environment. Ground Water 34(4): 618-627.
USGS, 2006. U.S. Geological Survey National Elevation Dataset. http://ned.usgs.gov/.
[Accessed March 2010].
USGS, 2008. U.S. Geological Survey Shuttle Radar Topography Mission data.
http://seamless.usgs.gov/products/srtm3arc.php. [Accessed March 2010].
USGS, 2009. U.S. Geological Survey NED release notes, 10 p.
http://ned.usgs.gov/downloads/documents/NED_Release_Notes_Feb09.pdf. [Accessed March
2010].
USGS, 2010. U.S. Geological Survey National Hydrography Dataset. http://nhd.usgs.gov/.
[Accessed March 2010].
Walton, W. C., 1992. Ground water modeling utilities. Lewis Publishers, Chelsea, MI, ISBN: 087371-679-5, 656 p.
Westjohn, D. B. and T. L. Weaver, 1998. Hydrogeologic framework of the Michigan basin
regional aquifer system. U.S. Geological Survey Professional Paper 1418, 47 p.
Zheng, C. and R. Mia, 2010. IGW/DL: A digital library for teaching and learning hydrogeology
and ground water modeling. Ground Water 48(3): 339-342.

172