 

.

vo..i.- < E. 32.0va
”4;,” .».K.j:.,r.e.m. a. M
an mm.
“a”. a _

       

.0..u1

\

    

   

O
,.J-

                 
   

‘Eu
.1"

on
' «

  

‘

      

..
v

inf,

                

21...; .
ml 13.0

"a

       

.1 Jinn-”‘47. .

¢. 1...: 1.01 «I ‘91.! d

 

      

  

. ‘ . .

8 Juanausu
‘ i . ..L—.v.p..§.,kl<ﬁmﬂ.
. . n..,§wl.-._u.i. if. r

‘ H MUN—Ni. ., In»? .souml'lali I.

 

b
a?) . It...
1!

   

 
 

, 87.4».th
‘ p..
t.

w“
%
gc

 

  

 

. , :7- , , v {.31. .ullL‘ . , .
‘dmﬂdlad; 1.4L- TV» slaughﬂlﬂiﬁjvu ‘33,}.1. u‘ ‘... ; Hr; ld' -Lu.. A. h l . : .
All. I. 1all... III . Illllﬂ I .

I.

 

\\
‘..\
m%\-

(\s

LIBRARY
Michigan State
University

This is to certify that the

thesis entitled

/n/erpre£m Hu Impum‘s % [and Use m (Ma/fr Qua/M
03mg Gm (#007sz F/ou/ and 77/01/73er SIM/61790775
lh géhe Grand Wat/4756 86:3 Mix/OTS/kd

presented by

DOV/d F Bauéé

has been accepted towards fulﬁllment
of the requirements for

 

H, S degree in GfO/ngf'fa/ 80/60/7sz

ﬂat/M

Major professor

Date 7’/é9717

0-7639 MS U is an Afﬁrmative Action/Equal Opportunity Institution

 

 

PLACE IN RETURN BOX to remove this checkout from your record.
To AVOID FINES return on or before date due.
MAY BE RECALLED with earlier due date if requested.

 

DATE DUE

DATE DUE

DATE DUE

 

 

MAY 1 9 2008

 

 

 

 

 

 

 

 

 

 

 

 

 

 

2/05 wacmm.smd-p.1s

INTERPRETING THE IMPACTS OF LAND USE ON WATER QUALITY
USING GROUNDWATER FLOW AND TRANSPORT SIMULATIONS IN THE
GRAND TRAVERSE BAY WATERSHED
By

David F. Boutt

A THESIS

Submitted to
Michigan State University
in partial fulﬁllment of the requirements
for the degree of

MASTER OF SCIENCE

Department of Geological Sciences

1999

ABSTRACT

INTERPRETING THE IMPACTS OF LAND USE ON WATER QUALITY
USING GROUNDWATER FLOW AND TRANSPORT SIMULATIONS IN THE
GRAND TRAVERSE BAY WATERSHED
By

David F . Boutt

Groundwater contributions are often overlooked in approaches that attempt to
quantify the impacts of land use on surface water quality. This is surprising because
it is known that surface water chemistry at base ﬂow is representative of the
groundwater chemistry in humid regions and that vegetation and land cover effects
the distribution of recharge/evapotranspiration throughout a watershed. The Grand
Traverse Bay Watershed (GTBW), which is undergoing large land use changes, is an
excellent area to examine these relationships. A groundwater ﬂow model linked with
GIS inputs was used to demonstrate how groundwater links land use to surface water
quality in the GTBW. It was demonstrated that a regional scale groundwater flow
and solute transport model could help separate the impact of different land use related
solute sources in a watershed. This work also shows that groundwater is a key
linkage between land use and surface water quality. Groundwater ﬂow processes are
also neglected in studies that examine the inﬂuence of land use/cover changes on the
environment. It has been shown that the integration of these changes is possible and
they do have effects on the distribution of simulated groundwater heads across the

GTBW.

ACKNOWLEDGEMENTS

First and foremost I would like to thank my parents and family for allowing and
supporting myself to choose whatever career I wanted. In the two years that I have been
working on Grand Traverse project I have had help from many different people. I would
like to thank my advisor, Dave Hyndman, for not only treating me very well as a student,
but also as a good friend. He has always stood up for me and supported me in many
ways. I would not be writing this now if I had not had his guidance and patience. I
would also like to thank Bryan Pijanowski for all his help and comments on my GIS
work, without his lack of help in some cases I would have never learned GIS. Dave Long
and Sheridan Haack have been tremendous help in organizing the collection of ﬁeld data
as well as serving in an unbiased role when problems arose on the project. Grahame
Larson was key factor in allowing myself to become a better student of geology and
sorting out the glacial deposits of the GTBW. I am indebted to Sarah Woodhams for her
work with myself on the Traverse project. I would also like to thank Tom Vogel and
Duncan Sibley for being good friends and role models for the years that I have known
them. Lastly my great friend Astrid deserves lots of credit for distracting me from my

work and allowing me to have time for other things than the Traverse project.

iii

TABLE OF CONTENTS

LIST OF TABLES ................................................................................... vi
LIST OF FIGURES ................................................................................. vii
CHAPTER 1 ........................................................................................... 1
INTRODUCTION .................................................................................. 2
STUDY AREA ...................................................................................... 5
MODELING APPROACH ........................................................................ 8
Data Collection and Analysis ............................................................... 8
Groundwater Model Development ...................................................... 12
RESULTS .......................................................................................... 18
Model Calibration .......................................................................... 18
Groundwater Flux Estimates .............................................................. 20
Delineating Groundwater Flow Divides ................................................ 21
Simulating Land Use Water Quality Impact Scenarios ............................... 23
DISCUSSION ...................................................................................... 31
CHAPTER 2 .......................................................................................... 33
INTRODUCTION ................................................................................. 34
APPROACH ....................................................................................... 37
STUDY AREA .................................................................................... 42
METHODS ......................................................................................... 46
RESULTS .......................................................................................... 49
DISCUSSION ...................................................................................... 56

iv

APPENDICES ....................................................................................... 58

APPENDIX A - Geology of the Grand Traverse Bay Watershed ........................... 59
Introduction ................................................................................. 59
Bedrock Geology ........................................................................... 60
Glacial and Recent Deposits .............................................................. 69
Hydrostratigraphic Interpretation ......................................................... 82
Conclusion .................................................................................. 87

APPENDIX B — Generation and Examples of Total and Differential Sourcesheds ...... 88
Introduction ................................................................................. 88
Tools and Preprocessing of Data ......................................................... 91
Methods ..................................................................................... 94
Sourceshed Applications ............................................................... 101
Conclusions ............................................................................... 103
ARC/INFO Function Deﬁnitions ....................................................... 104

LIST OF REFERENCES ......................................................................... 106

LIST OF TABLES

1 Evapotranspiration values. Literature derived values of evaopotranspiration for
different vegetation types ........................................................................... 38

2 Annual percentages of ET and recharge. Annual percentages of ET and recharge used
in the ﬂow simulations as well as the percent of each land use found in the GTBW ...... 38

3 Conductivity values used in ﬂow model. Horizontal and Vertical Conductivities for
both layers used in the groundwater ﬂow model ................................................ 84

vi

LIST OF FIGURES

Figure 1 Site map. The Grand Traverse Bay Watershed is located in the northwestern
portion of the Lower Peninsula of Michigan, which is home to the Traverse City one of
the fastest growing cities in the

Midwest ................................................................................................ 6

Figure 2 Well logs. Map of the roughly 11,000 residential water well logs and
approximately 7000 Oil and gas well logs were used for estimating geological properties
of the watershed ....................................................................................... 9

Figure 3 Bedrock Variogram. The bedrock elevations were kriged using an exponential
variogram with a variance of 1300, and a isotropic range of 3500 meters ................... 10

Figure 4 Elevation information. The surface elevation (a) was compiled from over 48
1:24,000 30-meter Digital Elevation Models available from the USGS. The bedrock
elevation (b) was calculated from the well logs shown in ﬁgure 2 and the glacial sediment
thickness (c) map was calculated by subtracting (a) from (b) ................................. 11

Figure 5 Surﬁcial geological map. Surﬁcal geology for the region was digitized to a
one-kilometer grid from the map by F arrand and Bell
(1982) ................................................................................................ 14

Figure 6 Simulated groundwater elevations. Predicted hydraulic heads varied between
378 masl and 177 masl, which is the elevation of the Grand Traverse Bay. Boundaries
for the ﬂow model are also shown ................................................................ 15

Figure 7 Simulated versus observed plot. Modeled hydraulic heads for the region was
plotted against observed data for the USGS wells in Grand Traverse
County ............................................................................................... 19

Figure 8 All wells calibration. The simulated versus observed plot for all the water well
information follows the trend of the one to one line very nicely. These wells were not
used in the calibration process ................................................................... 19

Figure 9 Groundwater divides. Surface water divides (Shown in white) for the watershed
were not always the same as the groundwater divides for the watershed. Note here that in
the top of the ﬁgure the groundwater ﬂow vectors ( shown in black) differ from the
surface water divide ............................................................................... 22

Figure 10a Chloride spatial distribution. Spatial plots of the residual concentrations of

chloride show that most of the underestimated values lie in the regions where other
sources of chloride are found .................................................................... 24

vii

Figure 10b Simulated versus observed chloride plot. The comparison between simulated
and observed concentrations of chloride show two distinct trends that can be explained by
separate sources of chloride, road salt, chloride from septic systems, and improper

disposal of brines .................................................................................... 27

Figure 11: Spatial distribution of chloride across the watershed. The advective transport
of water from chloride source areas is responsible for the differences in the simulated
chloride concentrations at 10 and 50 years. The 10 and 50-year simulations Show a high
chloride concentration surrounding urban areas where road density is

higher ................................................................................................. 28

Figure 12a Chloride ﬂux histogram. This histogram Show the differences in chloride
ﬂux between 10 and 50 years for cells bordering Grand Traverse
Bay .................................................................................................... 29

Figure 12b Chloride ﬂux. The ﬂux of chloride into The Grand Traverse Bay shows a
signiﬁcant lag time due to transport of chloride from distant source areas .................. 29

Figure 13 Streamﬂow distribution. Cells with a lot of impervious, urban land use/cover,
surfaces will direct more water via overland ﬂow to streams possibly causing higher
erosion rates .......................................................................................... 41

Figure 14 Current and predicted land-uses in the GTBW. The land use coverage on the
left is the 1990 land use derived from aerial photos. The two land use coverages on the
right hand side, 10 year and 50 year projections, of the ﬁgure are projected future land
use coverages generated using the Land Transformation Model. All three of these
coverages were used in the simulations .......................................................... 45

Figure 15 Recharge Sensitivity. The groundwater ﬂow model is sensitive to changes in
recharge in areas that are not green on the maps above. The ﬁgure shows how small
changes in recharge will effect the groundwater table elevations of the whole

watershed ............................................................................................. 50

Figure 16 Current land use - constant recharge simulation. The head diffemce map for
the variable vs. constant recharge simulations Show that an overall increase in head
occurs but localized decreases in head are found near the Traverse City region ........... 51

Figure 17 10 year land use projection - current land use simulation. The change in head
between the 10 year projected land use and the current land use Simulation Show changes
in the gradient developing as well as lower in the quickly urbanizing areas ............... 52

Figure 18 50 year land use projection - 10 year land use projection simulation. The

difference in heads developed in the 40 years between the two simulations still shows the
inﬂuence that the Traverse City area will have on groundwater elevations ................ 53

viii

Figure 19 Bedrock Geology. The bedrock in the GTBW, which consists mainly of shale,
provides a good barrier to vertical ﬂow in the region. The Devonian age Traverse
Formation is slightly fractured by appears not to have downward gradients ............... 61

Figure 20 Shale Outcrop. The picture here is of the Mississippian Ellsworth Shale that
outcrops at quarry just east of Central Lake. Originally this deposit was thought to be a
lacustrine clay ....................................................................................... 63

Figure 21 North to south cross-section of the GTBW. A typical North-South cross
section of the GTBW is shown above. This ﬁgure demonstrates how sediments with
different characteristics associate in the watershed ............................................. 66

Figure 22 Bedrock elevation. Bedrock elevations were interpolated using a kriging
technique from Oil and Gas wells, outcrops, and the few wells that penetrated the
bedrock ................................................................................................ 67

Figure 23 Driﬁ thickness. The drift thickness map for watershed was calculated by
subtracting the elevation of the surface for each point by the bedrock elevation for the
same point. Note that some prominent valleys exist in the watershed ...................... 70

Figure 24 Surﬁcial geology of the GTBW. The surﬁcial geology of the Grand Traverse
Bay Watershed is dominated by the high topography moraines separated by the broad
Mancelona Outwash Plain. The crossection line A-A’ is found in Figure 19 .............. 72

Figure 25 Outer Port Huron Moraine. The Outer Port Huron Moraine, which formed
between 13,500 and 13,000 years BP, rises to almost 200 meters over Grand Traverse
Bay .................................................................................................... 73

Figure 26 Inner Port Huron Moraine. The Inner Port Huron Moraine formed the
Mancelona Outwash Plain, which consists of over 150 meters of sand and gravel that is
home to the Boardman River valley .............................................................. 74

Figure 27 Manistee Moraine. The Manistee Moraine is responsible for forming the
southwestern edge of the Grand Traverse Bay Watershed, where it meets the Inner Port
Huron Moraine ...................................................................................... 77

Figure 28 Drumlins. Drumlins, like the ones seen here, dominate groundwater ﬂow and
stream drainage patterns in the Northeast portion of the GTBW ............................ 78

Figure 29 Smoothed hydraulic conductivity values. The hydraulic conductivity values
for each sediment zone were smoothed using a interpolation scheme to better represent
actual sedimentologic conditions ................................................................. 85

Figure 30 Total sourcesheds. Total sourcesheds for the watershed are generated by
calculating the total areas that contributes to a speciﬁed point .............................. 89

Figure 3] Differential sourcesheds. The differential sourcesheds shown here are
developed by calculating the drainage area between selected points, shown in red, that lie
within the same drainage basin ................................................................... 90

Figure 32 MiRIS and DEM generated rivers. This ﬁgure illustrates the major differences
between the digitized rivers used in the study and the DEM delineated rivers. Note that
points located on these two different sets could be off by more than 0.5 kilometers. .....95

Figure 33 Total sourceshed calculation from differential sourcheds. The site 63 and site
65 differential sourcesheds can be combined with the differential sourceshed for site 9 to
yield a total sourceshed for site 9 ................................................................ 100

CHAPTER 1: DEVELOPMENT OF A HIGH RESOLUTION
GROUNDWATER FLOW MODEL TO DESCRIBE POTENTIAL
LAND USE IMPACTS ON WATER QUALITY IN A LARGE
WATERSHED

INTRODUCTION

Groundwater contributions are often overlooked in approaches that attempt to
quantify the impacts of land use on surface water quality. This is surprising because it is
known that surface water chemistry at base ﬂow is representative of the groundwater
chemistry in humid regions (Modica et al., 1997). Groundwater chemistry in turn
represents time-weighted averages of anthropogenic inputs from the land surface. These
spatially and temporally varying impacts of land use on surface water quality can be
evaluated using sophisticated tools such as groundwater ﬂow models and geographic
information systems (GIS). In this paper we use a groundwater ﬂow model linked with
GIS inputs to demonstrate how groundwater links land use to surface water quality in
humid regions.

Establishing and quantifying linkages between land use and water quality is best
accomplished at a regional scale because factors driving land use change operate at a
variety of scales, and regional scale groundwater ﬂow models are more robust. For
example, true hydrologic boundaries can be represented at a regional scale with a
reasonable degree of certainty. These regional scales models can also provide estimated
ﬂuxes between surface water and groundwater, delineate groundwater source areas for
lakes, rivers and bays, develop boundaries for site speciﬁc transport models, and
demonstrate the inﬂuence of land use on water quality. Unfortunately, few cases have
been examined with large-scale high-resolution models.

Numerous ﬁeld studies have demonstrated that changes in land use/cover affect

groundwater dynamics and surface water quality. Pucci and Pope ( 1995) identiﬁed

signiﬁcant differences in groundwater systems between the developed and undeveloped
areas in the northern Coastal Plain of New Jersey. They found that development affected
the redistribution of recharge and discharge areas, and reduced groundwater discharge to
streams. Thorbum et a1. (1991) observed that forest clearing had a substantial effect on
groundwater recharge with average rates of 29 to 70 mm per year in two cleared
catchments, compared with 7 mm per year in an uncleared catchment. Researchers have
also found correlation between land use patterns and stream chemistry. For example,

F lintrop et a1. (1996) showed that land use patterns affect the major ion and nutrient
geochemistry in selected tributaries of the Rhine River. These observed changes
underscore the need to study the effects of land use on stream water quality, groundwater
quality, and related hydrologic processes.

Many researchers have examined the effect of topography and land use on surface
water processes using a popular code called TOPMODEL (Beven and Kirkby, 1979),
which is based on the concept that topography exerts a dominant control on ﬂow routing
through watersheds. This code, which uses the conservation of mass equations for
reservoirs within a watershed, is primarily designed to explore surface water runoff
processes.

Numerical models have been developed to examine groundwater and streamﬂow
response to a variety of vegetation changes. For example, Hookey (1987) developed a 2-
D model to simulate delays in groundwater response to clearing a catchment.
Mcphereson and Peck (1987) simulated an increase in stream water salinity after forest
clearing. Few studies have coupled groundwater ﬂow and solute transport models to

quantify potential inﬂuences of altered landscapes on solute ﬂuxes at a regional scale.

Dunn and Mackay (1995) and Dunn et a1. (1996) were one of the ﬁrst groups to use
coupled models, SHETRAN with a solute transport code, to examine the fate of nitrates
from agricultural sources in a large watershed. Coupled regional groundwater ﬂow and
solute transport models could provide a method of predicting land use related solute
transport through watersheds.

Groundwater models are commonly used to describe regional ﬂow patterns
(Anderson and Woessner, 1992), while solute transport models are less commonly
explored at this scale. Two-dimensional basin hydrodynamic models have been used to
explore a range of scales from 40 to 1000 km across, and up to 8 km deep (Person et al.,
1994, Bethke, 1986, Garvin, 1986). On a regional scale (10,000 — 100,000 kmz), the
USGS Regional Aquifer Systems Analysis (RASA) initiative supported model
development for aquifers throughout the United States (i.e., Mandle, 1992). However,
the RASA model grid parameterization was generally too large to effectively explore
solute transport processes. One example of a large transport model was presented by
Martin and F rind (1998), who developed a high-resolution groundwater model of a 750
km2 heterogeneous aquifer system to identify capture zones of municipal wells.

Our large-scale (3300 kmz) high-resolution (100x100 In cells) model enables us to
quantitatively explore the inﬂuences of land use and land cover on water quality. Our
regional scale groundwater ﬂow model simulates hydraulic heads using real hydrologic
boundaries. We integrate land use/ land cover changes into this model using
relationships between speciﬁc land use types and their associated chemical inputs to
groundwater; and describe the potential inﬂuence of non-point source solutes on surface

water quality.

STUDY AREA

Michigan's Grand Traverse Bay Watershed (GTBW) provides a good opportunity
to explore linkages between land use and water quality. The GTBW, located in the
northwestern portion of Michigan’s Lower Peninsula (Figure 1), is an important natural
and economic resource for the State of Michigan and the Great Lakes Region due to its
great scenic beauty, excellent water quality, and recreational opportunities. However, the
watershed is currently undergoing some of the most rapid population growth and land use
change of any region in the United States (Vesterby and Heimlich, 1991). From 1990 to
1998, the watershed increased its resident population by over 10%, making it one of the
fasting growing areas in the Midwest. Traverse City, with a resident population of
approximately 100,000 and a summer recreational population often exceeding 500,000, is
the largest city in the watershed. We are currently examining the effects of these land use
changes on water quality to help manage the region's superior water quality, which is of
considerable importance to watershed stakeholders (Pijanowski et al., unpub.).

The regional geology is characterized by the Quaternary glacial advances and
retreats, in which glaciers carved deep valleys into the shale and limestone bedrock and
deposited sediment accumulations as thick as 365 meters. Sediment characteristics vary
widely across the watershed, in some areas changing from a thick lacustrine clay to a
coarse grained moraine within 100 meters. The regional geology has been characterized
in a variety of reports and papers (e. g., Leverett and Taylor, 1915; Martin, 1957;
Eschman et al., 1973). Martin (1957) found that most of the glacial sediments in the
southern portion of the watershed are coarse-grained moraine and outwash deposits.

Blewett and Winters (1995) were the ﬁrst to conduct a detailed study of the Port Huron

moraine system, which crosses the GTBW. They found that the deposits of this system
were stratiﬁed and had more characteristics of an outwash setting, despite its earlier
classiﬁcation as a moraine. The northeastern and northwestern portions of the watershed
are highly drumlinized, which controls surface water ﬂow patterns in these regions.
Devonian limestone and shale outcrop along the east Grand Traverse Bay shore, and
Mississippian shale outcrops northeast of Torch Lake (Kesling, 1974).

The Boardman river is the main tributary draining the GTBW, which is one of the
few oligotrophic bays in the Great Lakes region. This watershed contains over 100 lakes,
including Torch Lake and Elk Lake (Figure l), which are likely quite sensitive to
environmental change. The region receives an average annual rainfall of 105 cm, of
which approximately 40 cm is recharged to the water table, 50 cm is evapotranspired, and
15 cm becomes overland ﬂow to streams (Cummings et al., 1990).

Relationships between land use and water quality have been documented in the
GTBW. Rajagopal (1978) examined land use impacts on groundwater quality near Old
Mission Peninsula and Traverse City (Figure l), and found high levels of nitrate in wells
surrounded by cherry orchards. From 1984 to 1986, a USGS study examined
relationships between land use and water quality in Grand Traverse County (Cummings
et al., 1990). Thirty-four water wells and 24 surface water sites were sampled throughout
their study, which showed possible correlation between nitrogen input (fertilizer, septic,
etc.) and groundwater nitrate concentrations. Although this study demonstrated linkages
between agricultural land use and water quality in the GTBW, they did not explore the
land use impact of road salt (halite) application on stream and groundwater Chloride

concentrations.

MODELING APPROACH

Data Collection and Analysis

One of the challenges to regional scale groundwater ﬂow modeling is compiling
the necessary data for an effectively parameterized model. The approach that we take in
this paper is to parameterize our groundwater ﬂow model using a small number of zones
based on the Shallow aquifer geology. We use Geographic Information Systems (GIS)
approaches, speciﬁcally ARC/INFO (ESRI 1998a) and ArcView (ESRI 1998b), to
manage, manipulate, and analyze the spatial data for our groundwater ﬂow and solute
transport modeling, because managing and analyzing these large spatial databases is
difﬁcult with traditional modeling tools.

A large geologic database was compiled for the GTBW region from oil/gas and
water-well drilling logs and regional geologic maps. A digitized kilometer-scale surﬁcial
geology map for the region (Farrand and Bell, 1984) provided a hydrogeologic
framework for the unconﬁned aquifer. To reﬁne this framework, we used more than
11,000 residential water-well logs from the Michigan Department of Natural Resources
and approximately 7,000 oil and gas well logs from the Michigan Department of
Environmental Quality - Geological Survey Division (Figure 2). The residential water-
well logs provided information about degree of aquifer heterogeneity, and were used to
develop initial hydraulic conductivity values for delineated aquifer zones. Many of the
water well logs were incorrectly located or had inadequate and unreasonable geologic
characteristics. As a result a series of quality checks were performed on the database
before any data was incorporated into our model. Water wells were removed from the

database if the recorded elevation did not match the regional 30-meter Digital Elevation

_|—r

N Model Boundary
all uses Wells A
. Water Wells
5 N
'1 Oil and Gas Wells
0 4 8 12 Kilometers
E:

 

 

 

 

 

 

2051M

 

 

 

 

2410000

 

 

 

/
g { 1?]
z .2 x! ,r/ .
if, ' w“: ', t."
I“ ‘
’ 1" ‘ ${‘.“

 

uroooo
r
1’
, .
it“

 

 

 

 

 

 

 

 

 

Figure 2: Well logs. Map of the roughly 11,000 residential water well logs and
approximately 7000 oil and gas well logs were used for estimating geological properties
of the watershed.

Model (DEM) within +/— one-meter, which removed about 6,000 of the 11,000 residential
wells. Oil and gas well logs were also examined using this method, but only 300 of these
7,000 wells were rejected because the location and elevation values were more reliable.
In addition, well logs were removed if the geologic characteristics were unreasonable,

such as bedrock on top of unconsolidated sand and gravel.

 

Bedrock Variogram

 

 

 

(Hulk-«Kl S lull-out Inn-«u Isa-«gut 21m“); 3 six-+04 tux-«H

 

 

 

Figure 3: Bedrock Variogram. The bedrock elevations were kriged using an exponential
variogram with a variance of 1300, and a isotropic range of 3500 meters.

The elevation and lithology of the bedrock surface in the watershed was mainly
provided from the oil and gas well logs. Additional bedrock elevation values were taken
from the few water wells that were deep enough to reach bedrock and the few bedrock
outcrops in the northeastern portion of the watershed. The bedrock elevation data were
interpolated to the groundwater model grid using ordinary kriging using an exponential
variogram with a variance of 1300, and an isotropic range of 3500 meters (Figure 3). The
interpolated bedrock elevations deﬁne the bottom of the glacial aquifer for the
groundwater model. The glacial sediment thickness was calculated by subtracting the

bedrock elevation grid from the 30-meter DEM (Figure 4). Notable features in these

.2: Sea A3 ucuombnsm >n noun—:23 mm? 9:: A3 3053::

30853 332m 05 can N ocswc E 5537. awn: :03 2: 60¢ toga—:23 33 3V coming—u
x0835 3:. .mOmD of 50¢ 0322;“ 2352 aosnixm 13mm:— SuoEém ooquN;

3. $30 Soc “.2388 33 A3 c3332» oowtzm 2:. .coszcomﬁ cosmic—m "v 9;»:—

 

 

 

D 92.9" D278. D 2...»: D Enema D
D Ev. 0* on on or o 378. l 8. .8 D Ev. 0' on ON or o 3.4:. I 92.3” D
I 3. .2. I R .8 I I
I . I I I
I D I I

  
  
 

 

mmoconu. E6 _m_om_m .0 85920 xoocoon .n cozm>o_m 895m .m

 

 

 

 

maps include bedrock valleys in the southeastern portion of the watershed, and the deeply
incised valley of the Boardman River in the lower central portion of Figure 4b.

The hydrophysiographic database for the GTBW was obtained from the Michigan
Department of Natural Resources MiRIS Ofﬁce. MiRIS is a statewide land use database
developed in the late 19703 from 1224,000 aerial photography. MiRIS line ﬁles of the
transportation network and locations of rivers, lakes, and the bay boundary were
integrated into our GIS database to identify model boundaries and potential solute

sources for water quality scenarios.

Groundwater model development

Our conceptual model for the GTBW was developed using all available
hydrogeologic information. Based on the available geologic maps and well log data, we
chose to develop a two layer model that accounts for ﬂow and transport through the
glacial sediments. The computational resources for data processing a model with more
than two layers were prohibitive. The bedrock underlying the glacial sediments is mostly
low-conductivity shale, which limits vertical movement of water in or out of the shallow
aquifer, although a small amount of limestone exists in the northern areas of the
watershed (Kelly, 1968).

We developed a 3-dimensional groundwater ﬂow model of the GTBW with a
IOO-meter by IOO-meter ﬁnite difference grid cells using MODF LOW (Macdonald and
Harbaugh, 1988). The model has over 1 million active cells, with more than 34,000 river
cells (including the lakes). The Groundwater Modeling System (BYU, 1994)

preprocessor allowed us to convert GIS based data into model input ﬁles. This ability to

12

incorporate GlS coverages into the model is necessary to examine the inﬂuence of land
use and land cover on groundwater ﬂow and solute transport at these scales.

The hydraulic conductivities for the region are parameterized into zones that are
concordant with the mapped glacial units (Figure 5). Both layers have the same
horizontal conductivities except in areas where the lacustrine sand and gravel high
conductivity units overlie low conductivity clays. Initial hydraulic conductivity values
for the zones are consistent with the measured lithologies recorded in reliable well logs
(Freeze and Cherry, 1979). For example, conductivity values in the sand and gravel
outwash deposits will be higher than values in areas with till or ﬁne-grained lacustrine
deposits. Assigning a single hydraulic conductivity value for each geologically
parameterized zone is one of several possible methods to incorporate geology into
groundwater models.

The hydrologic boundaries used in our model are shown in Figure 6. Constant
head boundaries are used along the perimeter of Grand Traverse Bay, Lake Michigan,
and Lake Charlevoix, with the head values set to known lake elevations. The remaining
boundaries are signiﬁcant surface water ﬂow divides, which are represented as no-ﬂow
boundaries. The northeast and northwest boundaries of the model were extended beyond
the surface water divides because of the uncertainty in the location of these divides, as
discussed in the "Delineating Groundwater Flow Divides" section. Uncertainty in model
boundaries commonly leads to poor prediction of hydraulic heads (Anderson and
Woessner, 1992).

Effective recharge for the model was estimated by taking the total annual

precipitation (105 cm) and subtracting losses due to estimated evapotranspiration and

13

          

  

Aquifer Materials A K

   
  
  
  
 

[:1 Dune Sand 10 mlday
Lacustrine Sand and Gravel 21 mlday
Glacial Outwash 24 rnlday
[I Glacial Till- Medium 1 mlday
{8] Glacial Till - Coarse 3 mlday
El End Moraines - Coarse 19 m/day

8 16 24 Km

      

Figure 5: Surﬁcial geological map. Surﬁcal geology for the region was digitized to a
one-kilometer grid from the map by Farrand and Bell (1982).

 

f \

/\/ Contours of Smulated Hydraulic Heads

N Constant Head

/ \ / No Flow Boundary

     

Sinulatod Hydraulic Heads

[:] 171 - 227

- 227 . 278 A
- 278 - 328 N
- 328 - 379

 

 

L: No Data

 

Figure 6: Simulated groundwater elevations. Predicted hydraulic heads varied between
378 masl and 177 masl, which is the elevation of the Grand Traverse Bay. Boundaries
for the ﬂow model are also shown.

15

overland ﬂow for the region (See Study Area Section). Annual precipitation varies only
by a few centimeters across the watershed, but soil type and land use/land cover types
could result in Spatially variable recharge. For example, in areas where the soils are
predominantly sand, recharge would be higher than in areas with clay rich soils (i.e., till).
Since the necessary soil data were not readily available at the time of this model
development, we used a single recharge value of 40 cm/yr for the entire region.
Evapotranspiration was not directly simulated in the model, since it is subtracted from the
effective recharge value.

River cells for the model were developed based on topography due to problems
with the regional MiRIS river coverage. Problems with this coverage included
incorrectly coded lines, such as roads, and inconsistent river gradients obtained by
overlaying the river coverage with the DEM. River elevations increased as well as
decreased downstream which is clearly unrealistic. This oscillatory effect was likely
caused by the spatial averaging of the DEM and slight errors in the locations of the
digitized MiRIS line locations. To address these problems, ARC/INFO was used to
delineate drainage networks based on a DEM with 30x30 meter cells. First, the
FLOWDIRECTION function was used to calculate the direction water would ﬂow on the
land surface based on this high resolution DEM, then the FLOWACCUMULATION
function was used to calculate the number of cells ﬂowing into each grid cell. Cells with
a high-calculated ﬂow accumulation are areas of concentrated ﬂow, and can thus be used
to identify likely stream channel locations. A threshold of 400 cells was placed on this
process to deﬁne river reaches that approximated the extent of the MiRIS ﬁles.

Occasionally rivers were generated that did not exist, so the method was compared with

1:24000 topographic maps and aerial photographs where possible. Generated river cells
that were not in the maps and photos were deleted to better represent actual hydrologic
conditions. The modiﬁed drainage networks were then upscaled to a 100 by 100 m

model grid using the lowest 30-m cell elevation within each cell as the river elevation.

17

RESULTS

Model Calibration

The seven hydraulic conductivity values were slightly modiﬁed in the model
calibration process to minimize the squared residuals between observed and simulated
heads (Figure 7). Conductivity values from pump tests in the region (Cummings et al.,
1990) compared very well to the values in the calibrated model. Several head datasets are
available for the region, with monthly readings between January 1985 and January 1986,
and another set in August 1998. The steady state calibration, using the 1998 values,
provides a reasonable match between Simulated and observed heads (Figure 7). The
simulated heads also match the trend in the large regional database of heads measured
during well installations (Figure 8). These regional data were limited to the well
locations that matched the DEM elevations within +/- one-meter, as discussed in the Data
Collection and Analysis section, although there is still a large degree of uncertainty in
these values. Even if the measurements were accurately recorded, they were measured
during a wide range of hydrologic conditions. As a result, the measured heads will vary
across a wide range relative to the steady state simulated values. For example,
Cummings et a1 (1990) found that the hydraulic head near Fife Lake in the Grand
Traverse Bay region varied by approximately 2 meters from 1976 until 1988, with
average seasonal changes of roughly 0.5 meters.

Hydraulic heads in the GTBW vary from 177 meters at the Grand Traverse Bay
boundary to over 350 meters in the eastern high topography areas (Figure 4a). It is

apparent from Figure 5 that signiﬁcant groundwater ﬂow divides exist on both the

18

Simulated

Simulated

Simulated vs. Observed Heads for USGS Wells
in Grand Traverse County
300.| rrrlrlywvlelrlvl11.1y¥+

 

 

280 T

L .2 k_.: ., .A _

260
240 ‘
220

200

 

180 ,

 

 

1 l A 1 LL i A I J

160 180 200 220 240 260 280 300

Observed
Figure 7: Simulated versus observed plot. Modeled

hydraulic heads for the region was plotted against observed
data for the USGS wells in Grand Traverse County.

Simulated vs. Observed Heads for All Wells

350 l" T' _T rrrrrrr T T“ V ﬁT I v 1 f 1 I 1 V’ r V
l- i .
F v
A
. -<
a

 

300 E .
250 L

200 .

 

 

 

100 150 200 250 300 350

Observed

Figure 8: All wells calibration. The simulated versus
observed plot for all the water well information follows the
trend of the one to one line very nicely. These wells were
not used in the calibration process.

19

northeastern and northwestern boundaries of Grand Traverse Bay. The drumlinized
topography in the northeastern portion of the watershed results in signiﬁcant variability in
these simulated heads. The Boardman River southeast of Traverse City (Figure 1) is a
signiﬁcant groundwater discharge area due to the large topographic changes. This can be

observed in the high simulated hydraulic gradients in this portion of Figure 6.

Groundwater Flux Estimates

We estimated groundwater ﬂuxes into Grand Traverse Bay using the groundwater
model because it is very difﬁcult to obtain measurements of groundwater ﬂuxes in the
ﬁeld. Approximately 8 cubic feet per second of direct groundwater discharge to Grand
Traverse Bay was calculated based on the calibrated ﬂow model. This estimate was
calculated by summing ﬂuxes into all the cells bordering the bay (Figure 6), accounting
for approximately 7% of the model outﬂow. Since groundwater is also the primary
source of river water, it also accounts for nearly all the streamﬂow into the bay.

Previous studies (such as Brater, 1972) have estimated direct groundwater inputs
into Grand Traverse Bay using mass balance approaches. They estimated the inputs and
outputs of the system, including the discharges of streams, total precipitation, and
evapotranspiration. The difference between inputs (rainfall) and the outputs (streams and
evapotranspiration) was assumed to be the groundwater discharge to the bay. This
approach estimated that between 0 and 4 cfs of groundwater was directly entering the
Grand Traverse Bay. In contrast, our method can estimate the spatial distribution of

discharges and can predict solute ﬂuxes into lakes and rivers.

20

Delineating Groundwater F low Divides

Groundwater divides commonly mimic surface water divides in large watersheds,
thus if Signiﬁcant topography exists surface water divides can often be represented as no-
ﬂow boundaries. An alternative method of deﬁning these divides is to predict their
locations based on ﬂow models that use true hydrologic boundaries. In our model, the
northeast and northwest boundaries were extended beyond the surface water divides
because of uncertainty in these groundwater divides. Since there is little topography in
these areas and the bottom the aquifer has a signiﬁcant slope (Figure 4b), we allowed the
model to delineate these groundwater divides. Differences in these divides can yield
discrepancies in the predicted groundwater and solutes ﬂuxes if signiﬁcant solute sources
are located above regions where groundwater and overland ﬂow paths do not drain to the
same location.

The simplest technique to develop groundwater divides is to plot simulated
groundwater ﬂow vectors and digitize the boundary between regions that ﬂow into the
watershed and regions that ﬂow out of the watershed (Figure 9). An alternative is to
import the simulated hydraulic heads into ARC/INFO and calculate ﬂow directions and
ﬂow accumulations, using the approach we used to delineate our river network. This
method is Similar to delineating watershed boundaries for surface water catchments
(Homberger, 1998). Figure 7 illustrates a few of the differences between the region’s
surface water and groundwater divides, with the white line representing the surface water
divide and the ﬂow vectors indicating the groundwater divide. In this western portion of
the watershed, the surface water divide clearly deviates signiﬁcantly from the

groundwater divide.

21

.oEZv 883 custom 05 Eat
5&6 283 E 550% V E90? Bo: 6333:on 05 05mm on“ no a8 2: E 85 So:
802 .wonﬂoﬁa 05 o8 823v ooﬁavezouw 05 mm 08.8 05 983—“ «o: 23» conﬁnes?
05 o8 ASE? E :Boawv 822v 633 coatsm .moEzv 58>»

   

mug

/l
H...

..

a. ... .opﬂwwl
.. . eff/56W?

22

Simulating Land Use Water Quality Impact Scenarios

A coupled groundwater ﬂow and solute transport model allowed us to infer
potential impacts of land use on surface water quality. The groundwater ﬂow model
provides ﬂow velocities for the solute transport simulation. The scenario we chose to
examine represented the fate of chloride, which is a common solute found throughout the
GTBW. Chloride is a stable conservative anion that is commonly derived from the
dissolution of halite applied to roads in the winter months and also from urban sources
such as septic ﬁelds. An additional source of chloride in the GTBW region is the
improper disposal of subsurface brines brought to the surface by the drilling of oil and
gas wells (see Figure 10a). This scenario illustrates the potential inﬂuence of the
dominant source of chloride, although other scenarios could explore alternative solute
sources. The goal of conducting this scenario was: (1) to determine if observed chloride
concentrations in stream samples are the related to the application of halite to roads, and
(2) to demonstrate that groundwater transport of halite is likely a signiﬁcant factor that
inﬂuences the observed stream chloride concentrations.

Road salting with halite is a common practice in the Midwestern United States to
manage icy roads. Many groups (e.g., D’Itn', 1992; Jones and Sroka, 1997) have studied
the inﬂuence of halite on groundwater quality. In this simulation, we used chloride
concentration as an analog for halite, which is primarily sodium chloride but can include
small percentages of magnesium chloride or potassium chloride. Since chloride is
conservative, no retardation factors or reactions were simulated.

The simulation was run for 50 years, with a constant recharge source along salted roads

in the GTBW. Local county road commissions provided the locations of roads with

23

 

 

 

 
  
 

 

 

 

   
  
    

 

 

A Reported contamination of groundwater by bnne‘s
Surface water sampling locations
. Potential Brine/Septic Influence
A Potential Hallie Inﬂuence
I Nee served by mu'ucpel sewers
Major OII ma Gas Irend : f
Roads
D Model boundu'y N I
0 7 14 21 Kllomete rs
A
. —\_I
A (fr ,3
(-
/J
Ii
A «. X
A ‘ ‘ i
‘5 l
3
J/ - r (
A J

 

0 .
A

Figure 10a: Chloride spatial distribution. Spatial plots of the residual
concentrations of chloride Show that most underestimated values lie in the
regions where other sources chloride is found.

24

applied halite. Although road salt is only applied in the winter months, some halite is
likely retained in the unsaturated zone and distributed to the water table throughout the
year, hence we simulate the average annual input.

Measurements of solute concentrations discharged to surface water bodies were
calculated from groundwater ﬂow and solute transport model. In the simulations,
chloride concentrations were saved every 5 days for cells along the shoreline of Grand
Traverse Bay, while chloride concentrations were saved every ten years for the entire
watershed. Simulated chloride concentrations in streams were calculated every ten years
for each of our sampling sites noted on Figure 10a. At each stream node the ﬂux
(m3/day) of water discharged from the aquifer and the concentration of chloride (mg/L)
were used to determine a total mass ﬂux (mg/day) of chloride at that stream node. Both
the mass ﬂux of chloride and the ﬂux of water were summed down stream paths to
sampling locations. Simulated concentrations (mg/L) of chloride in the stream at each
sample site were then calculated by dividing the total mass ﬂux of chloride (mg/day) by
the total streamﬂow (m3/day). Since we are comparing our simulated concentrations to
observed concentrations that were taken at low-ﬂow hydrologic conditions, we are only
concerned with groundwater inputs into the stream.

Figure 10a shows the location of sample points relative to potential solute
sources. The striped horizontal lines refer to an area where there is a great risk for
groundwater contamination by the improper disposal of subsurface brines and the hollow
triangles indicates sites with documented groundwater contamination by brines
(Skillings, 1982; Fryer, 1982). The striped vertical lines refer the Traverse City region,

which has potential septic and urban chloride sources.

25

Our simulation appears to identify road salt transported through groundwater as a
signiﬁcant chloride source to the stream water in the GTBW. The only simulated
chloride source was halite application to roads. Figure 10b Shows two primary trends in
the simulated versus observed chloride concentrations that support this concept. Our
initial concentration of 100 mg/L was reduced to 70 mg/L to provide a better ﬁt to the
lower trend (plus symbols) in simulated and observed concentrations, because these were
measured in regions where there are not known oil brines or urban sources (Figure 10a).
As a result, the most likely chloride source in these areas is halite application to roads.
The upper trend (solid black circles) represents simulated concentrations that are
underestimated with respect to observed chloride concentrations in streams. These
locations are all either in, or downstream from, the noted oil/gas trend region (horizontal
stripes) or the portion of the Traverse City without municipal sewer service (vertical
stripes). The decrease in the chloride concentration from the Boardman's upper reaches to
the region just upstream of Traverse city are likely due to dilution from groundwater with
little chloride input due to the Pere Marquette State Forest, which is a region with few
salted roads.

Figure 11 shows predicted chloride concentrations for 10 and 50 years after the
initiation of halite application. Areas with high simulated chloride concentrations are
mostly within urban areas, such as Traverse City, due to the high density of salted roads.
The change from the 10-year to the 50-year concentration proﬁle (Figure l 1) shows the
expected increase in chloride concentrations across the watershed through time.

The temporally varying ﬂuxes into the Grand Traverse Bay from 10 to 50 years is

illustrated in Figure 12a as a histogram of the difference of the total chloride ﬂux to the

26

Halite Simulation

 

 

A Halite Inﬂuence

0 Potential Brine/Septic Inﬂuence

 

 

 

 

 

 

 

 

 

 

 

C
.9
4a
m
b d l ‘ﬁ—IT' Ifi IIIII
C a) I : ‘
a) l- . u
U E ' 2
CC) _ o o ..

15 - ----------------------------------------------------- -
o - . , I

o

e ' -
.5 . :
E _ 2
0 10 - ----------- o" -------------- -
L- l- . E u
G) ' ‘0‘ :
.O :0 E
0 I- O u E
‘6' ' ". ‘l i
O P A
w 5 _ .......................... A. ........................... g .............. -
*- 5 :
a _ n :
E i ‘
g 0 aaaaa A a a l 111111111 i lllllllll i j a a4
.0 o 5 10 15 mell-
0

Simulated Chloride Concentration after 50 years

Figure 10b: Simulated versus observed chloride plot. The comparison between
simulated and observed concentrations of chloride Show two distinct trends that can be
explained by a separate sources of chloride, road salt, chloride from septic systems, and
improper disposal of brines.

27

.conwE mm base can: Duos? mecca 5.55 365556 :oumbcoocono DES—no awE a 26% mucus—:56
82?? van 3 2C. .88» om can 9 a 2235:0280 DES—no 38386 2: E moococmEE 2: no.“ Bewcoamﬂ mm mecca
850m DES—no 50¢ $83 no teams“: 0289? BE. dosage? 05 388 02.830 Mo 82316:. Guam "3 2:“:—

 

28

Difference in Fluxes of Chloride between 10 years and 50 years

~r ili'i ‘7 ‘V'iri h 7* 7pm,,771‘, r ‘1

l

l
1
J
l
l
1
J

4
if
1
1
°0 1
r:
2;» J1
«3
mm 1
ea 1
Ulu 1
g‘ :> .
02 .‘
33*" l
.013 1
ES. 1
ZED

MLJ 1L... 5..-.1. LL .11... ._..A . ..

o 2000 4000 6000 8000 1 10‘
Chloride Flux Differences

 

Figure 12a: Chloride ﬂux histogram. This histogram show the differences in
chloride ﬂux between 10 and 50 years for cells bordering Grand Traverse Bay.

Flux of Chloride into Grand Traverse Bay

3500 . _ . ... .
3000 '
2500
2000
1500 :

1000 '

Chloride Flux (kg/day)

 

Time (vears)

Figure 12b: Chloride ﬂux. The ﬂux of chloride into The Grand Traverse Bay
shows a signiﬁcant lag time due to transport of chloride from distant source
areas.

29

Grand Traverse Bay. According to this simulation, these ﬂuxes have increased
dramatically during this time period. Although a few cells have Signiﬁcant increases in
ﬂux, the majority of the ﬂux comes from a high percentage of cells with small ﬂux
increases. This simulation is based on steady state ﬂow rates and solute inputs; thus the
differences are entirely due to the solute transport times. Figure 12b shows the expected
time lag between halite application, and the corresponding asymptotic chloride ﬂux of

roughly 3640 kg/day to the bay due to an chloride input of 70 mg/l.

30

DISCUSSION

Groundwater ﬂow and transport processes need to be described to effectively
interpret the inﬂuence of land use on water quality. We have demonstrated that a
regional scale groundwater ﬂow and solute transport model can help separate the impact
of different land use related solute sources in a watershed. This work also shows that
groundwater is a key linkage between land use and surface water quality.

Groundwater models provide a critical component to evaluate the linkage between
land use/cover scenarios and water quality. The halite scenario simulates transport of
halite from roads through groundwater into streams. This method is extremely powerful
for two reasons. First it allows us to evaluate the magnitude of solute inputs by providing
simulated values of stream water concentration, that can be compared to measured
values. Second, it helps determine the relative importance of speciﬁc land use/cover
contributions to measured solute concentrations. In our example chloride from the
dissolution of halite applied to roads did not turn out to be the sole contributor to chloride
concentrations in streams. Chloride from urban sources, such as septic systems, and
subsurface brines likely contribute to the chloride concentrations in streams.

The methods developed in this paper use groundwater Simulations to predict long-
term solute loadings due to land uses/covers. We estimated groundwater ﬂuxes into
Grand Traverse Bay by summing all cells bordering this water body. This information
can be used to estimate contaminant loadings into the bay. Our calculated discharge of 8
cubic feet per second was somewhat higher than earlier estimates, which has signiﬁcant
implications about the potential solutes ﬂuxes into the Grand Traverse Bay. The halite

impact scenario estimates that signiﬁcant ﬂuxes of chloride will enter Grand Traverse

31

Bay using a reasonable value for halite input. In addition, the simulations Show the long
legacy of land use in predicted solute ﬂuxes to surface water. An ability to estimate the
impact of past and current land uses provides land use managers with a powerful tool to
examine the impact of past and future management decisions. Non-conservative solutes,
such as nitrate, can also be evaluated using the developed methods.

The model was also used to constrain watershed boundaries, which has
implications to the magnitudes of estimated solute ﬂux to the GTBW. Surface water
divides and groundwater divides often mimic each other in areas of high topography, but
in the GTBW, there are areas with relatively low topography where these divides differ.
If solute ﬂux predictions use the incorrect divides, the estimates would be biased by the
amount of solute applied to the region between the divides. Therefore proper
development of groundwater divides ultimately leads to more accurate water and solute
budgets for the watershed. This technique can also be used to decipher groundwater

source areas for to smaller regions of a watershed.

32

CHAPTER 2: COUPLING A THREE-DIMEN SION AL GROUNDWATER FLOW

MODEL TO A SPATIALLY EXPLICIT LAND TRAN SF ROMATION MODEL

33

INTRODUCTION

Groundwater ﬂow processes are often neglected in studies that examine the
inﬂuence of land use/cover changes on the environment. Likewise, regional-scale
groundwater ﬂow models rarely incorporate the effects of vegetation and land cover on
the distribution of recharge/evapotranspiration (ET) in their models, despite their
observed effects on groundwater and surface water dynamics (Thorbum et al. 1991;
Pucci and Pope 1995; Taniguchi 1997). This paper will address these issues by
simulating land use changes and their associated inﬂuences on groundwater, using
literature derived relationships between groundwater processes and land uses/covers.
The developed approaches will be able to predict the magnitude of future land use/cover
changes, as well as the resulting changes in hydraulic head distribution.

Land use/cover changes have caused various observable effects on watershed
processes. These include changes in: 1) streamﬂow (Bosch and Hewlett 1982; Hibbert
1971; Rothacher 1970); 2) chemistry of major ions (C1) and nutrients (N and P) in
streams (Herlihy et a1. 1998; Field et al. 1996; Flintrop et a1. 1996); 3) stream ecology
related to changes in stream ﬂow and temperature (Wiley et a1. 1997); and 4) the
magnitude of groundwater discharge and the distribution of recharge and discharge areas
(Pucci and Pope, 1995). Although these studies have documented the inﬂuences of land
use/covers on watersheds, few modeling studies have described the processes responsible
for these inﬂuences. Dunn and Mackay (1995) presented an example of a successful
incorporation of land use/covers into models. Their model, termed SHETRAN, evaluated

the inﬂuences of changes in evapotranspiration (ET) in selected catchments of the Tyne

34

Basin in Northeast England. They estimated ET rates using climatological data and
information on vegetation cover, which decreased ET and increased streamﬂow.

Groundwater ﬂow models are commonly developed both to examine regional
ﬂow patterns and to establish boundary conditions for smaller-scale site-speciﬁc models.
Such models commonly neglect the inﬂuences of speciﬁc land uses/covers on
evapotranspiration and recharge values. ET and recharge are commonly estimated by
examining the steady state mass balance of a watershed (Anderson and Woessner, 1992).
These methods use reasonable average inputs, but do not adequately represent spatial
variability in such ﬂuxes, which can be important for regional scale groundwater studies.

Spatial variability in ET and the resulting variations in effective groundwater
recharge have a direct inﬂuence on water table elevations. ET and effective recharge are
altered by land-use/cover changes (Taniguchi, 1997). Thus, land-use/cover changes must
inﬂuence water table elevations. Relationships between ET/recharge and groundwater
levels are well established. In order to test the hypothesis that land-use changes will
affect the distribution of water table elevations, we must infer how different land
uses/covers effect water table recharge.

Pijanowski et a1. (1995; 1996; 1997) have developed a conceptual modeling
framework called the Land Transformation Model (LTM) that describes how policy, as
well as socioeconomic and environmental factors, drive land-use change. This GIS-
based model includes a variety of land use change driving variables such as: land
ownership; local and regional population change; economics of land ownership;
transportation effects on urbanization; and local and farm-level agricultural economics.

The model calculates driving variable outputs, integrates all driving variable transition

35

probabilities, and infers future land use maps for the entire watershed. It projects land
uses in watersheds over roughly a 100-year period with annual time steps. The original
GIS-based Land Transformation Model (Pijanowski et al. 1995) was recoded in ArcView
Avenue and redesigned and integrated with a variety of spatial-temporal tools. This new
LTM includes a variety of historical land use/cover databases; a Land Use/Parcel Query
and Analysis; and a graphical user interface. The land use drivers are integrated by one
of several methods including a neural network simulation procedure and multi-criteria
evaluation. A set of spatial-temporal calibration tools is also integrated into the LTM
Toolbox, as will be discussed in a future paper.

Although other studies have modeled relationships between land use/covers and
hydrologic processes, no known studies have researched these effects using real or
predicted land-use changes. This study is unique because we have measured land use as
shown in the land use coverages as well as the ability to predict future land use changes
related in the watershed. The studies mentioned above, which have modeled these
relationships, have observed changes in the hydrology due to drastic land use changes,
such as the complete clearing of a forest to grassland. We can Show that changes in
recharge and ET, due to actual and predicted urbanization in the watershed, could have

considerable effects on water table elevations.

36

APPROACH

We have developed a straightforward way to relate land use/cover changes with
groundwater ﬂow processes. This involves assigning different percentages of the
region’s precipitation to groundwater recharge by land use/cover zones. This assumes
that the losses due to ET and overland ﬂow are associated with land use/cover, and the
remaining water becomes effective groundwater recharge. ET is likely the most
important of these factors because of its large range of variability for land uses/cover
types. ET incorporates both transpiration of moisture by plant activity and evaporation of
water from the soil and plant surfaces, which will both be affected by land use/cover
type. As a result, it is a very difﬁcult process to quantify for areas with diverse
vegetation types. For example deciduous and coniferous forests will have higher annual
ET values than pastures, because of their greater leaf cover.

Effective ET values can be estimated using mass balance techniques, remote
sensing methods (Mauser and Schadlich 1998; Sucksdorff and Ottle 1990), and modeling
studies. Our goal is to simulate the inﬂuence of spatially variable water table recharge on
groundwater levels by assigning ET percentages of the annual regional precipitation to
different land use/cover types throughout the watershed. Table 1 shows literature derived
annual ET percentages for land covers present in the GTBW. These percentages will be
applied to similar land use/covers found in the modeled region (Table 2).

Recharge is the process by which precipitation is transported to the water table.
The permeability of a soil and the magnitude of precipitation falling on the soil surface
govern how much water will be recharged to the water table. The relationship between

recharge and our land-use/land cover database depends on physical properties exhibited

37

 

 

 

r om. on carom

m - - l .223

F - ov :, 38:35 “memo“.

we on om 3:05:80 was $5289 23%“.

E. .3 om l @9333 23 5.4.0

om - cm on A390 26.. £2323 £29.89 23.39:?

m or om 53.3

3.; m9< Each 3 omcmcoom .x. .5 .macc< 569mm: 9.3
$356 on“

E 958 um: 22 some he «cocoon 2.: we :25 we moose—:86 26¢ 05 5 com: owe—£02
use Hm mo mowﬁcoeoa 3:52 .owoenooc 93 Pm mo mowsgocom 3:52. "N ~35.

 

 

 

 

00:23 .225 38: .288”. use 3.83. 8: @380 5;, $8QO 26280 V Geo”.
- 28.388222- Cam: .2938 Em elmm___mm 2... $380
cozngoem 389.... 88.85528» 38: ..m>mem:m can, zoom on: $5.0 50:8 .2355 890
cozﬁased omwgm Soc 85339 88: 5355 can zoom- owl - , 233a
EwEEoo. 858 as E .352 cozsomm>

.899

:ouﬁowg 38086 com cosmoamcmboaogo mo mos—g 33.8w 83885 .323 couﬁamcwbomgm “n 239—.

38

by the land use/cover, such as the impermeable area used by the land use/cover. An
example of a physical property that a land use/cover could exhibit is its ability to deﬂect
or prevent water from entering the soil proﬁle. For example, a parking lot will prevent
more water from entering the soil proﬁle than open grassland. Urban land-use/land
covers (which comprise 6% of the total area of the watershed) will have lower recharge
rates to the groundwater table with corresponding increases in overland ﬂow to streams
depending on the engineering applied to the area.

The recharge percentages presented in Table 2 were developed by considering
properties of the land use/cover discussed above. Recharge percentages for land
use/covers with higher impermeable areas are much lower than land use/covers with less
impermeable areas. Variable soil properties were not incorporated into the calculation of
these percentages due to incomplete soil coverages for the watershed. The recharge
percentages were determined assuming that relatively permeable soil underlies the land
use/cover. Multiplying recharge percentages by a factor that describes the soils ability to
transmit water can ultimately incorporate soil properties into these calculations. An
example of this would be a case where a very permeable sandy soil occurs where a high
recharge value from a land use/cover existed. This soil would transmit almost the
maximum amount of recharge, so 100% or 90% would be multiplied by the recharge
value for this area. On the other hand a clay loam, which has very low permeabilities,
may have a factor around 30% or 40%. Since we are concerned with just the effects of
varying land use/cover, the spatial variability of soils in the GTBW averaged over a 100-

meter cell were assumed to not have a great inﬂuence on the recharge percentages.

39

The percentage of the water that becomes streamﬂow, via interﬂow and overland
ﬂow, is determined by the amount of water left over from the annual ET and recharge
percentages. The ability of a land use/cover to direct water available to streamﬂow is
embedded in the recharge percentage value. Of the total precipitation falling on the
watershed the least amount annually becomes overland ﬂow, although these values can
change depending on how land use management practices. Our groundwater ﬂow model
does not simulate stream ﬂow but maps showing the potential source areas of streamﬂow,

such as Figure 13, may be developed.

40

 

 

 

 

Overland Flow High

Low

Streamflow distribution. Cells with a lot of impervious, urban land

use/cover, surfaces will direct more water via overland flow to streams possibly

causing higher erosion rates.

13

Figure

41

STUDY AREA

Michigan's Grand Traverse Bay Watershed (GTBW) is a good site to simulate the
inﬂuence of land use change on water table elevations by coupling the LTM with a
groundwater model. The GTBW, located in the northwestern portion of Michigan’s
Lower Peninsula (Figure 1), is currently undergoing some of the most rapid population
growth and land use change of any region in the United States (Vesterby and Heimlich,
1991). During the last 10 years, the watershed increased its resident population by over
10%, making it one of the fastest growing areas in the Midwest. Traverse City is the
largest city in the watershed, with a resident population of approximately 100,000 and a
summer recreational population often exceeding $00,000. We are also examining the
effects of these land use changes on water quality, which will be presented in a later
paper.

The watershed receives 105 cm of annual precipitation, of which approximately
40 cm is recharged to the water table, 50 cm is evapotranspired, and 15 cm becomes
streamﬂow (Cummings et al., 1990). The Boardman River, which drains about 60% of
the entire watershed, empties into Grand Traverse Bay near Traverse City. Two large
lakes, Elk Lake and Torch Lake, drain the rest of the watershed into Grand Traverse Bay.
Grand Traverse Bay is attached to Lake Michigan and receives most of its water from the
watershed’s streams. Most of the annual discharges to the watershed’s streams (95%)
come from groundwater sources. Groundwater ﬂow in the study area is from areas of

high topography towards the Grand Traverse Bay. In the southern portion of the

42

watershed, groundwater velocities tend to be high because of high conductivity sediments
and high gradients.

Land use in the GTBW is predominantly forest (49%) and agriculture (20%).
Much of the forested portions of this watershed are managed within the Pere Marquette
State Forest, which encompasses most of the upland reaches of the Boardman River
tributaries as shown in Figure 1. Urban land use comprises 6% of the total area of the
watershed, with the other main land covers are open herbaceous/shrub/grasslands (15%),
water (9%), and wetlands (1%).

A 3-dimensional groundwater ﬂow and transport model (See Chapter 1 for
complete description) has been developed and calibrated for the GTBW that incorporates
recharge and ET as constant values. The model simulates groundwater ﬂow through the
glacial sediments that lie on top of older Paleozoic bedrock as discussed in Chapter 1.
The Paleozoic rocks, which are mostly shale, act as good barriers to vertical ﬂow in the
region. The modeled region encompasses an area of 3300 km2 with 500,000 100-meter
horizontal cells divided into 2 layers. The computational resources for data processing a
model with more than two layers were prohibitive. The 2-layer model provides a
reasonable representation of measured head data (Figure 6). Important features to note
on this map include the large discharge region of the Boardman River and the high
gradients in the southern portion of the watershed. Overall the conductivity values in the
model are higher in the southern portion of the watershed and become lower northward.
The lower conductivity values in the northern portion of the watershed suggest slower

travel times.

43

The Land Transformation Model has been developed to forecast land use patterns
for the GTBW to the year 2030 (Figure 14). Our GIS-based Land Transformation Model
(LTM) uses neural network simulation algorithms that integrate our spatial-temporal
drivers with historical land use databases. This novel technology allows the computer to
“learn” how socioeconomic drivers affect documented land use changes. These patterns
are then used to forecast land use assuming a “business as usual” scenario. We have
found that the LTM predicts urban change with an 89% accuracy using 100 m cells. This
cell size was chosen to allow easy interaction with the groundwater ﬂow models. Areas
where land use change are not well predicted include new recreational areas (mostly golf
courses) and large residential subdivisions, which will be addressed by incorporating new

drivers that characterize these changes.

44

.meouﬂsemm 2: 5 68: 82> moms—38 82: mo 8:: =< .6on

gang—828% v53 2: maﬁa @8828» mowﬁo>ou 8: US: 23% 380.8%” 8m Saw—m 2: mo .228985
SPA on 98 So.» 2 .266 28: Emu 2: so muwﬁo>oo 8: v5: 93 2F .mouona Eton 88m Bigot om:

2:: 89 2: mm a2 2: no 03550 8: 9:2 BE. .BmHO on“ 5 33-32 38:65 can E056 u: 9.53%

 

 

z a 850 02 D 8.3.6”. I
<6! on or 0 or Show I ucﬂ :80 D

mucm=o>> voumocom D 233:? D
£28980: a .225 I ascogmom D

 
  
 

 

 

 

 

  

45

 

83-23 ovom , ,. . _ 33.23 88

 

METHODS

Studying the inﬂuence of land use on watershed hydrology requires the use of
large datasets that often need signiﬁcant preprocessing. In addition, managing and
analyzing these large spatial databases is difﬁcult with traditional modeling tools. As a
result, we used Geographic Information Systems (GIS), speciﬁcally ARC/INFO (ESRI
1998a) and ArcView (ESRI 1998b), to manage, manipulate, and analyze the spatial data
for our groundwater ﬂow model. F ortran 90 codes were also written to perform
operations on gridded datasets for model input.

Our land use database for the GTBW was obtained from the Michigan
Department of Natural Resources MiRIS Ofﬁce. MiRIS is a statewide land use database
developed in the late 19703 from l:24,000 aerial photography. Land uses were classiﬁed
at Anderson Level I, which include: urban, agriculture, nonforested vegetation
(grasslands, shrubs, herbaceous), forests, wetland, open water, and barren. Land use data
for the GTBW region were acquired by township and compiled using ARC/INFO.

The three land use/cover datasets, shown in Figure 14, used in the study were
gridded at 100 m horizontal cellsize. Each distinct land use/cover was given a unique
numerical identiﬁer which stayed was the same for the current, 10 year projection, and 50
year projected land use/cover datasets. A Fortran 90 routine was written to convert the
gridded land use/cover datasets into recharge grids. Each land use/cover dataset was
converted into an ASCII ﬁle format, where each land use/cover dataset was associated
with an ET / recharge percentage. A single ET percentage was assigned to each 100-
meter cell, which was considered to have a homogenous land use/cover. The ET

percentage was then multiplied by the total annual precipitation to yield an average

46

annual ET for the watershed (Equation 2.1). The ET value calculated above was then
subtracted by the average annual precipitation and multiplied by the recharge percentage
to yield a recharge value. The recharge supplied to the model is a percentage of water
leﬁ over aﬁer the ET processes occur. Lastly the streamﬂow component is the sum of

water left over aﬂer ET and recharge processes have occurred for that land use/cover.

 

ET = Pr ecipa*(E72/o)
Rec .—. [Pr ecipa — ET]* Re c% 2.1
Oﬂw = Pr ecipa — [ET + Recharge]

 

 

 

where:

ET — average annual Evapotranspiration

Precip — average annual Precipitation

E T a, - E T percentage factor for land use/covers

Rec — average annual Recharge

Rec% - recharge percentage factor for land use/covers

Oﬂw - average annual overland ﬂow to streams

The recharge values calculated were then input into the groundwater ﬂow model

and ran until a steady state solution was achieved for that time step. Since MODF LOW
does not simulate ET that takes place above the saturated zone, we were limited to
remove the calculated ET rate and only provide the model with the calculated recharge
value. This process was repeated for each land use/cover time step. Head maps for each
land use/cover time step were saved and differential head, or drawdown, maps were
calculated for current time against the preceding time. These maps were created by

subtracting the head values in a cell of one simulation by the preceding simulation. In the

case of the current land use/cover dataset, it was subtracted against the constant recharge

47

case to provide a map showing how the incorporation of variable recharge as the result

land use/change will change the groundwater heads in the watershed.

48

RESULTS

Positive numbers in the head difference maps correspond to increases in head
over the simulation period and negative numbers correspond to decreases in head. The
maps were calculated in the following order: the 1990 land use/cover simulation was
subtracted by the constant recharge simulation (Figure 16), the 10 year projected land
use/cover ﬂow simulation was subtracted by the 1990 land use/cover simulation (Figure
17), and the 50 year projected land use/cover ﬂow simulation was subtracted by the 10
year projected land use/cover simulation (Figure 18). Two maps were also developed to
show the sensitivity of the model to positive and negative changes in recharge. This map,
Figure 15, shows areas that will change in head due to overall changes in recharge in the
watershed. Areas on these maps that are not green show a change in head relative to the
constant recharge simulation. These areas are deciphered as regions that are extremely
sensitive to changes in recharge and should not be used in the interpretation of the land
use/cover simulations.

Incorporating land uses/covers into groundwater ﬂow models provides a powerful
approach for determining their inﬂuences on watershed hydrology. Recharge values are
often difﬁcult to measure and predict thus, constant values are oﬁen used in groundwater
ﬂow simulations. The incorporation of land uses/covers into a groundwater ﬂow model
allows users to specify a spatially variable recharge, based on real data. Figure 16 shows
the changes in head between the integrated 1990 land use/cover simulation and the
constant recharge simulation. This simulation shows the head differences expected due

to variable recharge values associated with different land use/covers (Figure 16). This

49

66:808.,» 2055 2: mo gouge—o 638
Houaacasohw 2: Soto EB owbﬁoE E mowamno =25 Bo: 952i 2st “EH .969“ ES: 2: :o 5on 8: 8m
:2: 3on E omega“: 5 mowcmno 9 03:38 E 38E 3o: 53323on 22. 523288 owcmaoom "mu 95w?—

3.3.5:—
ems—Evy—

c\oc —

 

50

[OW

 

-4 —2 O 2 4

Figure 16: Current land use - constant recharge simulation. The head difference map
for the variable vs. constant recharge simulations show that an overall increase in head
occurs but localized decreases in head are found near the Traverse City region.

[OW

 

 

 

Figure 17: 10 year land use projection - current land use simulation. The
change in head between the 10 year projected land use and the current land use
simulation show changes in the gradient developing as well as lower in the
quickly urbanizing areas.

52

 

-4 -3 -2 -1 0 l 2 3
Figure 18: 50 year land use projection - 10 year land use projection simulation.
The difference in heads developed in the 40 years between the two simulations
still shows the inﬂuence that the Traverse City area will have on groundwater
elevations.

53

simulation shows that heads are higher in most areas of the watershed. This result is
related to the large agricultural and non-forested areas that have higher recharge values
than the average simulated constant recharge value. In contrast, the variable recharge
heads are lower in areas where the recharge value is smaller than in the constant recharge
simulation, such as urban areas. Notice that heads around Traverse City (see Figure l
for location), and neighboring areas are lower in the variable recharge simulation due to
the high percentages of impervious surfaces in these urban areas (Figure 14).

The head difference map for the lO-year projection-current land use/cover (Figure
17) examines the potential effects of predicted future land use change on water table
elevations. The current version of the LTM only predicts changes from any land
use/cover to urban, so the changes in groundwater heads will show only the effects of
urbanization. These effects can be noticed very clearly in areas surrounding Traverse
City, which the LTM predicts increases in urban land uses in the next ten years. The rest
of watershed shows little or no decrease in head values, besides areas that are highly
sensitive to changes in recharge, over the same time period.

The changes in head between the 10-year projection and the 50-year projection
are slightly larger in their inﬂuence (Figure 18). Changes are once again found to occur
in the Traverse City area, but a new prominent decrease in head occurs directly east of
Traverse City. The heads here are much lower than directly beneath Traverse City. The
large changes east of Traverse City can be attributed to a large urban sprawl predicted by
the LTM (Figure 14). This reduction in recharge is directly within an area that is
somewhat conﬁned from the regional ﬂow ﬁeld. The Boardman River runs south of this

region and is of considerably higher head than this area, whose ground elevation is

54

considerably lower. This in essence isolates this region from the other areas and limits
the ability of this area to receive groundwater. Another reason that this area has lower
heads than in the previous simulation is the fact that this region of urbanization was not
as large in the 10-year prediction. This results in a large decrease in head due to less
recharge.

Changes in groundwater head between these two times have an effect on the
regional gradients to the west of Traverse City and the eastern portion of the watershed.
These changes in head are a direct result of the model’s sensitivity to recharge in these

areas.

55

DISCUSION

Simulating land use changes and their associated inﬂuences on groundwater
processes is an important step in determining the impact of these changes on watershed
hydrology. We have already shown (Part 1) that groundwater processes play an
important role in interpreting the inﬂuences of these changes on water quality. Therefore
in order to accurately describe groundwater processes, for the correct interpretation of
water quality data, land use/cover change must be integrated into the groundwater
simulation. In this last section we have shown that the integration of these changes is
possible and they do have large effects on the distribution of simulated groundwater
heads and gradients across the GTBW.

The variable versus constant recharge simulation demonstrated the impact of
variable recharge on simulated heads. This simulation also depicts changes in head as the
result of the redistribution of recharge areas due to the incorporation of high impervious
surfaces in urban areas. Changes in groundwater heads as the result of 10 years of ﬁiture
land use change show the effects of urbanization near the Traverse City region. The
lowered groundwater heads reﬂect that more water in these areas is directly ﬂowing into
Grand Traverse Bay through streams or sewer systems rather than groundwater.
Signiﬁcantly lower heads just east of Traverse City, in the 50 year versus 10 year
projected land use/cover simulation, are indicative of the potential consequences of urban
sprawl on wateshed hydrology. The reduction in recharge to these areas plus the
geographic location of this region causes groundwater heads to drop considerably

between the two time periods. This simulation also depicts the changes in regional

56

groundwater gradients as a result of land use/covers and the model’s sensitivity to small
recharge changes. Changes in regional gradients can disturb current ﬂow paths and
increase ﬂow velocities.

The land use/cover datasets that are being used in this simulation are observed
and predicted. The observed land use/cover dataset allows us to calibrate our model
using real data to obtain more accurate estimates of aquifer parameters, such as hydraulic
conductivity. The forecasted land use/cover datasets have allowed us to infer the future
changes on groundwater heads as a result of land use changes. This is a powerful tool
that can be used in unison with other approaches, such as solute transport models, that
attempt to delineate areas in watersheds that are sensitive to small or large land use/cover
changes.

Together these simulations underscore the need to better approximate and
understand anthropogenic inﬂuences on watersheds. We have shown that considerable
decrease in groundwater heads will occur if current land use practices continue. These
changes were predicted only from the consideration of varying land use/covers. The
simulations do not integrate the inﬂuences of groundwater pumping, which in its breadth
and scale may have a larger effect on heads. Overall the changes and redistribution of
recharge areas as a result of urbanization may have a considerable effect on the simulated

water table elevations and gradients in the Grand Traverse Bay Watershed.

57

APPENDICES

58

APPENDIX A: GEOLOGY OF THE GRAND TRAVERSE BAY WATERSHED

INTRODUCTION

The sole reason for writing this section of the thesis lies in the amount of time that
was spent researching and thinking about the geology of the watershed. The objective for
researching the geology of the watershed was to develop a hydrostratigraphy of the water
table aquifer for the input into a 3-dimensional ﬂow and transport model. This section
could serve as a guide for someone starting a project similar in scale to the GTBW, in
which they might learn something from mistakes and assumptions that were made
throughout the 2 years of research. Unfortunately there was little information on the
glacial deposits of the region with the exception of the drumlinized till plains in the
Charlevoix region. A lot of information was found on the bedrock of the region, because
of the exploration interests of oil companies, but this information was not relevant at this
point. Professor Grahame Larson proved to be an excellent resource on the glacial
geology of the region, showing just where to search. This section will begin with a
description of the oldest rock unit subcroping in the GTBW and progress forward into the
present depositional systems noting mistakes made along the way and what was done to

correct them.

59

BEDROCK GEOLOGY

The bedrock of the GTBW was not researched at great depth in this study because
we are only concerned with the ﬂow and transport through the water table aquifer. The
bedrock serves as a barrier unit to the ﬂow and transport of groundwater in the glacial
materials. Middle Devonian rocks, deposited on a carbonate platform, are the oldest
rocks found subcroping beneath the watershed, although older rocks are found beneath
these rocks. The Devonian Period in the GTBW was a time of widespread but ﬂuctuating
seas, extensive limestone and dolomite, and shale accumulation, and some evaporite
deposition (Dorr and Eschman, 1970). Near the end of the period black muds of the
Antrim Shale were carried into the Michigan Basin from uplift in the Appalachian region.

The ﬁrst unit subcroping in the GTBW is the Middle Devonian age Traverse
formation. The Middle Devonian Traverse group is mostly fossiliferous limestone and
dolomite that is commonly light olive-gray to brownish gray, ﬁne-grained becoming silty
upsection. The limestone is occasionally fractured and some slight deformation
structures are observed in the uppermost units, the units above mainly the Jordan River
formation (Kesling et al., 1974). Correlated limestones found in the Alpena area
commonly have karstic characteristics, including sinkholes and conduit features. The
Traverse Group is found outcropping and subcropping beneath both the Leelanau
peninsula and in the Charlevoix region (Figure 19). Along the shore of northern East
Grand Traverse Bay north and south of Norwood, Michigan (shown as the blue star on
Figure 19) the Middle Devonian Traverse Group and overlying Upper Devonian Jordan
River Formation (informally referred to as the “Norwood Shale”, by Ehlers, 1938) and

the type section of the Antrim Shale outcrop. Other exposures of the Middle Devonian

60

.mEBva Emankoc PE: 8 8a magma .3 @8308.“ bﬁmzm m_
:ocancom 3822p. own :m_:o>oQ BE. .aBmE 05 E 30¢ 32:? 9 SEE boom 5 825i
.29: mo 352: mEmaoo 533 5:50 05 5 £863 BE. .3280 xooﬁom "a 9...»:—

 

 

 

hrdr
1

9390 3.36;
5.6960 23.2

wlcw ESE
cscclo .233

2.5 55me
863902559353}

 

323:3 Ezra} .....
coaemmﬁmi Logo:

.3022} 9.3

 

 

 

 

61

Traverse Group in the watershed are found just south of the town of Charlevoix and
along the shore of the Leelanau Peninsula. The Whiskey Creek Formation, which lies
beneath one of our stream sampling sites, was examined at outcrop and found to be gray
to brown limestone, contained some dolostone and chert and some unknown shelly
fossils. The Jordan River, which lies between the Traverse Group and the Antrim Shale,
is a formation composed of very thin limestones and nonresistant shales, about 1 to 2
meters thick (Gutschick and Sandberg, 1991a).

The transition from carbonate platform to deep basin was a time of dynamic
change produced by the Acadian orogeny. The deep subsidence toward the end of the
transition period had completely changed the seawater and sediment regimes in the newly
formed Late Devonian Michigan Basin and produced a stagnant, black-mud sea ﬂoor
devoid of oxygen (Gutschick and Sandberg, 1991b). The type locality of the Antrim
Shale, deposited during this time, is 1.6 kilometers south of Norwood. The Antrim Shale
is mostly a black ﬁssile pyritic shale with thin fossiliferous limestone beds in the lower
Norwood Member, and transitions into a “light-colored calcareous unit” in the Paxton
and Lachine Members (Gutschick and Sandberg, 1991b). The Antrim Shale in the
GTBW is overlain by the Ellsworth shale whose type section is in the Petosky Portland
Cement Co. quarry pit , 2.4 km south of the town of Ellsworth, Antrim County, Michigan
(yellow star Figure 19). The Ellsworth is greenish gray (Figure 20), bluish gray, light
gray, and dark gray shale with some sandy and oolitic limestone interbeds in the
subsurface (Gutschick and Sandberg, 1991b). The base of the Ellsworth is deﬁned as the
change from brown or black Antrim Shale to greenish gray shale; and the top is an

unconformity overlain by Mississippian rocks, a red limestone at the base of the

62

cam—o 35302 a on 9 Ewsoﬁ ma? amount £5 bEEwtO .834 35:00 go ﬁne
53 Easy 3 388:0 35 29% 55325 cﬁqﬁmmumaz 05 mo mm 20: 230% 2t. .3830 oﬁsm 5N charm

 

63

Coldwater shale. The facies relationship between the Antrim and Ellsworth is an
intertounging or alternation of greenish gray deltaic Ellsworth with the basinal, brown or
black Antrim. The location of the Antrim and Ellsworth Shales are found in Figure 19.
These Shales underlie a signiﬁcant portion of the GTBW.

As was stated above, the top of the Ellsworth shale is marked by an unconformity
with the lowest member of the Coldwater Shale Group. The Mississippian Coldwater
Shale Group is only known through subsurface drilling in the GTBW (Figure 19). The
Coldwater is a gray shale with maximum thickness of up to 500 meters in the Michigan
Basin with local thickness in the GTBW up to 200 meters. The Marshall Sandstone,
which lies above the Coldwater, subcrops only slightly into the groundwater modeled
domain. The contact between the Coldwater and the Marshall is gradational, as the
Coldwater becomes coarser near its top. The Coldwater provides a nearly impermeable
barrier to vertical ﬂow throughout the whole state where the Marshall overlies it. It also
appears to be a good barrier to vertical ﬂow in the case where glacial deposits lie on top
of the Coldwater, which occurs in the GTBW. The Coldwater in the subsurface is often
confused with glacial lacustrine clay deposits. It some cases the shale and lacustrine clay
(that has the Coldwater Shale as it’s source) are indistinguishable. This provides a
problem when trying to deﬁne the elevation of the bedrock in the GTBW. Many times it
is just a blind guess as to where the bedrock starts. Occasionally the driller’s would note
a weathered zone that could be used to identify the top of the bedrock. This did not prove
to be too much of a problem since the hydrostratigraphy in these areas was deﬁned to

have thick low conductivity lacustrine clay above the Coldwater shale. The only problem

64

 

this raised is when trying to interpolate the bedrock elevation grid for the watershed this
caused small mounds, or an irregular surface, to develop in the output.

Overall four distinct bedrock units directly underlie the glacial sediments of the
GTBW. From bottom to top, or oldest to youngest, they are the Middle Devonian
Traverse Group, the Upper Devonian Antrim Shale, the Upper Devonian Ellsworth Shale,
and the Mississippian Coldwater Shale. A typical North-South cross-section of these
rocks in the GTBW is shown in Figure 21. All of these units provide a very good barrier
to vertical ﬂow, with exception where the Traverse Group is possibly fractured and this
area was quite small.

The bedrock elevations for these units were calculated from the 7,000 oil and gas
wells throughout the watershed and the few water well logs that reached the shallow
Traverse Group in the northern parts of the watershed as well as the few outcrops along
the shore of the bay and in Ellsworth, Michigan. All of the data for Oil and Gas wells
were obtained from Steven E. Wilson at the DEQ — Geological Survey Division.
Bedrock elevations were also obtained from hand drawn bedrock elevation maps
acquired from the Michigan DEQ — Geological Survey Division (Millstein, 1983; Reska,
1983; Elowski and Bricker, 1982; Bricker and Lilienthal, 1976). Water well logs and all
other related information were obtained from David F orstadt at the Michigan Department
of Natural Resources -Real Estate Division. The data for the elevation of the bedrock
(Figure 2) in the GTBW was interpolated using 2-dimensional ordinary kriging (Figure
22). The areas where the map has mounding represents kriging calculating a change
from the mean of the dataset to the value located at this spot. One reason for this odd

feature are not enough data points surrounding a local high in the bedrock causing a

65

.vonmcouaB on: E 0868?. 83380985
EthE HEB ﬂuoﬁmuom Bo: 885883 Paw—m we .33“ 550% mm 2350 05
me 8:08 $0.5 ﬁscmﬁtoz 3293 < .2»th 06 we 8:08-895 8:8 8 8.52 "mu 9...»:—

 

 

Sﬁsau. mania. m 5 in 5350 2285: a. oo

0.25 E35. nun—.2 05922 :95... an. .05.. , ., EC

296 seesaw m. o: 05222 9520 I 6

22m 53.8 I o: 23.0 258%.. D no

283:8 32.: 22.5 ea £28 05883 I and
Bio: .55: .E .35 .3“? I

 

 

 

 

66

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

   

 

l
'i
l
Bedrock Elevation
I 7412-32627 I 133483-158 703
I 32 627-57 843 I 158 703- 183918
I 51 843~83 058 E] 183 918-209 133
[:l 83 058408 273 |:] 209133-234 349
U 108 273-133488 ,
7 0 7 14 Kilometers A i
N

 

nu. C'-

 

Figure 22: Bedrock elevation. Bedrock elevations were interpolated using a

kriging technique from Oil and Gas wells, outcrops, and the few wells that

penetrated the bedrock.

67

 

mounding effect or a bad elevation value for this well causing a deviation from the mean
and other surrounding points. This ﬁrst explanation has some validity because quality
checked data on hand drawn maps show the same features as well as Melhorn (1954) and
Leverett and Taylor (1917) describe bedrock knobs common throughout the region. The
second explanation has little weight because most of the data from oil and gas wells are
quality checked though misinterpretation of a lacustrine clay for a shale could happen.
There are areas where the bedrock will have elevation changes of 100 meters in the
vertical and laterally in distances of less than a kilometer. In the southern portion of the
watershed, where more data exists, these features appear as bedrock valleys, which could
be either preglacial drainage patterns or syn-glacial erosional features. Some of these
valleys show the same trend as the large linear lakes in the watershed, notably Elk Lake
and Torch Lake (See Figure 22). Overall the bedrock elevations are much higher in the
northern portion of watershed and slope along the same trend as dip of the Michigan

Basin, which is toward the southeast.

68

GLACIAL AND RECENT DEPOSITS

The stratigraphic descriptions above indicate that within the region there are
formations of greatly different hardness and lithology. These initial variations probably
produced a preglacial topography in which a northward facing cuesta had been developed
on the Marshall Sandstone (Melhom, 1954). This cuesta faced a lowland belt to the north
developed on the softer Antrim, Ellsworth, and Coldwater shales. The Traverse, as well
as younger Devonian limestones, formed another upland district toward the north. The
presence of sinkholes in the Alpena area indicates chemical weathering had been initiated
on this upland. Glacial erosion and deposition modiﬁed this pattern but did not
completely obscure it, so that the topographic relationships are still shown on the present
surface. For example, the Port Huron moraine throughout the GTBW follows the
approximate position of the scarp fronting the Marshall Sandstone cuesta (Figure 21).
This cuesta may been have been a controlling factor in morainal deposition in earlier
times, acting as a focus of building for the Port Huron moraine during the Late
Wisconsinan. As was stated before, the Coldwater and even the Antrim and Ellsworth
were important contributors because many of the tills and clays of the region are known
to be derived from these shales.

Figure 23 shows this high accumulation of sediment in the southern portion of the
watershed whereas in the northern portion the drift is very thin. This ﬁgure can be
visualized as a wedge of sediment, where the sediment is thickest in the south and thin to
almost bare in the north (Figure 21). Sediment thickness in the watershed range from
over 250 meters in the bedrock valleys to no sediment in some spots of the northern part

of the watershed. It should be noted here that within the large lakes present in the region,

69

l'OW

 

col

0 50 100 150 200 250

Figure 23: Drift thickness. The drift thickness map for watershed was
calculated by subtracting the elevation of the surface for each point by the
bedrock elevation for the same point. Note that some prominent valleys exist in
the watershed.

70

speciﬁcally Elk Lake and Torch Lake, the drift thickness is not accurate and may even
overestimate the amount of sediment. This is because data for the depth of lakes was not
available to us and was not used in the calculation of sediment thickness. It is my guess
that these lakes are seated in the soft Ellsworth and Antrim Shales, and may contain some
lacustrine clays that were deposited in the Late Wisconsinan. If this statement is correct,
then the depression where these lakes now exist must have been present during the last
retreat of the glaciers and probably contains a large amount of sediment from the
deposition during the glacial lake stages of the Great Lakes Basin. With reported lake
depths in excess of 100 meters, this means that these depressions may have been well
below sea level initially, possibly remnants of the preglacial drainage network.
Throughout the Pleistocene there were close to 11 glacial advances in the Great
Lakes Basin. Unfortunately our sediment deposits in the GTBW only reﬂects the last of
these great glaciations. Figure 24 is a digitized 1 km scale map of the surﬁcial deposits in
the GTBW (Farrand and Bell, 1982). The deposits of the GTBW are mostly Late
Wisconsinan in age but the sand dunes found along the shores of Lake Michigan and
Grand Traverse Bay and within the Mancelona Outwash Plain are Holocene (recent) in
age. The retreat of the Wisconsinan ice sheet began about 20,000 years ago near the
northern Ohio River valley. 13,000 years BP the ice sheet stopped and started advancing
for a short period. As this ice sheet stood stagnating, this pause built the Outer Port
Huron Moraine that now forms the southern watershed boundary of the GTBW (Figure
25). The ice sheet again retreated but stopped again and formed the inner Port Huron
Moraine and the Mancelona Outwash Plain where the Boardman River has carved a huge

valley (Figure 26). In most of the watershed the Outer Port Huron is at an average 40

71

 

  
  
  
  
 
 
   
  

m.

Su'ﬁdal Geology (Bel and Fmand. 1982)
Other Areas

Peel and Muck

Poslglacml AlluVIum +

Dune Sand

Lacuslnne Clay and Sill

Lacustnno Sand and Gravel

Glacwl Outwash Said and Gravel
Ice Contact Outwash

Glacml Tlll- Fine

End Moralnos- Fine

Glamal Tlll- Medium

End Mommas - Medium

Glaclal Tlll - Coarse + %
End Mommas - Coarse
Thin to Discontlnous

Mlllcml Flll

I
E]
I
El
El
I
El
I
I
I
I
I
I
I
C]
I

Exposed Boamck

I No Date
a? Drumlins A

 

 

, .__.,I,__

30000 600-02 9%
Figure 24: Surﬁcial geology of the GTBW. The surﬁcial geology of the Grand
Traverse Bay Watershed is dominated by the high topography moraines
separated by the broad Mancelona Outwash Plain. The crossection line A-A’ is
found in Figure 19.

 

 

 

72

 

Formation of Outer Port Huron Moraine

 

(1'

 

 
 

W \ gimmw

% Ice Position
Glenwood M phase of Lake Michigan
- Saginaw phase of Lake Huron Lobe

 

 

 

 

Q

V

 

 

Q
\

, A 1.
W:

 

 

 

 

 

 

  
  

 

 

 

 

 

 

Stage 30: 13,000 BP
J l l

 

 

Figure 25: Outer Port Huron Moraine. The Outer Port Huron Moraine, which
formed between 13,500 and 13,000 years BP, rises to almost 200 meters over
Grand Traverse Bay.

73

 

Formation of Inner Port Huron Moraine

 

 

 

MW \. C W‘Y “3”"
W Ice Position
E Glenwood M phase of Lake Michigang

- Saginaw phase of Lake Huron Lobe

(1’

 

 

 

 

 

0Q i k

 

 

 

 

 

; 31.3%,.
If. '

 

 

 

l l l
Stage 32: 12,875 BP
. ,.

   

 

 

 

 

 

 

 

 

 

 

Figure 26: Inner Port Huron Moraine. The Inner Port Huron Moraine formed
the Mancelona Outwash Plain, which consists of over 150 meters of sand and
gravel that is home to the Boardman River valley.

74

meters above the Mancelona Plain and the Inner Port Huron is about 8 meters above the
Mancelona Plain. The Mancelona Outwash Plain once served as a major dicharge valley
for an enormous glacier that occupied the region just to the north and east of the plain.
Now the Mancelona Plain is home to the Boardman River, which has cut a roughly 50-
meter valley into the plain before cutting through the area where the Inner Pt. Huron
merges with the Manistee Moraine.

The crests of the Inner and Outer Port Huron Moraine trend northeast-southwest
in parallel fashion for about 40 km. In cross section (normal to trend) each moraine has a
pronounced asymmetrical proﬁle with steep northwest facing proximal slopes (Figure
21). The surﬁcial sediment in the moraines and outwash plain are predominantly
glacioﬂuvial sand and gravel, with lacustrine clays and tills limited to areas up—ice from
the moraines (Blewett and Winters, 1995). Subsurface data suggest that tills may be
more common, though thin and separated by stratiﬁed sand and silt. Blewett and Winters
(1995) suggest that the Inner and Outer Port Huron may not be actual moraines but more
likely associated with an outwash setting. Outwash plains tend to have more stratiﬁed
deposits and would more likely have higher conductivities than a morainal complex.
Leverett and Taylor (1915) have classiﬁed this complex as a moraine, even though the
calibrated hydraulic conductivities of the groundwater model suggest a higher conductive
material such as stratiﬁed sands and gravels. I have discovered thick clay and gravel
layers in some deep well logs that may suggest a morainal base to the. whole complex.

The Manistee moraine is found in the western portion, mostly southeastern
Leelanau County, of the GTBW where it appears to be closely blended with the main

Port Huron system (Figure 24). The Manistee ﬁrst appears as a massive ridge at

75

 

Manistee on Lake Michigan, and trends northward as a series of overlapping segments
parallel to the lake in the vicinity of Glen Lake in Benzie County (Melhom, 1954). Here
it turns at a right angle to reach a point near the head of Grand Traverse Bay. When it
crosses the Boardman River south of Traverse City, it again turns northeast so as to
almost intersect the Inner Port Huron ridge near Rapid City. The Manistee loses its
identity here, but it can be found extending northward in Antrim and Charlevoix Counties

marked by small isolated morainal segments separated by parallel southeastward trending

 

valleys. Melhom (1954) believes that the ice front held steady at the head of Grand
Traverse Bay while the Port Huron ice north of Rapid City was retreating westward to the r
position of the Manistee Moraine (Figure 27). The minor moraines of Leelanau County
and the Old Mission Peninsula appear to be directly related to the Manistee morainic
system. They could be interlobate ridges formed between ice bodies which occupied
Grand Traverse Bay on the East and Lake Michigan on the west. Highest points on the
Manistee moraine are just west of Traverse City between Grand Traverse and Leelanau
Counties where the moraine achieves a total width of four miles. The Manistee Moraine
shares many of the same characteristics as the adjacent Port Huron Moraines.

The best drumlins in the GTBW lie on ground moraine in northern Antrim and all
of Charlevoix Counties, which are more or less parallel to the axis of Lake Charlevoix (S
30° E) (Figure 28). A few drumlins are found east of the chain of lakes near Bellaire and
Ellsworth. In this area the forms are crowded into a narrow belt of ground moraine
adjacent to segments of the Manistee moraine. Where the ground moraine is very thin,
many of the ridges have been sculptured out of the Ellsworth Shale by ice. Between the

north end of Torch Lake and Traverse City, drumlins trend only a few degrees south of

76

Formation of the Manistee Morainic System

-n (\ n n Y” §
M Ice Position WVQA
- Lake Chicago
:1 Early Algonquin
s

V -.

 

 

(

 

 

 

{J1

 

 

 

 

 

4- u

 

 

 

1
Stage 46: 11,800 BP

 

 

 

 

 

 

 

Figure 27: Manistee Moraine. The Manistee Moraine is responsible for
forming the southwestern edge of the Grand Traverse Bay Watershed, where
it meets the Inner Port Huron Moraine.

-W

 

Figure 28: Drumlins. Drumlins, like the ones seen here, dominate groundwater
flow and stream drainage patterns in the Northeast portion of the GTBW.

east. Since this is subparallel to the depressions occupied by Torch and Elk Lake, it is
doubtful the same ice activity was the cause of the two features. Drumlins of Charlevoix
and Antrim counties are composed of very stony, sandy, buff or light brown till
(Melhom, 1954). The drumlins of the Grand Traverse region are partly composed of red
till. It has been hypothesized (Melhom, 1954; Leverret and Taylor, 1915) that drumlins
of the Grand Traverse region were formed during a minor readvance of ice onto ground
moraine laid down subsequent to building the Manistee moraine. At a later time, ice
moving out of the Lake Michigan basin advanced inland and plastered red till on older
drumlins and in some cases secondary drumlins were built on older forms.

In addition to morainal ridges of the GTBW, a large percentage of the surﬁcial
geology is made up of ground moraine (Figure 24). Ground moraine usually forms by a
quickly retreating glacier that is losing sediment by dislodgment or melting, which is
often referred to as an ablation till. Ground moraine is unstratiﬁed and unsorted and may
contain sand stringers and abandoned meltwater channels. The ground moraine in the
GTBW is very loosely compacted in some places and very heavily compacted in others.
This is due to the sand-clay percentage in the sediment composing the ground moraine.
In general the ground moraine closer to the major morainic ridges will have higher sand
content and thus have a looser texture. This could also play a large factor in the transport
of groundwater through these sediments. Old Mission Peninsula is a dissected high level
ground moraine, which reaches altitudes between 750 and 850 feet. The Old Mission
Peninsula contains a few drumlins with the same composition as those on the west or the

east.

79

 

The moraines found in the GTBW often exhibit high topographies and when
water levels were higher they were probably not submerged. This leads to the
development of highly different sediment characteristics on the moraine and on the ﬂanks
of the morainal complex. The lower areas are important to us because they can contain
different (ﬁner grained) sediments that can impede groundwater ﬂow. In examining
Figure 24 and Figure 4a areas can be chosen that have a high chance of accumulating
ﬁner grained sediments when lake levels in the Great Lakes basin were higher. These
areas are depicted on the surﬁcial geology map (Figure 24) as lacustrine sand and
gravels. But we must remember that this is only a surﬁcial geology map and does not
reﬂect the entire sediment proﬁle beneath the area depicted on the map. It is the case in
these areas that they have a large accumulation of lacustrine clays some up 100 meters in
the Traverse City region (Cummings et al., 1990). The normal shallowing upward proﬁle
in a basin ﬁlled with water has ﬁne-grained sediments at the bottom and coarsening
upwards. A good example of this is in the Traverse City area where thick lacustrine clays
grade into sands and gravels related to the progradation of a delta into the west arm of the
Grand Traverse Bay (Twenter et a1. 1985). In other areas the transition is similar from
lacustrine clay into beach and dune sands along most of the lowlands surrounding the
Grand Traverse Bay shore. I have interpreted the areas on the surﬁcial geology map
depicted with lacustrine sand and gravels to have this same transition from ﬁner to
coarser grained materials. Thus low areas (chosen from the Digital Elevation Model) that
are depicted as surﬁcial lacustrine sand and gravel likely have some amount of lacustrine

clays beneath.

8O

 

After the ice had left the Great Lakes basin completely (about 8,500 B.P.) the
level of water in the lakes varied due to the location of the lake outlet. Changes in the
level of the lakes have had a great effect on the types of sediments that were being
deposited in the GTBW. Sometimes the level of the lakes were 100 meters below the
lake levels today. This could cause changes in the relative base level of streams and
severe erosion of many of the deposits of the region. The Great Lakes basin has also
undergone isostatic rebound due to the weight of the glaciers being removed from the
system. This has caused a shift in the elevations of many of the beaches from previous
lakes in the basin. What hasn’t changed is the formation of sand dunes along the current
shore as well as some the high areas of the watersheds. Along the west arm of Grand

Traverse Bay and the east shore of the east arm dunes are quite common.

81

 

HYDROSTRATIGRAPHIC INTERPRETATION

The ﬁrst approach that was used to develop the glacial stratigraphy and
hydrostratigraphy of the GTBW focused on correlating a regional clay (till). We thought
if we could ﬁnd and correlate a regional clay, that this clay would act as the bottom of
our water table aquifer. After many weeks of working with Microsoft Access and
Arcview GIS trying to correlate water well drillers logs with the word clay, another
approach was taken. This approach entailed obtaining a digitized map of surﬁcial
sediments and assigning zones of hydraulic conductivities. From that point on focus was
shifted to creating a bedrock elevation for the GTBW and adding complexity in the form
of conductivity zones.

The bedrock that outcrops or subcrops beneath the glacial sediments of the
GTBW is mainly shales with some limestone and a very small amount of sandstone. This
provides us with a very good bottom of the aquifer. Shales of the Michigan basin have
over the years proven to be very good barriers to vertical ﬂow (Sibley, Personal Comm.)
The Ellsworth and Antrim shales reach thicknesses of up to 300 meters. The ﬁrst major
break came with the discovery of a l-kilometer scale digitized map of Farrand and Bell’s
Surﬁcial Geology Map of Michigan. This map gave us the framework for building a
discrete number of hydraulic conductivity zones based on the sediment characteristics
within those zones. We decided to use 6 zones for the GTBW; in no particular order they
are; End Moraines - Coarse, Glacial Outwash, Glacial Till - medium, Glacial Till —
coarse, Lacustrine Sand and Gravel, and Dune Sand. These zones were chosen to be
homogenous in the horizontal direction with respect to the model. Even though we know

this is truly not the case and even less the case in the glacial setting, this provided a

82

reasonable estimate of regional geology. Next we ground truthed these zones by
comparing them as best as we could to the water well logs and oil/gas well logs by
lithology. Obviously “Glacial Till medium and coarse” (refer to Figure 24) we could not
tell apart in the logs but we could see if there were identiﬁers in the drill logs such as
stony or pebbly. The glacial till - medium and coarse is the same as the ground moraine
as discussed in previous sections. The northward extent of the Manistee is also digitized
as glacial till.

After ground truthing the logs and the maps, which worked well, the next task
was to assign properties to each of the characteristic zones. In the ﬁrst model that was
developed the model was only 1 layer thick, which meant we could only specify one
conductivity value for the whole 100 by lOO-meter cell. There were numerous factors to
take into account such as clay layers not accounted for in the zones. We accounted for
various heterogeneities by adjusting the vertical conductivity in zones that may be more
acceptable to having a large continuous clay lenses. Table 3 shows the conductivity
values used for the ﬁrst model developed for the GTBW. This model proved to provide a
reasonable estimate of hydraulic heads for the GTBW (Figure 6).

A code was developed that split the one layer model into two layers and smoothed
the conductivity data for the outwash and moraine deposits. By smoothing, there was not
one value of hydraulic conductivity for the whole zone but a transition from a high
conductivity to a low conductivity much like what would exist in an actual outwash
setting (Figure 29). Proximal to the head of outwash a higher energy environment exists
and as you move distally the energy becomes less where there is ﬁner grained sized

material and hence a lower conductivity. The till and moraine zones were treated as

83

 

 

F ,F , F £2005.» __F .2090
m m m 00.600 -. _,_F .5030
N 9 3 09000 - 0050.02 05
m - em - vw 592:0 3.00.0
-Foono - 1 Sb l E - l .996 0:60:5033
50.0 3.0 or 0:3 0:30

 

w 33:26:00 .0033? 33:30:00 23:02.0: N 00:.— §>=o:0:0o 23:30.0: w .30..

0.8.?— ﬁon :0.“ 8030030000 30:08, 0:: 3:80:03 .3008 30: 5 00m: 002? 3300:0000 "m 030,—.

0:0N EoEEom
.3008 30c 03030038.: 05 E 08:

84

 

0.0 2.5 5.0 7.5 10.0 12.5 15.0

Figure 29: Smoothed hydraulic conductivity values. The hydraulic
conductivity values for each sediment zone were smoothed using a
interpolation scheme to better represent actual sedimentologic conditions.

85

having the same conductivities in both layers of the 2 layer model. This is because there
is no justiﬁcation to split these units into two distinct layers. But what we did was
incorporate vertical conductivities in these zones by adjusting the leakance, which
accounts for the anistropic behavior or sediments and rocks, factor between the two
layers to account for the lower conductivities in the vertical direction in these deposits.
The lacustrine sands and gravels as well as the dune sand zones were believed to
have a lower conductivity material beneath them, a thick clay probably lacustrine in
origin due to the lack of pebbles and stones. All of the water well logs in these areas also
showed relatively thick (sometimes up to 300 feet) clay, or less likely 3 till, beneath the
surﬁcial sands and gravels. This is integrated into the hydrostratigraphy by these units
having higher hydraulic conductivities in the upper third of their thickness. This was
checked with well logs in these cases and shows a very good correlation. The lower
lacustrine clay units had very low conductivities, which transmit very little groundwater
through these regions. This acts as the bottom of the aquifer and closely correlates with
what other researchers (Twenter, 1985; Cummings et al., 1990) have found in well logs.
These researchers only looked at small portions of the watershed and were not concerned
with the big picture. These initial reports of thick extensive clay, and the ﬁeld
observation of the Ellsworth Shale (which we thought was an outcropping of clay) led us

to think that it continued across the whole GTBW.

86

CONCLUSION

From these ideas and mistakes we have learned to take a step back from the ﬁne
detail and explore other ways of incorporating geology into groundwater models. We
realize that we cannot incorporate all of the spatial variability of glacial deposits in our
model, but we can break down the system into small manageable zones that will provide
us with the best possible deﬁnition of the system in order to answer the questions at hand.
If we were interested in modeling the transport of a contaminant plume at a small scale
then we would be worried about the occasional clay or sand lenses that might provide a
barrier or pathway to ﬂow. But since we are modeling at a large enough scale, these
heterogeneities will tend to average out. Thus it is justiﬁable to give one large zone of
glacial sediments a single conductivity value if the results are interpreted in the context of
your initial assumptions about the geology. In conclusion, the development of the
hydrostratigraphy of the GTBW proved to be less complicated than initially thought

when developed in the context of scales.

87

APPENDIX B: GENERATION AND EXAMPLES OF TOTAL AND

DIFFERENTIAL SOURCESHEDS

INTRODUCTION

A sourceshed is deﬁned as the total area that contributes to a selected drainage
point. Dr. Bryan Pijanowski ﬁrst coined the term “sourceshed” when he was explaining a
method on how we could visualize chemistry data taken in the ﬁeld. A sourceshed is
simply a watershed for any point within a drainage basin. Sourcesheds can be used for a
variety of purposes such as aggregating data (e. g. land-use percentage), visualizing
chemical data between two water sampling stations, and analyzing the inﬂuence of land-
based processes in a catchment context. Sourcesheds can be developed with respect to
surface water ﬂow over landscapes but can be easily transferred to watertable elevations
or piezometric surfaces. Two types of sourcesheds have been developed: the total
sourceshed and the differential sourceshed. The total sourceshed is the one deﬁned above
(Figure 30), which includes the total drainage area to a selected surface water sampling
location. The differential sourceshed (Figure 31) is the area between two chosen
drainage points of differing elevation that lie in the same drainage network. Each type of
sourceshed has its strengths and weaknesses that will be discussed later in this section.
The following section will include an in-depth discussion on the development of
sourceshed from theoretical and applied viewpoints. The sections will also discuss
sourceshed pitfalls and sourceshed applications. Throughout the rest of the section it will

be assumed that we are talking about surface water sourcesheds unless otherwise noted.

88

.050: 000608 a 8 8:55:00 3.5 808 :88 20
maze—:28 3 08088» 80 0002803 2: 00.0 00008058 000,—. 00008058 50... 5m charm

 

gas 3 Gene 0800 Ognm gm Osmt 080' OSmn 080». gm" OSON

 

 

928.32 w v n v

 

 

 

 

11951130
W561!

 

 

 

 

 

 

14W
;
Omit

 

 

1405000
4
W509}

 

 

 

14‘0“!)
0W0|63

 

 

OMS”!

zsisooo

 

 

 

d5
€
‘

 

1425M

 

 

 

 

 

 

 

 

 

 

 

89

 

 

_ 20900 40900

Differential Sourcesheds
for the GTBW N

Goundwwater Model Boundary
Sou’ce Points

:1
3|

 

0 9 18Km

 

 

26m

2049“”

24m

 

 

 

$1 + + + + «g
.. . . . .

20000 40000 60000 80000

 

 

Figure 31: Differential sourcesheds. The differential sourcesheds shown
here are developed by calculating the drainage area between selected points,
shown in red, that lie within the same drainage basin.

90

 

TOOLS AND PREPROCESSING OF DATA

Sourcesheds were initially developed for visualization purposes but have many
other purposes. The development began as an attempt to deﬁne polygons that
represented drainage areas for our surface water sampling locations. These polygons
were then to be associated with chemistry data at a downgradient sample site in an effort
to visualize spatially variable data. Data and tools needed to develop a sourceshed for a
sampling location are:
1. Arc/Info with GRID module
2. High resolution Digital Elevation Model (DEM) for region of interest (30 meters

or less)

3. Relative location of drainage or sampling point
4. ArcView 3.0+

Arc/Info 7.2.1 (ESRI, 1998a) was used throughout the generation of the
sourcesheds. At the time of writing, Arc/Info 8.0.1 is supposed to be released soon, so
any commands that don’t appear in the new version of Arc/Info refer back to the 7.2.1
version. Also included at the end of this section are lists of Arc/Info commands and
functions used in the sourceshed generation. ArcView 3.1 was only used for
visualization and in further processing of the sourcesheds after initial generation. It
therefore is not absolutely necessary to have, but it is convenient and time saving.

It is important to have an accurate DEM on which to perform the sourceshed
analysis on. Any gaps or errors in the DEM will cause severe problems when calculating
the ﬂow patterns through the drainage basin. This will in turn cause elevation values of

drainage points to be inaccurate with respect to surface water drainage divides and will

91

result in inaccurate catchment delineation. DEMs can be downloaded for free from the
USGS Geo Data website (http://edcwww.cr.usgs.gov/doc/edchome/ndcdb/ndcdb.html).
The lO-meter and 30-meter DEMs are available by 124,000 quadrangles only, referred to
as 7.5 minute DEM on the Geo Data website. In the spring of 1998 the USGS changed
its data storage standard and went to a storage type referred to as SDTS (Spatial Data
Transfer Standard). SDTS was such a new method of storing data that conversions in
Arc/Info 7.2.1 did not even exist. Arc/Info 8.0.1 will have an SDTS to GRID ﬁle,
probably sdtsgrid in the grid module, format command that will make the DEM
generation process a lot easier. For the example in this section, each 1:24,000 quadrangle
was downloaded in SDTS and converted into an Arc/Info ASCII grid format using a C
program written by Sol Katz. The ASCII grid ﬁles were imported into Arc/Info grids
using the asciigrid command in Arc/Info. About 1 out of every 4 USGS SDTS ﬁles had
the wrong coordinates in them so time had to be taken to make sure that the ﬁles have the
correct information and distribution. After that, quadrangles were aggregated together
using the merge command in the Arc/Info grid module. But before using the merge
command make sure that no gaps will occur between DEMs, as the USGS DEMs won’t
match up exactly with each other. You may have to use the gridshift command to move
grids with respect to each other. Any missing data between grids must be ﬁlled in,
possibly by using the resample command in the Arc/Info grid module. Data that falls on
a border outside of a DEM that has no other DEMs next to it won’t be a problem. Not all
of the DEMs, a total of 48, were lO-meter by lO-meter, some were 30-meter by 30-meter.
To get around this the 10-meter DEMs were all resampled to 30-meters before using the

merge command. Make sure that all of the DEMs download has the same units for the

92

vertical datum. All of the 10-meter DEMs were in meters and all of the 30-meter DEMs
were in feet. Alter merging the grids the last step in processing the elevation data is
reprojecting. Arc/Info’s projection command can be used to accomplish this task. The
process from downloading to merging 48 1:24,000 quadrangles took about one and a half

weeks of solid work.

 

The location of the drainage or sampling point, as you will read later, will change F
with the generation of the drainage pattern for your catchment. Since DEMs are spatially
averaged small errors occur in the exact location of smaller streams and rivers. Thus it is
necessary to move locations of sample points to the nearest stream cell based on the r

DEM.

93

METHODS

The next few paragraphs will discuss the direct generation of sourcesheds for
multiple points on the topography of the GTBW. The generation of differential
sourcesheds is much easier to do than total sourcesheds. This lies in the fact that when
making differential sourcesheds the area between two points is of concern. Whereas in
total sourceshed generation ﬁle preparation is much more difﬁcult, which lies in the fact
that you are concerned with total areas that overlap each, which Arc/Info polygons does
not support. It takes one source ﬁle for all of the sample points in the differential
sourceshed generation and it takes one ﬁle for each one of the sample points in the total
sourceshed generation. This makes any analysis procedure with total sourcesheds very
time consuming for any more than a few drainage points. A method will be discussed
later to generate data for total sourcesheds from data generated with differential
sourcesheds using a simple spreadsheet technique.

The ﬁrst step in the generation process is to identify drainage patterns in your
DEM. It is necessary to have a digitized stream network either from topographic
quadrangles or Digital Line Graphs (DLGs) on hand at this point to compare with the
delineated streams. Figure 32 shows streams as digitized from air photos (blue) and
streams delineated from the DEM (green) for the GTBW. Differences, as discussed
above, lie in the spatial averaging of the DEM. Note the red point and the yellow point
on Figure 30. The red point refers to the actual stream sampling location as determined
by GPS measurements. The yellow point refers to the location where the same point lies

on the DEM generated drainage network. If a sourceshed technique were used on the red

94

 

800080—0— md :05 0:08 .3 .mo 0: 0800 0000 80.00.20 030 088 :0 000000—
0080: 005 0002 .803: 00000800 Ema 08 0:0 080 08 8 00m: 802.. 000280 05 80300:
80:08.80 00.88 08 800.002: 088.: 83. 80>: 00:80:0w Ema 0:0 mag: "an 0.88%

 

 

oooim

oooim

 

 

Ev. no 0.0 o 0.0
>um 0326:. 0:80 \<

View 08:00 9a.! ...\/\
802m 03000.00 .00
2.550004 oBEow 002.00 wao 0
80:80.58 02 03.3 80:33 29:3

0.00%...523

080M

 

 

 

0000'»:

2"

 

 

 

 

95

point it would yield no drainage area because it lies off of the drainage path in the basin.
In turn the yellow point yields a sourceshed depicted as the shaded area on the ﬁgure,
which can be attributed to the real drainage area of the red point. It is very important that
sampling locations are moved to the same drainage spot as depicted on the DEM
delineated drainage pattern. In some cases where river or stream channels have valley
widths less than the DEM cell size, sourceshed generation will not work for those rivers.
This should not occur with natural stream systems but may occur with man made
channels or human inﬂuenced drainage patterns (geologically new channels). As with
the case in the GTBW many rivers are dammed where ﬂow is restricted but still occurs
through the river. This appeared not to affect sourceshed generation. Natural lakes, such
as Elk Lake and Torch Lake (see Figure 1), have natural ﬂow through them but the
surface elevations are identical on the DEM. This caused problems because when
ﬂowdirections were developed no gradient existed, therefore no ﬂow paths were
developed. This obviously is not the real case, but occurred because the lakes are so
large and in hydraulic connection with each other. The lake bottoms actually slope
toward the drainage points, much like a river, but the lakes are so deep and have so much
water they don’t show a surface gradient at the resolution of a DEM. The way to get
around this would be to obtain baythmetric data for the lakes and generate ﬂow paths and
sourcesheds for the points using that data.

In order to generate river or drainage networks with a DEM, small inconsistencies
and depressions must be ﬁlled to avoid all ﬂow paths from going towards a single point.
Each step of sourceshed development will generate a new ﬁle. Filling can be

accomplished using the ﬁll command in the grid module of Arc/Info with the sinks option

96

selected. This will generate a DEM with sinks ﬁlled to a ﬂat level. The next step is to
generate a ﬂowdirection grid with the newly ﬁlled DEM. The ﬂowdirection command in
the Arc/Info grid module is a cellular command that assigns each cell with a ﬂow
direction. It assigns the cell by examining the elevations of the 8 cells surrounding the
central cell, where a direction towards the lowest cell will be assigned to the central cell.
The ﬂowdirection command provides input into the ﬂowaccumulation command. The
ﬂowaccumulation command calculates where ﬂow will accumulate with the information
obtained from the ﬂowdirection command. The command will follow ﬂowdirection
values until they leave the edge of the grid or reach a ﬂat plane (such as Grand Traverse
Bay). Values in the ﬂowaccumulation grid refer to the number of cells ﬂowing into that
cell. A value of 1000 in the ﬂowaccumulation grid designates that 1000 cells are ﬂowing
into that speciﬁc cell. The ﬂowaccumulation grid serves as a drainage network for the
whole DEM area if it was ﬂooded with water. This grid will ultimately become your
drainage network. In order to adjust the ﬂowaccumulation grid to represent a realistic
drainage network use the can function in the grid module. The con Boolean function
separates all values above or below, depending on > or < sign, a chosen
ﬂowaccumulation value. A good value to start with is 100, but a value should be chosen
to approximate the actual spatial distribution of streams and rivers. A value of 400 was
the value that visually looked the best in the GTBW. If a value of 100 is used in the can
function it will generate a grid that has cells labeled with a 1 that have 100 or greater cells
ﬂowing into them all other cells will designated a O. This ﬁle will be used as a template
for further drainage network generation. The gridline command can be used to generate

this ﬁle into a line coverage for easier viewing. The most important thing about this step

97

is that the value that is chosen has no bearing on sourceshed generation. Remember that
the ﬂowaccumulation command is not used in the ﬁnal sourceshed generation. In this
case the value chosen is only for visualization purposes. When generating a river
network for input into a model it is necessary to calibrate the ﬂowaccumulation value.

Now that a drainage network is generated for the DEM, locations of where sample
points exist are needed to generate sourcesheds. The easiest way to do this is to pick the
x-y coordinate off of the generated drainage network for each point. Make sure to pick
the coordinate off of the middle of the cell of the drainage network. Any mistakes such
as picking a cell that is not part of the main drainage network will be readily apparent in
the size of the sourcesheds that are generated. The x-y coordinate will be used as input
into the selectpoint command in the grid module. Another way to input the points is to
generate an Arc/Info point ﬁle from the x-y coordinates. The selectpoint function takes a
ﬁlled DEM and sampling locations to generate a “source” grid that has elevation values
only for cells that sampling locations fall within. The source grid will look empty but
make sure that there are elevation values in cells where sampling points are located. The
source grid will be used as direct input into the watershed function.

With the newly generated source grid as well as the ﬂowdirection grid you are
now ready to create sourcesheds for your drainage points. The watershed function will
be used to create the areas of your sourceshed. This function calculates the contributing
area above your selected drainage point, or if another point lies upgradient of the ﬁrst
point in the drainage basin, then the area above the ﬁrst but below the second higher
point will be calculated. Areas are output to a grid ﬁle with the sourceshed for a point

designated as the elevation of that point. There is no way to have the watershed function

98

only generate differential or total sourcesheds. If any other points lie in the same
drainage basin then differential sourcesheds will be generated if not total sourcesheds are
generated. The grid of sourcesheds needs to be converted into polygons for further type
of analysis. The gridpoly function will take a grid and convert each unique group of grid
values into its own polygon. Be careful because if two drainage points are selected that
have the same elevation, which also abut each other, they could be combined into one
polygon. This step may need to be done many times to generate the correct sourcesheds
for each point, modifying your x-y values for drainage points each time.

As was stated above calculating total sourcesheds for more than a few sites is a
difﬁcult task. It is also very difﬁcult to visualize any information with total sourcesheds
because when dealing with multiple sites their drainage areas may overlap each other
thus obscuring information. Total sourcesheds tend to be better suited for data
aggregating and analysis. In a situation where differential sourcesheds are already
generated it is possible to calculate total sourcesheds from them very easily with a
spreadsheet. The upper right hand quarter of Figure 33 shows an enlargement of Mitchell
Creek with 3 sourcesheds for sample site 9, 65, and 63. The differential sourceshed for
sample 9 is the grey area in between the two colored areas. The sourceshed for sample 65
is yellow and for 63 its green. The total sourceshed for site 9 is the thick black line that
borders all of Mitchell Creek, or essentially the watershed for Mitchell Creek. If data that
is collected for site 65 and 63 is added to the site 9 differential sourceshed the data output
will be for the site 9 total sourceshed. This is a method that is simple to use with the aid
of a spreadsheet and can be done for most situations without redoing the total sourceshed

generation for each site.

99

.0 000 000 0000000000 0000
0 2000 00 0 000 00.0 0000000000 0000000000000 000 003 00:50:00 00 :00 00000000000 0000000000
3 000 0:0 8 000 00,—. 000000000 3000000000 00000 0000—00—00 0000000000 0000.0. “mm 000000

00 0:0 20500

00 000 0.0500 0 0:0 0.00.00

0000”.
00.02% _\/\
8% 0.0.08 0

 

              

 

 

 

 

 

 

 

 

 

 

mo 35 0_QEmw
000 00000058

0.00.0 .6500.

no 95 oEEmm
000 00000053

 

 

 

100

SOURCESHED APPLICATIONS

There are countless ways to utilize the extreme versatility of sourcesheds. First of
all, sourcesheds do not have to be generated for surface water studies. Groundwater
sourcesheds have been developed using simulated heads from a regional groundwater
ﬂow model. In this unconﬁned aquifer we were dealing with a system that is strongly
inﬂuenced by surface water discharges. Thus it was easy to develop sourcesheds in this
case because major ﬂow patterns were similar to those of the surface water drainage
network. Secondly, sourcesheds don’t have to just analyze or visualize the results of ﬁeld
sampling for stream chemistry, they can show where to sample for stream chemistry.
Sourcesheds are an extremely powerful tool in this respect because they have the
possibility to predict where and when a chemical of interest might be found if a known
source point existed. Therefore, the sourceshed will help deﬁne a better sampling
strategy to detect the occurrence of the chemical. Finally sourcesheds can be linked to
other tools to provide an even more powerful tool. For example if a time factor is added
to the sourceshed analysis sourcesheds will become dependent upon time. You can have
a 10 year, 20year, 30 year sourceshed for relatively slow systems such as groundwater or
you can have 2 hour, 6 hour, 10 hour sourceshed for faster systems such as a surface
water system undergoing a precipitation event.

Output from the high-resolution groundwater model can be used to develop
spatial sourcesheds for each of potential stream sample sites within the GTBW. Potential
stream sample locations can be determined by developing a spatial layer composed of
river/stream and road crossings at even interval spacing along river nodes. Two types of

sourcesheds could be developed with the land use and groundwater model data, “static”

lOl

 

 

 

groundwater sourcesheds and time-varying sourcesheds, thus enabling the deﬁnition of
regions that drain into stream segments for increments such as 1 year, 5 years, 10 years,
and 20 years.

Surface and groundwater sourcesheds can be used to assess the watershed’s
ecosystem health and studied spatially using ArcView. Each sourceshed (surface and
groundwater are evaluated independently) is assigned the value of the absolute
concentration of a parameter or concentration ratio of the water sample representing the
sourceshed. Both types of sourcesheds are studied because at times during the year the
dominant source for stream water changes (e. g., surface and groundwater for high ﬂow
and low ﬂow, respectively).

The sourcesheds also provide for modeling the relation of land use to
biogeochemical parameters. For example, each sourceshed can be characterized in terms
of associated land uses. Results of this analysis could be used to determine if types of
land use can be discriminated from one another in terms of the parameters.
Interpretations regarding the biogeochemical controls on the data can then be inferred. In
addition, the relative importance of the factors in describing the data can be used as an

indication of the relative signiﬁcance of each in controlling systems biogeochemistry.

102

 

 

CONCLUSIONS

Sourcesheds are powerful tools that allow researchers to spatially aggregate data
into regions that reﬂect properties of the watershed of interest. Sourcesheds are more
than just polygons, they are the area that contributes to a selected sample location.
Therefore they provide researchers with the ability to visualize and analyze point or
regional data within the context of possible source region. Overall, sourcesheds reﬂect
the complexity and simplicity of the emerging tools available to researchers working in

watershed hydrology.

103

ARC/INFO FUNCTION DEFINITIONS

FILL <in_grid> <out_grid> {SINK | PEAK} {z_limit} {out_dir_grid}

Arguments

<in_grid> - a grid representing a continuous surface.

<out_grid> - the output grid which has had its sinks ﬁlled or peaks cut.

{SINK | PEAK} - keyword indicating whether sinks will be ﬁlled or peaks will be cut.

SINK - all sinks which are less than the {z_limit} lower than their lowest adjacent neighbor will be ﬁlled
to the height of their pour point. A sink is a cell with undeﬁned drainage direction; all cells around it are
higher. The pour point is the boundary cell with the lowest elevation for the contributing area of a sink. If
the sink were ﬁlled with water, this is the point where water would pour out.

PEAK - all peaks less than the {z_limit} higher than their highest adjacent neighbor will be cut down to the
height of that neighbor. A peak is a cell that is higher than all adjacent cells.

{z_limit} - the maximum difference between a sink and its pour point or a peak and its highest adjacent
neighbor. If the difference in 2 value between a sink and its pour point is greater than {z_limit}, that sink
will not be ﬁlled. If the difference in 2 value between a peak and its highest adjacent neighbor is greater
than {z_limit}, that peak will not be removed. The default is to ﬁll all sinks (or remove all peaks)
regardless of depth.

{out_dir_grid} - an optional output grid showing the ﬂow direction from each cell of the <out_grid> to the
steepest downslope neighbor. The output is the same as that of the FLOWDIRECTION function.

Example Expression:
ﬁll DEM DEM_ﬁll sink 100 ﬂowdir

FLOWDIRECTION (<surface_grid>, {o_drop_grid},
{NORMAL | FORCE})

Arguments

<surface_grid> - a grid representing an elevation surface.

{o_drop_grid} - an optional output grid showing a ratio of the maximum change in elevation from each
cell along the direction of ﬂow, to the path length between centers of cells, expressed in percents.

{NORMAL | FORCE} - controls the direction of ﬂow at the edges of the surface grid.

NORMAL - if the maximum drop on the inside of an edge cell is greater than zero, the ﬂow direction will
be determined as usual, otherwise the ﬂow direction will be toward the edge. Cells that should ﬂow from
the edge of the surface grid inward will do so.

FORCE - all cells at the edge of the surface grid will ﬂow outward from the surface grid.

Example Expression:
ﬂowdir = ﬂowdirection (DEM)

FLOWACCUMULATION (<dir dgrid>, {weight _grid})

Arguments

<dir_grid> - a grid showing direction of ﬂow out of each cell. This can be created using the
FLOWDIRECTION function.

{weight _grid} - the weight to be assigned to each cell. If no {weight_grid} is speciﬁed, a default weight of
one will be applied to each cell. For each cell in the output grid, the result will be the number of cells that
ﬂow into it.

Example expression:
ﬂowacc = ﬂowaccumulation (ﬂowdir)

104

 

SELECTPOINT (<grid>, <x, y>, {INSIDE | OUTSIDE})
SELECTPOINT (<grid>, <point_ﬁle>, {INSIDE | OUTSIDE}
SELECTPOINT (<grid>, <point_cov>, {INSIDE | OUTSIDE})

SELECTPOINT (<grid>, <*>, {INSIDE | OUTSIDE»

Arguments

<grid> - an input integer or ﬂoating-point grid. The <grid> identiﬁes the cell values which will be selected
inside or outside a speciﬁed point. The grid can also be created in an expression.

<x, y> - identiﬁes the coordinates of the point to be selected. The coordinates are speciﬁed in map units
and are in the same units as the <grid>.

{INSIDE | OUTSIDE} - identiﬁes which cells should be selected, those contained in, or external to, the
cell containing the selected point.

INSIDE - a keyword specifying that the cell in which the selected point falls will be written to the output
grid. All cells outside the box will receive NODATA on the output grid.

OUTSIDE - a keyword specifying that the cells outside the input box should be selected and written to the
output grid. The cell within which the selected point falls will receive NODATA on the output grid.
<point_file> - the name of an ASCII text ﬁle containing coordinates of points to be selected.
<point_cov> - the name of the input coverage with point topology. The locations and attributes of the
points from this coverage will be transferred to the respective cells in the output grid.

<*> - allows for the interactive graphical input of the input point which will be used to identify the
selection criteria. The graphics environment (i.e., the display window and the map extent) must be set.
The grid speciﬁed by <grid> should be displayed.

Example Expression:
outgrid = selectpoint (DEM_ﬁll, sample _points, INSIDE)

WATERSHED (<dir_grid>, <source_grid>)

Arguments

<dir_grid> - a grid showing direction of ﬂow out of each cell. This can be created using the

F LOWDIRECTION function.

<source_grid> - a grid representing cells above which the contributing area, or catchment, will be
determined. All cells that are not NODATA will be used as source cells.

Example expression:
wsheds = watershed (ﬂowdir, source)

References for Arc/Info Functions:

Greenlee D. D. 1987. Raster and Vector Processing for Scanned Linework, Photogrammetric Engineering
and Remote Sensing. Vol. 53, No. 10, October 1987, pp. 1383-1387.

Jenson S. K. and J. O. Domingue. 1988. Extracting Topographic Structure from Digital Elevation Data for
Geographic Information System Analysis, Photogrammetric Engineering and Remote Sensing. Vol. 54, No.
11, November 1988, pp. 1593-1600.

Tarboton D. G., R. L. Bras, I. Rodriguez-Iturbe. 1991. On the Extraction of Channel Networks from Digital

Elevation Data, Hydrological Processes. Vol. 5, 81-100.

105

 

 

REFERENCES

106

 

REFERENCES

Anderson, MP, and W. W. Woessner, 1992. Applied Groundwater Modeling -
Simulation of ﬂow and advective transport. Academic Press, San Diego, 381 pp.

Bethke, C.M., 1986. Inverse hydrologic analysis of the distribution and origin of Gulf
Coast-type geopressured zones. Journal of Geophysical Research (91): 6,535-
6,545.

Beven, K]. and M.J.Kirkby, 1979. A physically based, variable contributing area model
of basin hydrology. Hydrological Sciences Bulletin, 24: 43-69.

Belijaars, A.C.M., and F .C. Boseveld, 1997. Caubauw data for the validation of land
surface parameterization schemes. J. Climatology, 10: 1172-1193.

Blewett, W.L., and Winters, HA, 1995. The importance of glaciofuvial features within
Michigan’s Port Huron Moraine. Annals of the Association of American
Geographers, 85 (2): 306-319.

Bosch, J .M., and J .D. Hewlett, 1982. A review of catchment experiments to determine
the effect of vegetation changes on water yield and evapotranspiration. J.
Hydrology, 55: 3-23.

Brater, ER, 1972. A hydrological model for estimating the inﬂows to and outﬂows from
Grand Traverse Bay. 32, Michigan Sea Grant Program, Ann Arbor.

Bricker, S.M., and Lilienthal, R.T., 1976, Grand Traverse County bedrock topography:
Michigan Geological Survey Division, Map 377ZBT.

Brigham Young University, 1994. D.O.D. Groundwater Modeling System (GMS), Salt
Lake City, UT.

Canter, L.W., 1996. Nitrates in Groundwater. CRC Press, New York, 263 pp.

Cummings, T.R., J .L. Gillespie, and N. G. Grannemann. 1990. Hydrology and land use
in Grand Traverse County, Michigan. U.S.G.S. Water-Resources Investigations
Report 90-4122, Lansing, MI.

Doe, W.M., Saghaﬁan, B., and Pierre, Y.J. (1996). Land-use impact on watershed

response: the integration of two-dimensional hydrological modeling and
geographical information systems. Hydro. Proc. 10: 1503-1511.

107

 

Dorr, J .A., D.F. Eschman, 1970. Geology of Michigan. The University of Michigan
Press, Ann Arbor, 476 pp.

Dunn, S. M., R. Mackay, 1995. Spatial variation in evapotranspiration and the inﬂuence
of land use on catchment hydrology. J. Hydrology, 171: 49-73.

Dunn, S.M., R. Mackay, R. Adams, D.R. Oglethorpe, 1996. The hydrological component
of the NELUP decision support system: an appraisal. J. Hydrology 177, 213-235.

Ehlers G.M., 1938, Eighth annual ﬁeld excursion, May 28-29, 1938: Michigan Academy
of Sciences, Arts and Letters, Section of Geology and Mineralogy, mimeographed
handout.

Ellis, B.G., 1982. Nitrate contamination of groundwater on the Old Mission Peninsula:
Contribution of land reshaping and septic drainﬁelds: Department of Crop and
Soil Sciences, Michigan State University, East Lansing, unpublished report, 17 p.

Elowski, RC, and Bricker, S.M., 1982, Kalkaska County bedrock topography:
Michigan Geological Survey Division, Map 3654BT, scale 1:62,SOO.

Eschman, D.F., William R. F arrand, and Edward B. Evenson, 1973. Pleistocene geology
of the northwestern quarter Southern Peninsula, Michigan, Geology and the
Environment; Man, Earth and Nature in Northwestern Lower Michigan; Michigan
Basin Geological Society Annual Field Excursion; Geological Papers, The Earth.
Mich. Basin Geol. Soc., Lansing, Michigan, pp. 19-25.

(ESRI)a, E.S.R.I., 1998. ARC/INFO. ESRI, Redlands, CA.
(ESRI)b, E.S.R.I., 1998. ArcView. ESRI, Redlands, CA.

Farrand, W., and Bell, D. 1982. Quartemary Geology of Southern Michigan. (1:500,000
scale map). Ann Arbor: University of Michigan.

F eick, G. et al., 1972. Release of Mercury from Contaminated Freshwater Sediments by
the Runoff of Road Deicing, Science, 175(3): 1142-43.

Field, C.K., Sover, RA, and Lott, A. (1996) estimating the effects of changing land use
patterns on Connecticut Lakes. J. Env. Qual. 25: 325-333

F lintrop, C., B. Hohlmann, T. Jasper, C. Korte, O. Podlaha, S. Scheele, and J. Veizer,
1996. Anatomy of pollutionzstreams of North Rhine- Wesphalia, Germany. Am. J.
of Science, 296(296): 58-98.

Freeze, R.A., and John. A. Cherry, 1979. Groundwater, 1. Prentice Hall, Englewood
Cliffs, 604 pp.

108

 

 

 

Fryer, W. M. 1982. An investigation of salt contaminated groundwater at the Shepard #1-
21/2-21 lease site Kalkaska Township, Kalkaska County. Michigan Geological
Survey Division OF R 82-2.

Garven, G., 1985. The role of regional ﬂuid ﬂow in the genesis of the Pine Point deposit,
western Canada basin. Economic Geology, 80: 307-324.

Gutschick, RC, and CA. Sandberg, 1991a. Upper Devonian biostratigraphy of
Michigan Basin. In: P.A. Catacosinos, and RA. Daniels Jr. (Editor), Early
sedimentary evolution of the Michigan Basin. Geological Society of America,
Denver, pp. 155-180.

Gutschick, RC, and CA. Sandberg, 1991b. Late Devonian history of Michigan Basin.
In: P.A. Catacosinos, and PA. Daniels Jr. (Editor), Early sedimentary evolution
of the Michigan Basin. Geological Society of America, Denver, pp. 181-202.

Herlihy, A.T., Stoddard, J .L., Johnson, CB, (1998) The relationship between stream
chemistry and watershed land cover data in the Mid-Atlantic Region, U.S. Water,
Air, Soil Pool. 105: 377-386.

Hibbert AR, 1971. Increases in streamﬂow after converting Chaparral to grass. Water
Resources Research. 7: 71-80.

Hookey, GR, 1987. Prediction of delays in groundwater response to catchment clearing.
In: A.J.a.D.R.W. Peck (Editor), Hydrology and Salinity in the Collie River Basin,
Western Australia. J. Hydrology, pp. 181-198.

Homberger, G.M, Raffensperger, J .P., Wiberg, PL, and K.N. Eshleman, 1998.
EleOments of Physical Hydrology. Johns Hopkins University Press, Baltimore.
302 pp.

Jones, A.L., and Bernard N. Sroka, 1997. Effects of highway deicing chemicals on
shallow unconsolidated aquifers in Ohio, interim report, 1988-93. 97-4027, US.
Geological Survey, Colombus, Ohio.

Kelly, R.W., 1968 (reprinted 1977), Bedrock of Michigan: Michigan Geological Survey
Division Geologic Map GM 1, scale 1:2,500,000.

Kesling, R. V., Segal, R.T., and Sorensen, H.O (1974). Devonian strata of Emmet and
Charlevoix Counties. Lansing, Michigan, Michigan Basin Geological Society.

Leverett, F ., and Taylor, F. 1915. The Pleistocene of Indiana and Michigan and the
history of the Great Lakes. United States Geological Survey Monograph 53.
Washington, DC: US. Government Printing Ofﬁce.

Macpherson, D.K., and Peck, A.J., 1987. Models of the effect of clearing on salt and
water export from a small catchment. In: A.J.P.a.D.R. Williamson (Editor),

109

 

Hydrology and Salinity in the Collie River Basin, Western. J. Hydrology, pp. 163-
179.

Mandle, R.J. and AL. Kontis., 1992. Simulation of regional ground-water ﬂow in the
Cambrian-Ordovician aquifer system in the northern Midwest, United States. 93-
0226, US. Dept. of the Interior, US. Geological Survey, Denver, CO.

Martin, H. M., 1957. Outline of the geologic History of the Grand Traverse region.
Michigan Department of Conservation, Geological Survey Division, 8 p., ﬁgs.

Martin, P.J., and E.O. Frind, 1998. Modeling a complex multi-aquifer system: The
Waterloo Moraine. Ground Water, 36(4): 679-690.

Mauser W., and S Schadlich, 1998. Modelling the spatial distribution of
evapotranspiration on different scales using remote sensing data. J. Hydrology,
212-213: 250—267.

McDonald, MG. and Harbaugh, AH, 1988. A modular three-dimensional ﬁnite-
difference groundwater ﬂow model. in Techniques of Water-Resources
Investigations of the United States Geological Survey, Book 6, Chapter A1.

Melhom, W.N., 1956. Valders glaciation of the Southern Peninsula of Michigan,. Ph.D.
Dissertation, Department of Geology, University of Michigan, Ann Arbor.

Millstein R.L, 1983, Charlevoix County bedrock topography: Michigan Geological
Survey Division, Map 3773BT, scale 1:62,500.

Modica, E., Thomas E. Reilly, and David W. Pollock, 1997. Patterns and age distribution
of ground-water ﬂow to streams. Ground Water, 35(3): 523-537.

Person, M.A., and Garven, G., 1994. A sensitivity study of the driving forces on ﬂuid
ﬂow during continental-rift basin evolution. Geological Society of America
Bulletin(106): 461-475.

Pijanowski, B. C., S. H. Gage, D.T. Long and WE. Cooper. 1998. A Land
Transformation Model for Michigan’s Saginaw Bay Watershed: An Application
of a Geographic Information System, in Landscape Ecology: A Top Down
Approach, eds. Larry D. Harris and James Sanderson, CRC Press/Lewis
Publishers.

Pijanowski, BC. 1997. Issues of space and time scales in changing landscapes: systems
approach and the use of geographic information systems. In Lecture Series in
Landscape Ecology. Larry Harris, ed.

Pijanowski, B. C., T. Machemer, S. Gage, D. Long, W. Cooper and T. Edens. 1996. The
use of a geographic information system to model land use change in the Saginaw

110

 

Bay Watershed. Proceedings of the Third International Conference on GIS and
Environmental Modeling, Sante Fe, New Mexico, January 21-26, 1996. GIS
World Publishers on CD-ROM.

Pijanowski, B. C., T. Machemer, S. Gage, D. Long, W. Cooper and T. Edens. 1995. A
land transformation model: integration of policy, socioeconomics and ecological
succession to examine pollution patterns in watershed. Report to the
Environmental Protection Agency, Research Triangle Park, North Carolina. 72

PP-

Pucci, Jr., A.A., and Pope, DA, 1995. Simulated effects of development on regional
ground-water/surface-water interactions in the northern Coastal Plain of New
Jersey. J. Hydrology, 167: 241-262.

Rajagopal, R., 1978. Impact of land use on ground water quality in the Grand Traverse
Bay region of Michigan. Journal of Environmental Quality, 7(1): 93-98.

Reska, CR. 1983, Antrim County bedrock topography: Michigan Geological Survey
Division, Map 3767BT, scale 1:62,500.

Roberts G., and A.M. Roberts, 1992. Computing the water balance of a small agricultural
catchment in southern England by consideration of different land-use types. 11.
Evaporative losses from different vegetation types. Agricultural Water
Management, 21: 155- 166.

Rothacher J ., 1970. Increase in water yield following clear-cut logging in the paciﬁc
northwest. Water Resources Research, 6:653-658

Scott RR, and RA. Sudmeyer, 1993. Evapotranspiration from agricultural plant
communities in the rainfall zone of the southwest of Western Australia. J.
Hydrology, 146: 301-319.

Sibley D.F., Personal Communication, July 5, 1999: East Lansing, MI.

Skillings, D. M., 1982. Report on groundwater contamination Mancelona and Cold
Springs Township, Michigan. Michigan Geological Survey Division OFR 82-3.

Sucksdorff, Y., and C. Ottle, 1990. Application of satellite remote sensing to estimate
areal evapotranspiration over a watershed. J. Hydrology, 121: 32-333.

Taniguchi, M., 1997. Subsurface water responses to land cover/use changes: An
Overview. In M. Taniguchi ed. “Subsurface hydrological responses to land cover
and land use changes” Kluwer Academic Publishers, 1997 226 pp.

Thorbum, P.J., Cowie, B.A., and Lawrence, PA, 1991. Effect of land development on

groundwater recharge determined from non-steady chloride proﬁles. J.
Hydrology, 124: 43-58.

111

 

Twenter F.R., T.R. Cummings, and N.G. Grannemann, 1985. Groundwater
Contamination in East Bay Township, Michigan. U.S.G.S. Water—Resources
Investigations Report 85-4064, Lansing, MI.

Vesterby, M and R. Heimlich. 1991. Land use and demographic change: results from
fast-growing counties. Land Economics 67(3):279-291.

Wiley, MJ, Kohler, S.L., and Seelbach PW. (1997) Reconciling landscape and local
views of aquatic communities: lessons from Michigan streams. Freshwater Biol.
37:133-148.

Wilson, J. P. 1996. 1996 GIS-based land surface/subsurface modeling: new potential for
new models? Proceedings of the Third International Conference on GIS and
Environmental Modeling, Sante Fe, New Mexico: GIS World Publishers on CD-
ROM.

Wilson, J .L., and J. Liu, 1995. Backward tracking to ﬁnd the source of the pollution. In
Waste management: From risk to reduction (Ch. 10), ed. R. Badha. Albuquerque,
New Mexico: ECM Press.

Zheng, (1992) MT3D: A modular three-dimensional transport model, S.S. Papadopulos
& Assoc, Inc.

112

      

MICHGA ST TE RSI

mum/mm in:

3 1293

L BRAHIES

{girl/ill” mm

36 5844