EVALUATING THE INFLUENCE OF HISTORICAL AND CURRENT
LAND COVER ON SEDIMENT AND NUTRIENT TRANSPORT
By
Christopher L. May

A THESIS
Submitted to
Michigan State University
in partial fulfillment of the requirements
for the degree of
MASTER OF SCIENCE
Geological Sciences
2011

ABSTRACT
EVALUATING THE INFLUENCE OF HISTORICAL LAND COVER
ON SEDIMENT AND NUTRIENT TRANSPORT
By
Christopher L. May
Alterations to historical land cover along with the spatial variability of
current land cover can directly impact groundwater flow and the water quality of
streams. As human population grows, land cover will continue to be altered,
enhancing its importance. This study examines how changing land cover affects
recharge to groundwater and groundwater transport to streams. The impact of
historical and current land cover on hydrology in the Jordan River Watershed
was modeled. The geometry of the sediment deposited within the Jordan River
was studied using Ground Penetrating Radar to correlate changes in
sedimentation to past land cover and climate change. Ground Penetrating Radar
combined with sediment cores aided in the interpretation of certain sediment
packages. The last part of the study examines total dissolved Phosphorous and
Nitrate groundwater inputs from land cover by modeling the transport of nutrients
through the Sycamore Creek Watershed. Nutrient loading for different land
covers and in-stream processing for Total Dissolved Phosphorous and Nitrate
were optimized to match measured stream data.

ACKOWLEDGEMENTS

I would like to thank my advisor, Dr. David Hyndman, for his help and
guidance through this process. I am also thankful for the insight and help that
have been provided by other two committee members, Drs. Remke van Dam and
Warren Wood. I am very grateful for the help and advice given by my colleague
Dr. Anthony Kendall. A thanks also goes out to Dushmantha Jayawickreme,
Abby Norton, Sherry Martin, and Mine Dogan for their help at various stages of
this study. This work was made possible by the funding from the Friends of the
Jordan River, National Science Foundation, and the Environmental Protection
Agency. A thanks also goes to Dr. Jan Stevenson, Seth Hunt, and Nick Welty for
their contributions to the data collection and sample analysis. I finally want to
thank my wife Sara for her continued support and inspiration.

iii

TABLE OF CONTENTS

LIST OF TABLES ................................................................................................ vi
LIST OF FIGURES ............................................................................................. vii
CHAPTER 1: EVALUATING THE IMPACTS OF LAND COVER CHANGES ON
STREAM FLOW BEHAVIOR................................................................................ 1
INTRODUCTION............................................................................................... 2
Site Description.............................................................................................. 6
METHODS ...................................................................................................... 16
Model Inputs ................................................................................................ 16
Steady-state Model ...................................................................................... 26
Transient Field Data..................................................................................... 32
Recharge Calculations for Land Cover ........................................................ 40
Transient Model ........................................................................................... 44
RESULTS AND DISCUSSION ........................................................................ 47
Model Results .............................................................................................. 47
Implications for Sediment Transport ............................................................ 49
CONCLUSIONS .............................................................................................. 52
BIBLIOGRAPHY.............................................................................................. 55
CHAPTER 2: USING RADAR STRATIGRAPHY TO EXAMINE THE RESPONSE
OF HISTORICAL SEDIMENT TRANSPORT TO LAND COVER
ALTERATIONS................................................................................................... 59
INTRODUCTION............................................................................................. 60
Site Description............................................................................................ 62
METHODS ...................................................................................................... 67
Ground Penetrating Radar........................................................................... 67
GPR Data Collection.................................................................................... 69

iv

GPR Processing .......................................................................................... 71
Auguring ...................................................................................................... 72
RESULTS AND DISCUSSION ........................................................................ 74
Borehole Data .............................................................................................. 74
GPR Signatures ........................................................................................... 76
GPR Interpretation ....................................................................................... 80
CONCLUSIONS .............................................................................................. 81
BIBLIOGRAPHY.............................................................................................. 83
CHAPTER 3: EVALUATING THE INFLUENCE OF LAND COVER ON
NUTRIENT TRANSPORT AND IN-STREAM NUTRIENT PROCESSING ......... 85
INTRODUCTION............................................................................................. 86
Site Description............................................................................................ 89
Source and Sinks......................................................................................... 91
METHODS ...................................................................................................... 93
GIS Data ...................................................................................................... 93
Steady-state Model ...................................................................................... 99
Groundwater transport using MT3D........................................................... 102
Optimization using a simple 1-D In-stream Processing Description........... 107
RESULTS AND DISCUSSION ...................................................................... 109
Optimization Results .................................................................................. 109
Final MT3D Model Results......................................................................... 110
CONCLUSIONS ............................................................................................ 118
Acknowledgments...................................................................................... 119
BIBLIOGRAPHY............................................................................................ 121

v

LIST OF TABLES
Table 1.1: Values for hydraulic conductivity after optimization for geologic zones
shown in Figure 8B. The starting hydraulic conductivity value for all the geologic
zones for the optimization was 5 m/d ................................................................. 30
Table 1.2: Total rain and runoff data collected by Onda on hill slopes near Point
Reyes, California immediately following a forest fire. This study used rainfall and
runoff data to calculate runoff ratios which was the basis for the precipitationrunoff relationship. Values in red indicate the collection tank was full and did not
actually record all the rain. Event 1 recorded 6.51 inches, but the actual rainfall
was 7.96 inches. Event 6 recorded 1.43 inches but the total rainfall was 2.46
inches. Thus, 7.96 and 2.46 inches were used in the runoff-precipitation
relationship. Modified from Onda (2008) ............................................................ 42
Table 2.1: Various materials and their relative permittivity and velocity, Modified
from Baker (1991)............................................................................................... 69
Table 3.1: Simulated flow and averaged measured flow at three
locations ........................................................................................................... 101

vi

LIST OF FIGURES
Figure 1.1: The Jordan River Watershed (JRW) is located in lower Northwest
Michigan (MIGDL) ................................................................................................ 7
Figure 1.2: Hydrograph, with precipitation, of the Jordan River at the United
States Geological Survey (USGS) Webster gage # 04127800. Flows and
precipitation plotted from 01/01/07-09/01/08 ........................................................ 8
Figure 1.3 A map of percent slope from the digital elevation model (A) for cells in
the Jordan River model area. Slopes in the Jordan River Valley (B) (shown by
the blue rectangle) contain the highest relief within the drainage basin.
Numbered locations refer points on the elevation profile for the Jordan River
(Figure 1.4) (Source, USGS SEAMLESS) .......................................................... 10
Figure 1.4: Elevation Profile along the Jordan River (blue line), relative to a
simple equilibrium profile (black line). The elevation profile was generated using
the points along the river that were extracted from the DEM (Figure 1.3A). The
red are locations marked on the river in Figure 1.3 ............................................ 11
Figure 1.5: Old railroad beds can be found in the Jordan River valley as linear
paths of bare ground .......................................................................................... 13
Figure 1.6: Agriculture as percent land area from 1900-1950 for the three
counties within the model area (See Figure 1.1 for county locations) (USDA).... 14
Figure 1.7: Land Cover Statistics (A) for the three different land cover data sets
and maps for pre-developed (B) and for current, 1978 MIRIS, (C) land covers
(MIGDL).............................................................................................................. 17
Figure 1.8: GIS maps for soils (A) and surface geology (B) (MIGDL)................. 20
Figure 1.9: Lake coverage of water, oil, and gas wells within the model
boundary shown in black (MIGDL). The red outline is the surface watershed.... 22
Figure 1.10: Drift thickness map (A) created from the interpolated bedrock
elevation (B) (MIGDL)......................................................................................... 24
Figure 1.11: Annual precipitation used to calculate recharge and runoff for the
Jordan River model area, includes snowfall (NCDC and MAWN) ...................... 26
Figure 1.12: Simulated map of heads (A), along with a comparison between
simulated vs. observed heads after PEST optimization (B)................................ 30

vii

Figure 1.13: Shows the JRW and groundwatershed (shaded) for the Jordan River
within the model area ......................................................................................... 32
Figure 1.14: Map showing gage locations within the Jordan River Watershed
(highlighted in red) (A) along with their stage-discharge relationships: USGS
Webster (B) Graves (C), Fisheries (D), and Deer (E) ......................................... 34
Figure 1.15: Spatial channel cross channel geometries for all the sites where
flows were conducted ......................................................................................... 37
Figure 1.16: Temporal cross channel geometry at Graves (A) and Deer Creek (B)
gage locations. See Figure 1.14A for gage locations ........................................ 38
Figure 1.17: Discharge at all four stream gages in response to Precipitation
Events. Arrows indicate lags (MAWN, NCDC, and USGS) ............................... 39
Figure 1.18: Precipitation: Runoff relationship used to calculate infiltration and
runoff. Modified from Onda (2008)..................................................................... 44
Figure 1.19: Approximate recharge for estimated for different land cover
scenarios within the model area ......................................................................... 45
Figure 1.20: Simulated Flow for Current Land cover within the Jordan River
model area plotted with stream flow and base flow for the USGS Webster gage
(USGS) ............................................................................................................... 48
Figure 1.21: Simulated stream flow at USGS Webster for all four tested
landscapes ......................................................................................................... 48
Figure 1.22: Hjulstrom’s diagram showing velocities required for the entrainment
and deposition of different sediment sizes. Modified from Prothero and Schwab
(2004) ................................................................................................................. 50
Figure 2.1: The Jordan River Drainage basin is in the northern portion of
Michigan’s Lower Peninsula (MIGDL) ................................................................ 62
Figure 2.2: Lake Level data for the NOAA gage for Lake Michigan and Lake
Huron. The dotted line shows the approximate level of Lake Charlevoix before
the opening of the channel between Lake Michigan and Lake Charlevoix
(NOAA). Lake levels are monthly averages........................................................ 64
Figure 2.3: Study area showing GPR lines, auger boreholes, and Vibracoring
locations adjacent to the Jordan River and Lake Charlevoix. The numbered
locations will be referred to in the text (MIGDL).................................................. 66

viii

Figure 2.4: GPR antenna setup during data acquisition in the winter of 2008.
Photo by Christopher L. May .............................................................................. 70
Figure 2.5: GPR processing algorithm. Arrows indicate the progressive
processing steps applied to the GPR data ......................................................... 71
Figure 2.6: Borehole data collected in the summer of 2007. All surface sediment
was a very thin (<1cm) layer of organics from current riparian vegetation and was
excluded from this graph. Boreholes are highlighted in Figure 2.3 .................... 75
Figure 2.7: A) Raw 2-D GPR data profiles shown with line numbers that
correspond to the GPR grid shown in Figure 2.3. Letters S, C, and M represent
sand, clay, and muck respectively, found in cores that coincided
with GPR lines .................................................................................................... 77
Figure 2.8: GPR profiles shown with 1 type of interpretation (light blue lines) .... 78
Figure 2.9: GPR profiles with interpretations of concave-up
reflectors (red lines)............................................................................................ 79
Figure 3.1: The model area along with Sycamore Creek Watershed (SCW,
highlighted in pink) is located in South Central Michigan (from MIGDL)............. 90
Figure 3.2: Map of Current Land cover for the Sycamore Creek model area
(MIGDL).............................................................................................................. 92
Figure 3.3: A 9 meter resolution DEM plotted with the model area (bold outline)
and the SCW (light outline) (SEAMLESS) .......................................................... 94
Figure 3.4: Surficial Quaternary Geologic Map used for geologic facies
description for MODFLOW and MT3D (Farrand and Bell, 1982) ........................ 95
Figure 3.5: Wells used to construct the map of the aquifer bottom, for water level
observations, and to provide septic well inputs. Wells used inside the model
boundary (yellow line) were used to simulate nutrient inputs from septic systems
(MIGDL).............................................................................................................. 96
Figure 3.6: Variogram model (using a Circular model, 255 m range, 164.14
nugget, 100 sill, 100 m lag size, 12 lags, and no anisotropy) (A) incorporated into
the interpolation of the bedrock elevation map (B) ............................................. 98
Figure 3.7: Simulated heads plotted relative to observed heads (A) and mapped
(B)..................................................................................................................... 101
Figure 3.8: A) Nitrate observations, B) TP observations (Welty, 2005) ............ 105

ix

Figure 3.9: Optimization results for Simulated and Observed total and differential
fluxes for Nitrate (A, B) and TDP (C, D) with in-stream processing
coefficients........................................................................................................ 111
Figure 3.10: Simulations based on the final optimized nutrient concentrations
from Agriculture, Pasture, Septic Tanks, and urban areas ............................... 113
Figure 3.11: Simulated and Observed TDP Flux from different land cover for
Sycamore and Mud Creeks .............................................................................. 114
Figure 3.12: Simulated and Observed Nitrate Flux from different land cover for
Sycamore and Mud Creeks .............................................................................. 115
Figure 3.13: Simulated and Observed TDP Concentration from different land
covers for Sycamore and Mud Creeks ............................................................. 116
Figure 3.14: Simulated and Observed Nitrate Concentrations from different land
cover for Sycamore and Mud Creeks ............................................................... 117

x

CHAPTER 1: EVALUATING THE IMPACTS OF LAND COVER
CHANGES ON STREAM FLOW BEHAVIOR

1

INTRODUCTION

Land cover alterations can affect recharge (Jayawickreme and Hyndman,
2007), stream flow (Mutie et al., 2006), as well as fluvial geomorphology and
stream ecology (Poff et al., 2006). It is important for management agencies to
understand the effect of land cover changes on the environment as the human
population grows. Alterations to the landscape such as logging, forest fires,
landslides, and farming can have significant impacts on stream flow and
sediment transport within a drainage basin. Significant changes in stream flow or
sediment transport can adversely affect stream ecosystems and degrade the
water quality, all of which can impact public water supply and recreation.
Examples of land cover alterations that can lead to changes in stream flow and
sediment transport include: forest fires, agriculture development, logging, and
landslides.
Studies have shown that logging can impact stream flow and sediment
transport thus leading to degraded water quality. Logging in the drainage basins
increased peak base flows, and overall sediment concentrations during high
flows (Waterloo et al, 2007). Cassie et al. (2002) found no changes to annual or
seasonal water yields in the middle reach or upper reaches after timber
harvesting of 2.3 and 23.4 % of the middle and upper reaches of the drainage
basin respectively, but peak flows increased in the upper tributary of the drainage
basin. Other studies (Hibbert, 1967; Bosch and Hewlett, 1982) showed that peak
flows, as well as annual and seasonal water yields only changed if greater than

2

20% of the total area was harvested. Other implications of logging include rising
water tables, and higher snow accumulation rates (Dube and Plamondon, 1995;
Winkler et al., 2005).
Landslides substantially change stream flow and sediment transport in a
drainage basin. Timber harvesting can increase the risk of landslides because
the absence of canopy cover leads to increased soil moisture and decay of roots,
which can lead to rapid soil saturation and the loss of soil strength (Johnson et
al., 2007). Other events that can cause landslides include earthquakes (Yang et
al., 2010) and heavy rain on slopes that contain only a thin layer of sediment over
steeply dipping rock (Saito et al., 2009). Landslides in turn are a major source of
sediment from steep valleys (Mikos et al., 2006). Kasai (2006) showed that
landslides can also change the fluvial geomorphology of a river system for long
periods of time by modifying the path of the river.
Forest fires can change channel morphology by decreasing bank stability,
physically changing the characteristics of the soils, and reducing
evapotranspiration (ET) (Dickman and Leefers, 2004). Batalla (2001) showed
forest fires can change stream flow behavior through large runoff events related
to hydrophobic soils created during the fire. Various studies have looked at
effects of forest fire on soils, runoff, erosion, and nutrient loss (Neary et al., 1999,
Cerda, 1998, Pierson et al., 2008, Onda et al. 2008, Tessler et al., 2008).
Because changes in stream flow due to land cover are being evaluated, the
biggest effect that a post-fire land cover would induce is a change in runoff.
Runoff is only considered in the post-fire land cover because this is one of the

3

few instances where runoff could occur within the JRW due to its course-textured
sediments. Water repellency (WR) by soils is heavily dependent on the
temperatures during a fire. Neary et al. (1999) discussed the expected
temperature ranges of fires occurring within different land covers (i.e. 200-300° C
for forests, 300-700° C for shrubs). The two types of wildfires that likely occurred
within the JRW would be slash fires of material left behind after tree harvesting,
which can have temperature ranges between 500-700 C, and forest fires in predeveloped landscape. Intense water repellency would have only occurred during
forest fires in the JRW and not from the slash that tree harvesting left behind
because the repellent layer tends to be destroyed at temperatures greater than
280° C (DeBano, 2000).
Agricultural development has been shown to affect several aspects of the
water cycle, which in turn change both stream flow and sediment transport. For
example, differences in tilling can affect different runoff behaviors during
precipitation events (Baker, 1983). Warnemuende, et al. (2007) showed that
fields with conventional tilling have 50% of the runoff than fields with no tilling,
however, they also show that sediment loss rates were greater for fields with
conventional tilling than those with no tilling. Wiley et al. (2010) found that if
current urban sprawl patterns, in the Muskegon Watershed (located in Michigan),
continued at their present rate that agriculture would be reduced along with some
loss of forested land cover. This would have significant effects of the hydrologic
cycle in such basins. Drainage basins with intense agriculture in Michigan

4

(greater than 60%) had less late-summer stream base flow than drainage basins
with less agriculture (Jayawickreme and Hyndman, 2007).
Models have been used to evaluate linkages between land use alterations
and changes in percolation and recharge (Selle et al., 2008; Thomas and Tellam,
2006), and to help solve water sustainability issues (Schuol et al., 2007). This
study examines the effects of land cover alterations on recharge and stream flow
within a modeling framework, and evaluates two major alterations to land cover
within the Jordan River Watershed (JRW), agriculture and forest fires. Forest fire
and agriculture represent the extreme range of expected changes that have
affected stream flow in the JRW in the past 150 years. Logging would only
remove the native canopy and not have the WR of post-forest fire landscapes
that would create large fluxes of runoff to streams (Batalla, 2001). A logged
landscape would also not have the high evapotranspiration of agriculture
landscapes (Kendall, 2009; Hyndman et al., 2007). Thus agriculture and post
forest fire landscapes would show more extreme effects on hydrology than
logging. Landslides were not included because no data could be found with
records of timing and location of historical landslides that could add sediment to
the Jordan River. The conceptual model for a forest fire that occurred in the JRW
could have created large amounts of runoff that would have also changed
recharge and stream flow. A 100% agricultural land cover in the model area is
used to test how much recharge and stream flow could have changed with the
absence of any native canopy. The native canopy was assumed to have been
completely destroyed in the conversion from forest to agricultural fields in this

5

scenario. Along with the two altered land cover scenarios, pre-developed land
cover (primarily forest) and current land cover (urban, agriculture, and forest) are
also simulated for comparison. Percent land cover for pre-developed and current
land covers were calculated using available GIS datasets. PrecipitationRecharge ratios derived from the Integrated Landscape Hydrology Model (ILHM)
(Kendall, 2009; Hyndman et al., 2007) were used to estimate how much
precipitation becomes recharge for each land cover type. The recharge, derived
from the precipitation-recharge ratios in ILHM, were then input into MODFLOW
(McDonald and Harbaugh, 1988) to simulate stream flow for each land cover
scenario discussed above. The flow simulations results are then discussed in
the context of how such modeled changes in stream flow would potentially affect
sediment transport and stream channel morphology.

Site Description
The Jordan River Watershed (JRW) (Figure 1.1) is an excellent natural
laboratory to observe the impacts of land cover change on hydrology. The JRW
has a documented history of logging, forest fires, farming, and possible
landslides since its early settlement in 1840. Residents of the basin recall a
gravel bottom river with diverse fish habitat. However, sand transported in the
river channel has reportedly filled in much of the gravel bottom, negatively
impacting fish habitat. Fish, such as Brook Trout rely on gravel for spawning,
thus this ecological attribute is highly valued (Klingeman et al., 1998).

6

Figure 1.1: The Jordan River Watershed (JRW) is located in lower Northwest
Michigan (MIGDL). For interpretation of the references to color in this and all
other figures, the reader is referred to the electronic version of this thesis.
The Jordan River drainage basin consists of outwash deposits and
course-textured end moraines. As a result, much of the precipitation within the
drainage basin infiltrates into the sediment and the vast majority of stream is
derived from groundwater rather than overland flow (Hay and Meriwether, 2004).
The peak flows in the River are caused by direct precipitation on the water or
near stream runoff after snow melt events (Figure 1.2). The Jordan is unlike
many rivers as it carries a stable flow, which varies little throughout the year.
However, the Jordan River can flood during spring thaw as ice blocks temporally
dam the river (Fuller, 2002; Sui et al., 2005). The Jordan River also has a very

7

large discharge to drainage basin area ratio, which indicates groundwater
capture from adjacent surface water drainage basins.

450

0

400
0.5

Discharge (cfs)

300
1

250
200

1.5

150
100

Precipitation (inches)

350

2

50
0
1/1/07

2.5
4/1/07

6/30/07 9/28/07 12/27/07 3/26/08 6/24/08
Date

Figure 1.2: Hydrograph, with precipitation, of the Jordan River at the United
States Geological Survey (USGS) Webster gage # 04127800. Flows and
precipitation plotted from 01/01/07-09/01/08.

The JRW has two distinct segments, the upper and lower Jordan River,
distinguished by changes in stream gradient and sinuosity. The upper part of the
Jordan River contains very steep terrain with a stream gradient of 0.76% or 7.61
m/km resulting in flow velocities of up to 1.5 m/s. The lower part of the Jordan
River has significantly lower gradients of approximately 0.080% or 0.80 m/km
and slower average measured velocities of 0.8 to 1 m/sec. The upper part of the
Jordan River has much higher sinuosity than the lower part of the Jordan River
(Figure 1.3), which is related to disequilibrium in the river elevation profile (Figure
8

1.3). Most river elevation profiles for systems that are in equilibrium can be
characterized with a smooth exponential curve. In contrast, the elevation profile
for the Jordan River has a sharp drop in gradient between 1 and 11 km from its
headwaters that could be caused by either a landslide or by the bedrock
topography. A landslide could have been initiated after removal of vegetation
during logging or after a forest fire along the steep slopes adjacent to the Jordan
River. Bedrock that is resistant to erosion can prevent the Jordan River from
rapidly down cutting, which could cause the river to migrate laterally with
significant lateral erosion of the sides of the valley. As will be discussed later, the
bedrock topography varies greatly across the model area and appears to be
influence the path of the Jordan River, especially in the upper part of the Jordan
River and in its largest tributary, Deer Creek.

9

Figure 1.3 A map of percent slope from the digital elevation model (A) for cells in
the Jordan River model area. Slopes in the Jordan River Valley (B) (shown by
the blue rectangle) contain the highest relief within the drainage basin.
Numbered locations refer points on the elevation profile for the Jordan River
(Figure 1.4) (Source, USGS SEAMLESS).

10

400

3

2

300

250

1
200

40000

30000

20000

10000

0

River Elevation (m)

350

150

Distance (m)
Figure 1.4: Elevation Profile along the Jordan River (blue line), relative to a
simple equilibrium profile (black line). The elevation profile was generated using
the points along the river that were extracted from the DEM (Figure 1.3A). The
red are locations marked on the river in Figure 1.3.
Although Deer Creek is the largest tributary to the Jordan River,
approximately half of the area of the JRW, it has relatively low flow for its surface
watershed. At the confluence of Deer Creek and the Jordan River, the flow of
the Jordan River averages approximately 180 cubic feet per second (cfs) where
Deer Creek has only 20 CFS. The base level of the Jordan River is the level of
Lake Charlevoix where it deposits its sediment.
Before settlement of the Jordan River Valley in 1840, the area was dominated
by coniferous forests (e.g. white cedar, white pine, hemlock) along the river and
deciduous forest in the upper parts of the valley (Hay and Meriwether, 2004).
The Great Chicago Fire of 1871 created a demand for timber which greatly
11

increased timber harvesting across the northern peninsula of Michigan
(Dickmann and Leefers, 2004). Up to eight logging companies operated out of
East Jordan including: Vince’s Sawmill, B&C Saw Mill, South Arm Lumber
Company, and Nelson Reddington & Co. Lumber Mill, after logging started in the
lower part of the drainage basin in 1879 (Hay and Meriweather, 2004). During the
logging period, the lower part of the Jordan River was used to carry logs into the
southern arm of Lake Charlevoix (Hay and Meriwether, 2004). By the 1900’s,
logging moved into the upper part of the drainage basin where conifers were cut
first in the lowlands and deciduous trees were cut last in the highlands. The
upper Jordan River did not have sufficient stream flow to transport these logs,
thus railroads and teams of horses were used to transport logs to Lake
Charlevoix (Hay and Meriwether, 2004). Evidence of the old logging railroads
can still be seen today (Figure 1.5). Logging continued until 1930 and took
several years to transport to Lake Charlevoix. Logs that sunk in the old harbor
can still be observed today in the vicinity where the Jordan River flows into Lake
Charlevoix.

12

Figure 1.5: Old railroad beds can be found in the Jordan River valley as linear
paths of bare ground.

The landscape of the JRW was also likely altered by forest fires (Dickmann
and Leefers, 2004). It is not known for certain that forest fires occurred within the
JRW but such events would likely have had serious implications on the
hydrology. The surface and subsurface hydrology would be changed after a
forest fire as transpiration ceases and soils become water repellent, causing an
increase of surface runoff and erosion and decreasing infiltration (Batalla, 2001).
Widespread forest fires occurred in Michigan during the 1870’s, some of
which swept across the shores of Lake Michigan to the shores of Lake Huron.
Significant forest fires were the Bad Axe Fire of 1881, the Metz Fire of 1908, and
the Au Sable-Oscoda Fires of 1911 (Dickmann and Leefers, 2004). One major
13

fire that likely occurred within the JRW is the 1911 Au Sable-Oscoda fire that was
observed in Antrim (Dickmann and Leefers, 2004). Hay and Meriwether (2004)
also reported that several fires occurred within Antrim County in the 1930’s after
logging ended, but the JRW was not mentioned.
From 1900 to 1920, data from the USDA agriculture census showed that
farming increased in the model area from 0% to more than 50% of the drainage
basin area as forest landscapes were converted to open landscapes during
timber harvesting (Figure 1.6). From 1920-1930 the agricultural area decreased
by 5% within the model area. An estimated 650 of the 1,750 original farms within
the JRW stopped operating between 1920-1930 due to nutrient poor sandy soils
(Hay and Meriwether, 2004).

60

Percent of Land Area

50
40
30
20
10
0
1900

1910

1920

1930

1940

1950

Year
Antrim

Charlevoix

Otsego

Figure 1.6: Agriculture as percent land area from 1900-1950 for the three
counties within the model area (See Figure 1.1 for county locations) (USDA).
14

Since 1930, the upper part of the JRW has returned to being mostly forest. In
1972 the Jordan River was declared the first Natural River in Michigan, which
resulted in a large amount of the land in the drainage basin being converted to
state forest (Hay and Meriwether, 2004). Currently coniferous forests are only
15% and deciduous forests are the dominant land cover type at 55% (Figure
1.7). The other major land cover types are open grass lands (15 %) and
agriculture (16%).

15

METHODS

Model Inputs
MODFLOW within the Ground Water Modeling System (GMS) was used to
construct steady-state and transient flow models of the Jordan River model area,
which was expanded from the drainage basin to include part of both arms of
Lake Charlevoix to the North, the Boyne River to the east, and the Cedar River to
the west. It was also determined that the Jordan River groundwatershed extends
well beyond its surface drainage divide to the south. A digital potentiometric
surface map was created by interpolated water levels in wells and lakes that
were connected to the groundwater using Wellogic data, from the Michigan
Geographic Data Library (MIGDL), which was then used to show direction of
groundwater flow and no-flow divides. The groundwater flow divide to the south
was represented as a no-flow boundary in the model.
Several GIS (Geographic Information System) data sets are used as
inputs for a groundwater model. Two land-cover maps, from the MIGDL, were
used to represent conditions before settlement and for the current time period.
The 1800-predevelopment map illustrates the native land covers, primarily forest
and wetlands, before large scale human disturbances (Figure 1.7). The
predevelopment map was generated from historical surveys conducted by the
MFI (Michigan Features Inventory, 2008). The current land cover used the 1978
MIRIS (Midwave Infrared Imaging Spectrograph) coverage, which was created
from aerial photographs. More recent land cover data is available with the 1992
16

IFMAP (The Integrated Forest Monitoring, Assessment, and Prescription)
coverage; however this dataset was not classified in the same manner as the
earlier data, and thus was not used.
100
90
80
Landcover (%)

70
60
50
40
30
20
10
0
Ag.

CON.

DECID.

OPEN

URBAN

WETLAND

Anderson Land Cover Type
1800 Land Cover

1976 MIRIS

1992 IFMAP

A
Figure 1.7: Land Cover Statistics (A) for the three different land cover data sets
and maps for pre-developed (B) and for current, 1978 MIRIS, (C) land covers
(MIGDL).

17

Agriculture
Coniferous Trees
Deciduous Trees
B

Open
Open Water
Rangeland
Urban
Wetlands

C
Figure 1.7 continued

18

Another major input for the groundwater model is the lithology of soils and
shallow sediments through which the water flows. Soil lithologies were used in
the calculation of precipitation-recharge ratios, which will be discussed in more
detail later, to help estimate the recharge into the glacial aquifer. The main
sediment types within the model area are sand, sandy loam, and loamy sand
(Figure 1.8). The map was generated from a digitized State soil maps by
STATSGO which was then merged with information on soil parent material
(Farrand and Bell, 1982). Most of the geologic units are glacial outwash and
course-textured tills (Figure 1.8). The highly permeable soils result in the stable
flow of the Jordan River because a large proportion of the precipitation in the
system becomes recharge, which eventually reaches the stream as groundwater
discharge.

19

Figure 1.8: GIS maps for soils (A) and surface geology (B) (MIGDL).

The Digital Elevation Model (DEM) is used to create a series of data sets
(Figure 1.4A), and as a result, is one of the most important data sets. The 1/3
arc second DEM, with a 10x10 meter cell size, was created from the National
Elevation Data set, derived from aerial photos and USGS quadrangle maps. The
20

data, which was processed by the USGS for missing data and artifacts in the
dataset, was extracted from the National Map SEAMLESS Server (non-boundary
constrained data sets).
The sources/sinks and boundary conditions for the MODFLOW model are
the Jordan River along with any of its tributaries and lakes (i.e. Lake Charlevoix
and Deer Lake). The lakes and rivers within the model were downloaded from
MIGDL (Figure 1.9). Most of the small lakes that were connected to the river
were used as draining features instead of constant head boundaries. Lake
Charlevoix was the only lake to be specified as a constant head boundary, at 177
m above MSL, in the model due to its size relative to the other lakes and its
connection with Lake Michigan. A Digital Raster Graphic of the United States
Geological Survey (USGS) 7.5 Quadrangle Map (DRG), dating 1997 (MIGDL,
2007), was used to assign the elevation of Lake Charlevoix.

21

Figure 1.9: Lake coverage of water, oil, and gas wells within the model boundary
shown in black (MIGDL). The red outline is the surface watershed.
Water wells were used to calculate starting heads along and observed
water levels were used for model optimization. Oil and gas wells were used to
calculate the bottom aquifer geometry. The water, oil, and gas well logs, which
were digitized by the Michigan Department of Environmental Quality (Figure 1.9),
provided measured water levels for a dataset of observed heads to be compared
against simulated heads from the steady-state MODFLOW model (Harbaugh et
al., 2000). Wells were filtered spatially by displaying the wells by county and any
well located outside of its county were deleted from the data set. Then the DEM
was used to filter out well logs that had differences between the DEM elevation
and the recorded elevation for ground surface greater than plus or minus one
22

meter. Well logs were also filtered by the hydrologic unit that the wells were
completed in. For example, if a well was completed in bedrock, then head
observations were not used for model comparison. A total of 190 wells were
used for water level residual analysis in the simulations.
The bottom aquifer geometry is an important boundary for the
groundwater model. A bedrock elevation map was constructed within ArcGIS
using the depth to bedrock values within oil and gas well logs, which penetrated
bedrock. The oil and gas wells were filtered based the ground surface elevation
being within 1 meter of the DEM elevation. The bedrock elevations were
estimated by subtracting the depth of bedrock from the recorded surface
elevation for each well site. Outliers, where bedrock elevations were inconsistent
with those from surrounding oil and gas wells, were removed before interpolation.
3,315 oil and gas wells were used in the creation of the bedrock elevation map
using ordinary kriging within the ArcGIS Geostatistical Analyst. Other maps such
as the USGS map created by Soller (1992) of the Glaciated United States do not
have the resolution as the maps shown in Figure 1.10 or detail on the sub-

23

kilometer scale.
Bedrock elevation (m)
176
30

A
Glacial Drift Thickness (m)
360
50

B

Figure 1.10: Drift thickness map (A) created from the interpolated bedrock
elevation (B) (MIGDL).
24

In order to evaluate stream flow response to land cover alterations,
precipitation was used as an input to estimate recharge for the groundwater
model. Hourly precipitation data was downloaded from National Climate Data
Center for the Bellaire NOAA (National Ocean Atmospheric Administration) gage
from January 1997 through June 2008. Due to the unavailability of hourly
precipitation from June through December 2008, data for this period was from
the Michigan Automated Weather Network (MAWN) at the Eastport gage,
approximately 15 km west of the model area, was used to fill in the missing data
gap. MAWN was not solely used because the gages use non-heated buckets
which lead to poor estimates of monthly precipitation during the winter. Standard
climate was used in this study, which assumes that the climate observed today
was relatively the same from the late 1800’s through the early 1900’s.
Unfortunately standard climate does not exist from the 1800’s to the early 1900’s,
thus the ten years of precipitation data from NCDC and MAWN was used as a
standard climate for the flow models with different land cover scenarios (Figure
1.11).

25

100
90

Precipatation (cm)

80
70
60
50
40
30
20
10
0

1998 1999 2000 2001 2002 2003 2004 2005 2006 2007
Year

Figure 1.11: Annual precipitation used to calculate recharge and runoff for the
Jordan River model area, includes snowfall (NCDC and MAWN).
Steady-state Model
Simulations of transient groundwater flow require starting heads and
boundary conditions. A 3-D steady-state groundwater model was constructed to
estimate starting heads for the transient groundwater model. A transient
groundwater model is able to simulate the groundwater fluxes when initial
conditions and transient data (i.e. daily precipitation) is available. A 3-D
groundwater flow model grid was constructed within GMS. The model grid
contained 724 columns by 729 rows, with approximately 50x50 meter cells. The
surface elevations from the 10x10 meter DEM were interpolated to the center of
the 50x50 m flow model cells using linear interpolation. Cells inside the

26

boundary were assigned boundary conditions (elevation, flow boundaries) and
hydraulic properties (hydraulic conductivity). Bottom elevations for each cell
were derived from the bedrock elevation map discussed earlier.
Initial conditions for MODFLOW are also needed to solve the groundwater
flow equation. MODFLOW uses the numerical solution for the partial differential
equation of Darcy’s Law for each cell (McDonald and Harbaugh, 1988). The
DEM was used to assign surface elevations as the starting heads for
MODFLOW. This is a quick method of assigning starting heads along with the
higher resolution inputs, as opposed to water wells, but it results in a longer initial
simulation time to simulate heads from this initial estimate.
A GIS coverage representing rivers was mapped over to the map module
from the GIS module in GMS. Within the map module of GMS, editing was done
to account for discrepancies where the river elevations from the DEM did not
match those from the USGS Quadrangle digital raster graphic (DRG) for any
given area. After the hydro coverage was edited to represent the Jordan River
and its tributaries, elevations were assigned from the DRG where the Jordan
River crossed elevation contours.
An average stream channel bottom conductance was assigned. A stream
width of 10 m was chosen as it was based on the average of all the cross section
data collected at stream gage sites. A value of 30 m/d, a hydraulic conductivity
value representative of sand and gravel (Freeze and Cherry, 1979; Schwartz and
Zhang, 2003), was assumed for the conductivity of the river bed sediments.
Thickness of river bed channel sediments was assumed to be 1 meter, limited to

27

the depth of the probing instrument. Equation 1 was then used to calculate the
conductance for all segments of the river (i.e. arcs of the polyline file in GIS).
The length of an arc is used by MODLFOW to assign conductance values in m2/d
for each arc.

C=

w *k *l
t

(1.1)

where;

C

is the conductance (L2/T)

w is stream width (L)
k

is the hydraulic conductivity (L/T) of stream bottom sediments

l

is the length segment (L)

t

is the thickness of the stream bottom sediments (L).

Recharge was estimated by choosing a value of stream flow from the
USGS gage that represented base flow for the Jordan River. The value of base
flow was then divided by the area of the groundwatershed for the USGS gage to
estimate a steady-state recharge value. The groundwatershed of a specific point
(i.e. the USGS gage) refers to the area from which all of the groundwater flows to
that specific point. Only one groundwatershed was used because there was only
one USGS gage available since 2007 that could be used to estimate recharge
from base flow. Using the groundwatershed for the USGS gage resulted in a
value of 0.001132 m/d for recharge, which is approximately 40 cm/yr.

28

The coverage representing the surface Quaternary geology of the area
(Figure 1.8) was mapped from a GIS to the map module of GMS. An average
hydraulic conductivity of 5 m/d was assigned as a starting value and average
value for all the geologic zones, which vary from lucustrine sand to glacial
outwash (Freeze and Cherry, 1979; Schwartz and Zhang, 2003). Parameter
estimation was used to optimize hydraulic conductivity values for the main
lithologic zones across the region.
After allowing MODFLOW to run and converge on a solution, optimization
was used to estimate the most important parameters within the groundwater
model. Residential water well levels (Figure 1.9) for locations within the model
area were plotted against simulated heads in those model cells. Preliminary
ground water simulations were run separately for specific yield and hydraulic
conductivity to test how of each of these hydraulic properties changed the
simulated heads with respect to the observed heads. Simulations that only
changed specific yield had little effect on the improving the differences between
the observed and simulated head. However, simulations that changed only
hydraulic conductivity improved the differences between the simulated and
observed head thus only hydraulic conductivity was used in the optimization of
the simulated heads. Four PEST (Parameter Estimation) simulations were run,
within GMS, beginning with a hydraulic conductivity for all the geologic
lithologies. Additional values of hydraulic conductivity were then added, starting
with 2 and then 3, until the sum-square residual error between simulated and
observed heads was no longer reduced (Table 1). Figure 1.12 shows a spatial

29

plot of heads and a graph of simulated versus observed heads for the best PEST
simulation.
Table 1.1: Values for hydraulic conductivity after optimization for geologic zones
shown in Figure 8B. The starting hydraulic conductivity value for all the geologic
zones for the optimization was 5 m/d.
Geologic Zones
Glacial outwash sand and gravel and postglacial alluvium
Lucustrine sand and gravel
Coarse-textured glacial till
End moraines of coarse-textured till

Hydraulic Conductivity (m/d)
6.34
1.12
0.55
2.81

A
Figure 1.12: Simulated map of heads (A), along with a comparison between
simulated vs. observed heads after PEST optimization (B).

30

400

Simulated Head (m)

350

300

250

200

y = 0.973x + 6.95
R² = 0.984

150
150

200

250
300
Observed head (m)

350

400

B
Figure 1.12 continued

The groundwatershed (Figure1.13) was drawn from the spatial plot of
heads from the steady-state model (Figure 1.12) to represent the capture zone of
the Jordan River system down to the headwaters of Lake Charlevoix. As
discussed in the in the methods section, the southern boundary was drawn along
a flow divide which was represented as a no-flow boundary.

31

Model Boundary
Surface Watershed
Boundary
Ground Watershed

Figure 1.13: Shows the JRW and groundwatershed (shaded) for the Jordan River
within the model area.

Transient Field Data
Transient flow data collected from summer 2007 through fall 2008 were
used to characterize the transient groundwater discharge and stream flow for the
JRW. Flows were collected across the JRW to measure changes in channel
geometry and flow for the Jordan River, and Deer Creek.
Gages installed at three locations, two on the Jordan River and one on
Deer Creek, provide nearly continuous estimates of stream flow to characterize
the detailed transient hydrology in the drainage basin. The discharge from these
gages combined with the USGS gage was used to analyze the stream flow
response to precipitation and snow melt events in different reaches of the JRW.
The gages at each site are pressure transducers that convert pressure created
32

from the water column above the transducer into water column height which
represents the stage of the river. Stream discharge was measured at multiple
stages and seasons, starting in 2007 and ending in 2008, for each site to create
a stage-discharge relationship (Figure 1.14). An important assumption for these
relationships is that the stream cross-channel geometry does not change over
the period with which the data is collected. Stage-discharge relationships are
commonly logarithmic because a river ultimately overflows its banks, reducing
the effect of increasing stage on discharge. Ideally, a wide range of flow values
would be available to establish this type of relationship, but the gage at Graves
and Fisheries gages only have relatively low flows observed. Deer Creek has a
well established stage-discharge relationship with stream flow observed from low
to high values, being the gage with the best accuracy. High flow events have not
yet been captured for the Jordan River because of the short lag between
precipitation and peak flows, on the order of only a few hours.
As a part of the flow measurements, cross channel geometry was
measured at multiple sites, to characterize how stream flow and bed geometry
changed across the watershed (Figure 1.15). Flow and bed geometry was also
measured at the gage sites multiple times to see how bed channel geometry
changed with time (Figure 1.16). These data indicate likely scouring between
flow measurements at two of the gage sites, Deer Creek and Graves (Figure
1.16).

33

Deer Creek

USGS 04127800

Fisheries

Graves
A

Figure 1.14: Map showing gage locations within the Jordan River Watershed
(highlighted in red) (A) along with their stage-discharge relationships: USGS
Webster (B) Graves (C), Fisheries (D), and Deer (E).

34

700

Discharge (cfs)

600
500
400
300

y = 27.831e

1.904x

2

R = 0.893

200
100
0.8

1

1.2

1.4

1.6

1.8

Stage (m)

B
170

Discharge (cfs)

160
150
140

y = 301.24x + 106.18

130

2

R = 0.336

120
110
0.06

0.08

0.1

0.12
Stage (m)

C
Figure 1.14 continued
35

0.14

0.16

0.18

75
70

Discharge (cfs)

65
60
55

y = 402.88x - 34.35
2

R = 0.905

50
45
40
0.19

0.2

0.21

0.22

0.23

0.24

0.25

0.26

Stage (m)

D
180
160
140

Discharge (cfs)

120
100
80
y = 5.407e

60

5.606x

2

R = 0.992

40
20
0
0.2

0.3

0.4

0.5
Stage (m)

E
Figure 1.14 continued
36

0.6

0.7

Rogers Rd.

Wesbter Br.

Old State Rd.

Graves Crossing

Headwaters

Penny Br.

Position (m)
0

3

6

9

12

15

18

0

Depth (m)

0.2
0.4
0.6
Headwaters
Penny Br.
Graves Crossing
Old State Rd.
Webster Br.
Rogers Rd.

0.8
1
1.2

Figure 1.15: Spatial channel cross channel geometries for all the sites where
flows were conducted.

37

South
0
1
0

2

Position (m)
3
4
5
6

7

8

North
9 10

Depth (m)

0.1
0.2
0.3
0.4
0.5
A

0.6
07/01/07

09/01/07

09/17/07

04/11/08

07/06/08

09/29/08

Position (m)

East
00

01/12/08

2

4

6

8

10

12

14

West
16

0.1
0.2

Depth (m)

0.3
0.4
0.5
0.6
0.7
0.8
B 0.9

Figure 1.16: Temporal cross channel geometry at Graves (A) and Deer Creek (B)
gage locations. See Figure 1.14A for gage locations.

38

Through the course of this study flow from all the gages showed rapid
response to precipitation events for Jordan River, while Deer Creek and the
USGS Webster site can have relatively long lags of a few days (Figure 1.17). It
should be noted that spike in precipitation around 8/6/07 is related to a rain gage
malfunction (indicated by large precipitation events) on that date (NOAA, 2008).
400

0.4

300

0.8
200
1.2
100

1.6
2

Flow (cfs)

Precipitation (inches)

0

8/6/07 8/6/07 8/6/07 8/6/07 8/6/07 8/6/07 8/6/07
Date

0

Precipitation
NOAA Precipitation
Gauge Flow

Measured Flow
USGS Webster
Graves
Fisheries
Deer

USGS Webster
Graves
Fisheries
Deer

Figure 1.17: Discharge at all four stream gages in response to Precipitation
Events. Arrows indicate lags (MAWN, NCDC, and USGS).

39

Recharge Calculations for Land Cover
It was important to characterize how recharge to the aquifer changed as a
result of land cover. Several studies, previously mentioned, have been able to
show how land cover affects groundwater recharge indirectly. Kendall (2009)
used ILHM to simulate the direct effects of land cover on groundwater and
stream flow in the Muskegon River groundwatershed. However, the
development of an ILHM model, including gathering the necessary data, is very
time intensive. Due to time constraints, a simplified approach that used
precipitation-recharge ratios from ILHM was used to calculate recharge from the
Jordan River precipitation data. Precipitation-recharge ratios allow an indirect
estimate of how much precipitation in a given period becomes recharge. The
calculation of the precipitation-recharge ratios accounts for sediment type
(infiltration), and land cover (interception and evapotranspiration). The
precipitation data set for the Muskegon drainage basin consisted of monthly
precipitation for 27 years. ILHM then used the precipitation data for the
Muskegon drainage basin to simulate monthly recharge for 27 years. In order to
apply the precipitation-recharge ratios to the Jordan River precipitation data, a
few steps were needed to estimate the ratios from the precipitation data and
recharge simulated from ILHM. First, the effects of the snow pack accumulation
and the large amount of melt that occurs in the Jordan River had to be examined.
A simple way to simulate the effect of snow-pack accumulation in excel was to
assign a constant-averaged value of precipitation to December, January,
February, and March of every year. The constant-averaged value of precipitation

40

was calculated by summing all the precipitation for those 4 months in that year
and then dividing by 4. The same value was then applied to all four months. The
next step was to take the recharge simulated for the Muskegon Drainage basin
by ILHM for a given month and divide by the precipitation for the Muskegon
drainage basin for that same month. After all of the precipitation-recharge ratios
were calculated for the 27 year data period, each value was averaged together to
calculate monthly precipitation-recharge ratios that were independent of
precipitation data. This process was repeated for each land cover type in the
Anderson Classification System: open, agriculture, coniferous forest, deciduous
forest, urban, and wetland.
The precipitation-recharge ratios were then used to calculate recharge
from the precipitation data for the model area. Monthly precipitation was used
instead of hourly, daily, or weekly because of the lack of information and data
that would be required to simulate of those time scales (i.e. evapotranspiration
estimates). Monthly precipitation was multiplied by the precipitation-recharge
ratio to calculate recharge for each month from January 1, 2007 through
September 1, 2008.
Due to the expected water repellent nature in soils after forest fires, an
approach was developed to evaluate the effect of repellent soils that might have
existed after forest fires, ultimately affecting the recharge. A runoff relationship
which assumed high severity of burned soils (i.e. maximum WR), modified from
Onda (2008), was used to evaluate runoff and infiltration for the simulation of
post forest fire conditions in the model area. Onda studied the effects of runoff

41

after large storm events the season after a severe forest fire in terrain similar to
the model area (steep slopes with sandy sediment), and devised devices that
could measure runoff during rainfall events. The amount of runoff was then
compared to the amount precipitation that a single rainfall event generated
(Table 1.2).

Table 1.2: Total rain and runoff data collected by Onda on hill slopes near Point
Reyes, California immediately following a forest fire. This study used rainfall and
runoff data to calculate runoff ratios which was the basis for the precipitationrunoff relationship. Values in red indicate the collection tank was full and did not
actually record all the rain. Event 1 recorded 6.51 inches, but the actual rainfall
was 7.96 inches. Event 6 recorded 1.43 inches but the total rainfall was 2.46
inches. Thus, 7.96 and 2.46 inches were used in the runoff-precipitation
relationship. Modified from Onda (2008).

Rainfall Event
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17

Total Rain (inches)
7.96
0.72
0.57
0.70
0.76
2.46
0.20
0.42
0.34
0.28
0.13
0.38
0.33
0.22
0.24
0.32
0.28

42

Runoff (inches)
0.74
0.44
0.17
0.24
0.21
0.74
0.04
0.15
0.27
0.09
0.02
0.06
0.04
0.07
0.09
0.02
0.04

Run off Ratio
0.19
0.62
0.3
0.35
0.27
0.52
0.18
0.36
0.81
0.32
0.17
0.16
0.12
0.31
0.37
0.07
0.15

The relationship from Onda (2008) was modified by removing rainfall
events larger than 2.5 inches, as Jordan River receives considerably less in a
rainfall event than the Point Reyes Peninsula in California where the study took
place. The total rainfall versus runoff in Table 1.2 was then plotted relative to
each other (Figure 1.18) to estimate a relationship between precipitation and
runoff. The relationship was used to first filter out runoff from the precipitation
data used in this study to simulate runoff in a post-forest fire land cover. The
difference between daily rainfall and daily runoff was termed daily infiltration.
The infiltration was then summed up for every month in the time-interval for this
study (January 2007-September 2008). To estimate the amount of recharge
from monthly infiltration, precipitation-recharge ratios from open land cover,
based on ILHM, were used to simulate recharge for the post-forest fire land cover
scenario.

43

0.8
0.7

Runoff (inches)

0.6
0.5
0.4
y = 0.244Ln(x) + 0.396
R2 = 0.771

0.3
0.2
0.1
0
0

0.5

1
1.5
Precipitation (inches)

2

2.5

Figure 1.18: Precipitation: Runoff relationship used to calculate infiltration and
runoff. Modified from Onda (2008).
Transient Model
This model describes how base flow from groundwater changes as the
land cover changes. Any alterations in the land cover change the percent of
precipitation that becomes recharge.
The transient groundwater model was developed by taking the head
solution from the steady-state MODFLOW simulation and using it as the starting
heads of the transient MODFLOW simulation. A general value of 0.2 was used
as the specific yield for the near surface aquifer, as no data was available that
would enable estimates for specific yield. Sensitivities were also analyzed by
scaling different parameters. Two sensitivities scenarios were examined: 1)
Decreasing the hydraulic conductivities by 50%, and 2) changing the specific
44

yield (S y ) from 0.2 to 0.35. Minimal sensitivity was found when increasing S y but
all the simulated stream flows decreased significantly when hydraulic
conductivities were decreased.
For each land cover scenario, monthly recharge values were used as
inputs into MODFLOW. The transient model was set up to run with monthly
recharge, calculated from the precipitation data, from 01/01/07 through 09/01/08
so data from the four gages within the Jordan River could be used as
observations. As a result, 21 stress periods were used with only 1 time step so
that only monthly groundwater fluxes would be written out.

A MATLAB script

was then be used to read in the groundwater fluxes from MODFLOW to sum all
the stream flow at all the gage locations within the JRW.
5
4.5
Recharge (inches/month)

4
3.5
3
2.5
2
1.5
1
0.5
0
01/01/07

04/11/07
Agriculture

07/20/07

10/28/07
Date
Pre-developed

02/05/08
Current

05/15/08

08/23/08

Post Forest Fire

Figure 1.19: Approximate recharge for estimated for different land cover
scenarios within the model area.
45

Recharge for the land cover scenarios of interest: pre-developed (1800),
post forest Fire (PFF), agriculture (AG), and current, where used as inputs into
the transient MODFLOW model. Simulated base flow, calculated using the
scripts and groundwater fluxes into the stream, was then matched with gage
data. Because the transient model averaged recharge over the majority of the
model area, the gage with the largest catchment, the USGS gage, was used for
the observed stream flow.
The simulated flow for the current land cover for the model area was
plotted with the gage flow at the USGS Webster gage and base flow. The base
flow was calculated using the local minimum base-flow separation method within
the Web-based Hydrograph Analysis Tool (WHAT) (Lim et al., 2005) with a
statistical comparison of the Eckhardt and BFLOW filter for 50 USGS gaging
stations in Indiana. Lim et al., 2005, found a Nash-Sutcliffe coefficient over 0.91
for both filters applied to the data for all 50 gaging locations, concluding that
WHAT will typically match most other traditional techniques that estimate base
flow from stream flow. Matched base flow can lead to error, based on the
method, where WHAT provides consistent results within minutes and is linked
into the USGS gage server (Lim et al., 2008).

46

RESULTS AND DISCUSSION
Model Results
The simulated base flow for the current land cover within the model area
provide a general representation of the base flow calculated from the USGS
Webster gage using WHAT (Figure 1.20). Stream evaporation was not taken
into account in the hydrologic simulations, which is part of the reason the
simulated stream flow is higher than the base flow separated or gage stream
flow. Other explanations for the simulated flow being higher than the observed
stream flow are using a recharge for the model area that is greater than the
catchment area for the USGS gage. Using MODFLOW, stream flow at the
USGS Webster gage was simulated for the other land cover scenarios (Figure
1.21).

47

400

Discharge (cfs)

350
300
250
200
150
100
01/01/07

04/11/07

07/20/07

10/28/07

02/05/08

05/15/08

08/23/08

Date
USGS Gage Flow

Base Seperated Flow

Simulated Flow

Figure 1.20: Simulated Flow for Current Land cover within the Jordan River
model area plotted with stream flow and base flow for the USGS Webster gage
(USGS).
230
220

Discharge (CFS)

210
200
190
180
170
160
150
140
1/1/07

4/11/07

7/20/07

10/28/07

2/5/08

5/15/08

8/23/08

Date
Current Landcover

Post Forest Fire

Agriculture

Pre-developed

Figure 1.21: Simulated stream flow at USGS Webster for all four tested
landscapes.
48

The post forest fire landscape showed the most significant deviation in
simulated stream flow from current land cover. The large difference in stream
flow between the current and post forest fire landscapes is the result of the
calculated runoff being added to the drainage basin from the method described in
Onda (2008). This study appears to be the only documented use of this
relationship, besides Onda, and this is only a gross estimate as such
relationships are not readily transferable to other watersheds. Simulated stream
flow for all the other landscapes (Ag, 1800) did not vary significantly from the
current land cover scenario.

Implications for Sediment Transport
The ability to simulate stream flow was important to provide insights on
how stream flow could have changed bed load transport. The stream flow affects
the size of transported particles as well as the volume of sediment transported.
Thus, a higher mean stream flow would have more capacity to transport
sediment downstream (Figure 1.22). This is called potential and is controlled by
the stream power, which is calculated using the following equation (Bagnold,
1966):

ω= ρgdSu
where;

ρ is fluid density,
g is the acceleration of gravity,

49

(1.2)

d is flow depth,
S is the energy gradient of the flow (stream gradient), and
u is the mean velocity of the flow.

Figure 1.22: Hjulstrom’s diagram showing velocities required for the entrainment
and deposition of different sediment sizes. Modified from Prothero and Schwab
(2004).
The Jordan River clearly has the potential to transport sediment, thus
there must be a source of sand that is supplying the sediment at high discharge
events. Two of the most likely scenarios, for the JRW are landslides (due to
steep slopes and deforestation) and runoff occurring soon after a forest fire.
After examining the stream flow from different land cover scenarios, forest
fires appear to be the most likely scenario to have large changes in stream flow
and sediment transport within the JRW. Massive amounts of runoff and high
50

stream flow from large precipitation events after a forest fire could have delivered
large loads of sediment to the river and throughout the drainage basin. The
sediment left from the unusually high runoff and stream events could still be
entrained at high stream flow and then transported at base flow. Unfortunately,
runoff from the hydrophobic soils from the post-forest fire landscape was not
directly evaluated as data within the region did not exist.

51

CONCLUSIONS

A simple modeling approach using MODFLOW was used to simulate
stream flow for the different land cover scenarios: pre-development, post-forest
fire, agriculture, and current. These four land cover scenarios were chosen,
based on the history of the area, to have the most potential impact on stream
flow and sediment transport. The simulations over estimated stream flow
relative to the USGS gage flow and base flow using base-flow separation. The
resulting groundwater shed for the Jordan River was found to be about 30%
larger than the drainage basin, making changes in land cover outside the
drainage basin just as important as those within the basin. Simulations of
changes within the past 150 years, such as deforestation, agriculture, and
urbanization did not show significant change in stream flow. Although changing
recharge via land cover alterations does not seem to have a large effect on the
sediment transport of the Jordan River, runoff in the post-forest fire land cover
did show the largest change in stream flow and thus a large change in sediment
transport. The runoff generated from high precipitation events after a forest fire,
along with impervious surfaces that have been added since forest fires (i.e.
roads), could be significant enough to have added large amounts of
sedimentation to the river which are still being transported today.
Future work will include developing a detailed ILHM model to directly
analyze land cover effects on stream flow and sediment transport. ILHM
simulates the terrestrial hydrologic cycle using hourly precipitation and provides
52

estimates of hourly recharge, runoff, and stream flow for different land cover
scenarios. Codes such as Hydrologic Engineering Centers River Analysis
System (HEC-RAS) can then model sediment transport and sediment source
potential using the hourly stream flow from ILHM. The ability to model hourly
runoff and stream flow will help answer the question of how much sediment
would get eroded and transported in high precipitation events after a forest fire in
the Jordan River groundwatershed. HEC-RAS will also be used to identify
potential sources of sediment loading (i.e. landslides).

53

BIBLIOGRAPHY

54

The Census of Agriculture. United States Department of Agriculture. 2010. July
2008 < http://www.agcensus.usda.gov/>
Bagnold, R. “An approach to the sediment transport problem from general
physics.”
United States Geological Survey Professional Paper 4221 (1966).
Baker, J.L., “Agriculture areas as nonpoint sources of pollution in Environmental
Impact of Nonpoint Source Pollution.” Ann Arbor, MI: Ann Arbor Science
Publishers: 275-310.
Batalla, R. “Hydrological Implications of Forest Fires.” Proceedings from 3rd
International Summer School on the Environment. (2001).
Bosch, J, Hewlett, J. “A review of catchment experiments to determine the effect
of vegetation changes on water yield and evapotranspiration.” Journal of
Hydrology 55 (1982) 3-23.
Dickmann, D., Leefers, L., The Forests of Michigan, The University of
Michigan Press, 2004.
S., Dube, A., Plamondan, R., Rothwell. “Watering Up After Clear-cutting on
Forested Wetlands of the St-Lawrence Lowland.” Water Resources Research 31
(1995) 1741-1750.
Cassie, D., Jolicoeur, S., Bouchard, M., Poncet E. “Comparison of Stream flow
between pre and post timber harvesting in Catamaran Brook (Canada).” Journal
of Hydrology 258 (2002) 232-248.
Farrand, W., and Bell, D., Quaternary geology of Michigan. Michigan Department
of
Natural Resources, Geological Survey Division map, Lansing. 1982.
Freeze, R., Cherry, J. “Groundwater.” Prentice Hall, Inc. Upper Saddle River, NJ.
(1979).
Fuller, D., Jordan Valley Voices, Resource Book, The River, 2002.
Harbaugh, A.W., Banta, E.R., Hill, M.C., and McDonald, M.G., MODFLOW-2000,
The U.S. Geological Survey Modular Ground-Water Model-User Guide to
Modularization Concepts and The Ground-Water Flow Process (2000).
Hay, R. L., and M. Meriwether. 2004. Jordan River Assessment. Michigan
Department of Natural Resources, Fisheries Special Report 28, Ann Arbor.

55

Hibbert, A. 1967. Forest treatment effects on water yield in Forest hydrology in
proceedings of a National Science Foundation advanced science seminar.” New
York: Pergamon Press: 527-543.
Hyndman, D., Kendall, A., and Welty N. “Evaluating Temporal and Spatial
Variations in Recharge and Stream flow Using the Integrated Landscape
Hydrology
Model (ILHM).” Subsurface Hydrology: Data Integration for Properties and
Processes,
AGU Geophysical Monograph Series 171 (2007) 121-141.
Jayawickreme, D., Hyndman D. “Evaluating the Influence of Land Cover on
Seasonal Water Budgets using NEXRAD and Stream flow Data.” Water
Resources
Research 43 (2007).
Johnson, A. “Ground-water Response to Forest Harvest: Implications for
Hillslope Stability.” Journal of the American Water Resources Association. 43
(2007) 134-147.
Kasai, M. “Channel processes following land use changes in a degrading steep
headwater stream in North Island, New Zealand.” Geomorphology 81 (2006)
421–439.
Klingeman, Peter C., Beschta, Robert L., Komar, Paul D., Bradley, and Jeffery
B.,
Gravel-Bed Rivers in the Environment, Water Resources Publications, LLC.,
1998.
Lim, K., Engel, B., Tang, Z., Choi, J., Kim, K., Muthukrishnan, S., Tripathy, D.
“Automated Web GIS Based Hydrograph Analysis Tool, WHAT.” Journal of the
American Water Resources Association (2005) 1407-1416.
Michigan Automated Weather Network. Michigan State University. 2010.
September 2008 < http://www.agweather.geo.msu.edu/mawn/>
Michigan Geogrpahic Data Library. State of Michigan. 2010. 1
December 2008 < http://www.mcgi.state.mi.us/MIGDL/>
Mikos, M., Fazarinc, R., Ribicic, M. “Sediment production and delivery from
recent large landslides and earthquake-induced rock falls in the Upper Soca
River Valley, Solvenia. Engineering Geology 86 (2006) 198–210.
Munro. 1869. Munro Family notes. Unpublished hand written transcript, East
Jordan
Library.
56

Mutie, S., Mati, B. Home, P. Gadain, H. and Gathenya, J. “Evaluating Land Use
Change Effects on River Flow Using USGS Geospatial Stream Flow Model in
Mara
River Basin, Kenya.” Exploration Geophysics 22 (1991): 19-22.
National Climatic Data Center. National Oceanic and Atmospheric
Administration.. 2010. September 2008 <
http://www.ncdc.noaa.gov/oa/ncdc.html>
National Map SEAMLESS Server. United States Geological Survey. 2010. July
2007.
< http://seamless.usgs.gov/>
National Water Information System. United States Geological Survey. 2007. 1
December 2007 < http://waterdata.usgs.gov/nwis/uv?04127800>
Pierson, F., Robichaud, P., Moffet, C., Spaeth, K., Williams, C., Hardegree, S.,
Clark, P. “Soil Water Repellency and Infiltration in Course-textured Soils of
Burned and Unburned Sagebrush Ecosystems.” Catena 74 (2008) 98-108.
Poff, N., Bledsoe, B., Cuhaciyan, C. “Hydrologic variation with land use across
the contiguous United States: Geomorphic and ecological consequences for
stream ecosystems.” Geomorphology 79 (2006) 264–285.
Prothero, D., Schwab, F. Sedimentary Geology, W.H. Freeman and Company,
2004.
Saito, H., Nakayama, D., Matsuyama, H. “Relationship between the initiation of a
shallow landslide and rainfall intensity-duration thresholds in Japan.”
Geoinformation 12 (2009) 487-495
Schuol, J., Abbaspour, K. “Using Monthly Weather Statistics to Generate Daily
Data in a SWAT Model Application to West Africa.” Ecological Modeling 201
(2007) 301-311.
Schwartz, F.W., Zhang, H., “Fundamentals of Groundwater.” John Wiley and
Sons, Inc. New York, NY. (2003).
Selle, B., Lischeid, G., Huwe, B. “Effective Modeling of Percolation at the
Landscape Scale using Data-based Approaches.” Computers and Geosciences
34 (2008) 699-713.
Sophocleous, M. “Groundwater recharge and sustainability in the High Plains
aquifer in Kansas, USA.” Hydrology Journal 12 (2005) 351-365.

57

Sui, J., Karney, B. Fang, D. “Variation in water level under ice-jammed conditionfield
investigation and experimental study.” Nordic Hydrology 36 (2005) 65-84.
Thomas, A., Tellam, J. “Modeling of Recharge and Pollutant Fluxes to Urban
Groundwaters.” Science of the Total Environment 360 (2006) 158-179.
USGS Real-Time Water Data for the Nation. United States Geological Society.
2010. December 2008. < http://waterdata.usgs.gov/nwis/rt>

Waterloo, M., Schellekens, J., Bruijnzeel, L., Rawaqa, T. “Changes in
catchment runoff after harvesting and burning of a Pinus caribaea plantation in
Viti Levu, Fiji.” Forest Ecology and Management 251 (2007) 31–44.
Warnemuende, E., Patterson, J., Smith, D., Huang, C. “Effects of tilling no-till soil
on losses of atrazine and glyphosate to runoff water under variable intensity
simulated rainfall.” Soil & Tillage Research 95 (2007) 19-26.
Wiley, M., Hyndman, D. Pijanowski, B.,Kendall, A., Riseng, C., Rutherford, E.,
Cheng, S., Carlson, M., Tyler, J., Stevenson, J., Steen, P., Richards, P.,
Steelbach, P., Koches, J., Rediske, R. “A multi-modeling approach to evaluating
climate and land use change impacts in a Great Lakes River Basin.”
Hydrobiologia 657 (2010) 243-262.
Winkler, R., Spittlehouse, D., Golding, D. “Measured differences in snow
accumulation and melt among clearcut, juvenile, and mature forests in southern
British Columbia.” Hydrological Processes 19 (2005) 51-62.

58

CHAPTER 2: USING RADAR STRATIGRAPHY TO EXAMINE THE
RESPONSE OF HISTORICAL SEDIMENT TRANSPORT TO LAND
COVER ALTERATIONS

59

INTRODUCTION

Landscape alterations through human disturbance can have a major effect
on sediment accumulation in streams and rivers (Liebault et al., 2005; Noel et al.,
2001; Lach and Wyzga, 2002; Kasai et al., 2005). Decisions about managing
channels and slope erosion can be greatly improved if the link between sediment
transport response and changes in both climate (i.e. precipitation) and human
disturbance (i.e. land use) are clearly understood. Kasai et al. (2005) showed
that major changes in sediment accumulation via deforestation altered channel
morphology after aforestation occurred.
In order to understand how climatic and anthropogenic forcing will affect
sediment transport, the effects of climate and land cover on past sediment
transport was evaluated. Dearing and Jones (2003) demonstrated the
importance of knowing the sediment transport responses to climate and
anthropogenic forcing. Deforestation and agricultural expansion can be a
primary cause of changes in fluvial sediment transport rates, while climatic
events play only a secondary role (Diodato, 2006). Foster et al. (2003) showed
that over long time scales, landscape alterations by human disturbance can
cause a larger shift in sediment movement within a stream, relative to climatic
forcing. Syvitski and Milliman (2007) were able to correlate modeled changes in
sediment yield due to landscape alterations, but their effort was limited by the
need for higher resolution record of historic sediment transport behavior.

60

Deltas can provide detailed records of fluvial sediment transport dynamics,
as sediments are eroded from the drainage basin and transported to the delta
(Prothero and Schwab, 2004). If a specific volume of sediment within a delta can
be dated, sediment transport rates and timing of deposition can be estimated.
High resolution geophysical techniques, such as ground penetrating radar
(GPR), can be a valuable tool for mapping deltaic structures. GPR can image
variations in subsurface structure in the meter to centimeter scale (Knight, 2001),
which is necessary to characterize deltaic structures. GPR has been used in a
number of deltas (e.g., Jol and Smith, 1991) to map out 2-D and 3-D depositional
structures that may then be correlated to changes in base-level (Bersezio et al.,
2007) and fluvial geomorphology (Beres et al., 1999). Pelpoa and Hicken (2004)
used GPR combined with sequences of historical aerial photographs to estimate
sediment volume and mass flux rates over a 52 year period.
The objective of this study is to characterize subsurface features related to
fluvial or deltaic sedimentary environments using GPR. Another goal of this study
is to provide information, such as locations of sediment types, and depths of
different sediment layers to help estimate deposition chronology and
sedimentation rates. The information from this study can also help estimate the
total volume of sediment in the delta and historic sediment transport rates in the
Jordan River.

61

Site Description
The Jordan River (Figure 2.1) has likely undergone changes in position
and sediment transport over the past 150 years due to changes in land cover,
climate, and dredging of the channel between Lake Michigan and Lake
Charlevoix in 1876 (Lake Charlevoix Association, 2008).

Figure 2.1: The Jordan River Drainage basin is in the northern portion of
Michigan’s Lower Peninsula (MIGDL).
The dredging of the channel between Lake Michigan and Lake Charlevoix
caused a drop in baselevel, which would have caused a change in the position of
sediment deposition where the Jordan River flows into the lake. Before the
channel was dredged in 1876, Lake Charlevoix level was approximately at 177.5
meters above sea level (Lake Charlevoix Association, 2008; see Figure 2.2).
After the channel was fully dredged, the Lake Charlevoix level was equal to the
level of Lake Michigan, and has since fluctuated due to changes in climate and
62

human activities. Figure 2.2B shows how the delta appears during periods of
high and low lake levels. Lake level affects sediment deposition by changing the
elevation where a velocity drop occurs, which causes it to drop much of its
sediment load. When the water level for Lake Charlevoix was higher, prior to the
dredging of the channel, the Jordan River would have reached near zero velocity
at a point upstream from the current zone of deposition.

63

Figure 2.2: Lake Level data for the NOAA gage for Lake Michigan and Lake
Huron. The dotted line shows the approximate level of Lake Charlevoix before
the opening of the channel between Lake Michigan and Lake Charlevoix
(NOAA). Lake levels are monthly averages.

64

Other disturbances that could have caused changes in the evolution of the
delta are alterations to land cover during the late 1800’s and early 1900’s
resulting from logging, forest fires, possible landslides, and agricultural
development.
Before geophysical measurements were conducted, the study area was
surveyed with a Total Station at the origin of the grid (500, 500). Points that were
not surveyed with the total station were surveyed manually by extending a
measuring tape between the stakes that were previously surveyed with the total
station. Then, additional stakes were placed every 20 meters until the northern
section of the study area was surveyed. 26 lines were laid out along the northern
section of the delta as shown in Figure 2.3. The section of the delta was chosen
based on preliminary GPR data collected in the winter of 2007, where good
reflection data was collected in the northern part of the survey area.

65

Figure 2.3: Study area showing GPR lines, auger boreholes, and Vibracoring
locations adjacent to the Jordan River and Lake Charlevoix. The numbered
locations will be referred to in the text (MIGDL).

66

METHODS

Ground Penetrating Radar
Ground Penetrating Radar (GPR) is a useful method for studying the
geometry and relationships between certain stratigraphic facies (Jol and Smith,
1991). GPR works by sending electromagnetic (EM) waves from a transmitting
antenna into the subsurface. These waves then reflect off interfaces between
layers with contrasting dielectric properties and back to the receiving antenna,
with some energy refracting into deeper layers. Electrical properties of soils and
rocks are controlled by factors such as porosity, clay content, dissolved minerals,
and water content (Beres and Haeni, 2000).
The reflection of a wave can be described using the reflection
coefficient, as mentioned by Davis and Annan (1988):

R=

Z 2 − Z1
Z 2 + Z1
(2.1)

Z =(

jωμ 1/ 2
)
σ + jωε

ε = ε 0ε r

(2.2)

(2.3)

where;
Z is impedance

ε 0 is the permittivity of a vacuum (8.85419x10-12 F/m)
67

ε r is the relative permittivity.
σ

is electrical conductivity (specified as a constant)

Another control on the penetration of GPR reflections is the magnetic
permeability of a material:

μ = μ0 μ r

(2.4)

where;

μ0 is the magnetic permeability of a vacuum (4 π x10-7 H/m)

μr is the relative magnetic permeability.
EM wave velocities can be used to estimate the depths for the GPR
reflections and can help identify the material (Table 2.1) through which the waves
propagated (Baker 1999). EM wave velocities can be calculated using:

Vr =

c

(2.5)

1/ 2

(ε )

where;

Vr

is the EM propagation velocity of the media,

ε is the dielectric constant, and
c is the velocity of EM waves in a vacuum (3x108 m/s).

68

Table 2.1: Various materials and their relative permittivity and velocity, Modified
from Baker (1991).
Materials
Relative Permittivity Velocity (ft/ns) Velocity (m/ns)
Fresh Water
80
0.11
0.03
Soil, sandy wet
15-30
0.18-0.25
0.05-0.08
Sand, wet
20-30
0.18-0.31
0.05-0.09
Silts
5-30
0.18-0.44
0.05-0.13
Clays
5-40
0.16-0.44
0.05-0.13
Clay, wet
15-40
0.16-0.25
0.05-0.13
Soil, loamy wet
10-20
0.22-0.31
0.07-0.09
Soil, clayey wet
10-15
0.25-0.31
0.08-0.09
Shales
5-15
0.25-0.44
0.08-0.13
Snow
80
0.28-0.35
0.09-0.11
Soil, sandy dry
4-6
0.40-0.49
0.12-0.15
Soil, loamy dry
4-6
0.40-0.49
0.12-0.15
Soil, clayey dry
4-6
0.40-0.49
0.12-0.15
Sand, dry
3-5
0.40-0.57
0.12-0.17
Clay, dry
2-6
0.40-0.70
0.12-0.21
Fresh Water Ice
3-4
0.49-0.57
0.15-0.17
Air
1
0.98
0.30

GPR Data Collection
The GPR data was collected in March 2008, with 163 files over 38 lines.
This GPR data was collected in the winter as a hard snow pack covers up
uneven terrain, leaving a flat surface to collect GPR, and we could collect data
over the water and other areas that are normally only accessible by boat.
The data was collected by pulling a pair of 200MHz bistatic antennas
mounted on a sled behind a snowmobile (Figure 2.4). The 200 MHz frequency
was based on preliminary GPR data collected in the winter of 2007 where 100
and 200 MHz antennas were tested The 200 MHz data showed more contrast in
the subsurface reflections than the 100 MHz data.

69

The system was triggered every 10 cm by a wheel that followed the sled.
Due to the low friction surface of the snow pack, wheel slip occurred. This was
accounted for by marking traces during data collection when the center of the
antenna array is at marked stakes in the field. The locations were then corrected
using a processing step called “rubbersheeting”. Some lines through thick
vegetation were collected manually by holding the bistatic-antenna at a fixed
distance of 50 cm. Data was collected using a 200 ns time window with 32 stacks
and a 0.01 meter step size.

Figure 2.4: GPR antenna setup during data acquisition in the winter of 2008.
Photo by Christopher L. May.

70

GPR Processing
After GPR data were collected, processing was done to remove artifacts in
the data and to calculate depths so the data can be analyzed and interpreted.
REFLEX was used for all the processing as described in Figure 2.5.

Figure 2.5: GPR processing algorithm. Arrows indicate the progressive
processing steps applied to the GPR data.
71

After importing the data into REFLEX, a 1-D Dewow filter (moving-average
filter in z direction) was applied to all the GPR data profiles. Rubbersheeting was
then used on all the GPR traces to move traces to more accurate locations. The
zero time for the GPR data was chosen to be after the direct air waves,
effectively removing these waves from the data.
A series of 1-D and 2-D filters were applied to the GPR data to reduce
noise. Background filtering was applied to all of the GPR data, before files were
merged. Band pass, deconvolution, and average subtraction filters were then
applied after merging.
After data sets were merged and filtered to isolate good reflectors, a
velocity was assigned to the GPR data to evaluate depths of reflectors. The
average velocity used for all 26 GPR lines was 0.06 m/ns, which is
representative for wet sand (Table 2.1).

Auguring
In the summer of 2008 the sediments across the study area were
characterized using peat and bucket sampling augers. Several of the auger
locations were selected on GPR line intersections to help link the GPR reflectors
to the site stratigraphy (Figure 2.3). Most auguring locations were chosen to
delineate zones where sand was deposited. For this purpose, a peat sampler
with an extension rod, both measuring 1 meter in length, were used. The 1st
meter of the auger was pushed in the ground sampling just the top meter. Notes
were taken on “feeling/sound” for clay or sand, and then these steps would be

72

repeated for the 2nd meter. Any material recovered in the top two meters was
described and sampled. Five detailed boreholes drilled in a northwest-southeast
transect perpendicular to the river and boreholes north of Sportsman’s Park were
used as a control on the spatial variation of sand outside the GPR grid (Figure
2.3). Other detailed borehole logs were conducted to characterize the
stratigraphy in the vicinity of the GPR grid. The contact between the bottom clays
and the sand were difficult to characterize, as the augers disturbed the sediment
to the point that it become difficult to decipher the contact. This uncertainty in the
horizon locations complicates the linkage between the borehole data and the
GPR profiles.

73

RESULTS AND DISCUSSION

Borehole Data
Sand was observed at borehole locations in the northwestern section of
the delta on both sides of the river. The thickness of the sand varied between 1
and 2 meters, and varied between fine-coarse sand, sometimes with abundant
organic material (Figure 2.6). The material below the sand is generally peaty or
mucky clay, likely of lucustrine origin. Silt was observed in some locations but
was not as abundant as the clay. In locations where no sand was found, the
peaty and mucky clay extended from the surface to at least 4 m depth. There
was abundant organic matter in all borehole cores, which mostly started between
0.75-1 meter depths. The organic matter was often present with sand, which
was logged as organic sand. The organic material primarily consisted of reeds
and woody debris, and seeds were present in some samples. Calcareous shells
were present in all the boreholes between 0.8 and 1 m. In boreholes where the
contact existed, the contact between the sand and the clay was often gradational
with finer sediments at depth.

74

Figure 2.6: Borehole data collected in the summer of 2007. All surface sediment
was a very thin (<1cm) layer of organics from current riparian vegetation and was
excluded from this graph. Boreholes are highlighted in Figure 2.3.

75

GPR Signatures
The GPR data was attenuated to the point where only a few significant
reflectors could be identified and interpreted. The degraded data quality could
have been from various reasons such as the antenna being close to the
snowmobile, rapid mode of data acquisition, high organic matter content, and
diffractions from roots of the vegetation. Out of all the GPR data, 5 GPR lines
were chosen, based on data quality, for interpretation (Figure 2.7).
Two types of reflectors are observed in the GPR data. The first can be
characterized as smooth low angle-concave up reflectors, while the second dip
toward the north and northwest. Most of the first types are concave-up reflectors
that occur in the East-West GPR profiles, roughly parallel to the current channel
of the Jordan River and are within 2 m of the ground surface and dips toward the
north and northwest (Figure 2.8). The second type of reflector (Figure 2.9)
occurs at depths between 1 m to 3 m and occurs as smooth low angle-concave
up reflectors.

76

Figure 2.7: A) Raw 2-D GPR data profiles shown with line numbers that correspond to the GPR grid shown in
Figure 2.3. Letters S, C, and M represent sand, clay, and muck respectively, found in cores that coincided with
GPR lines.
77

Figure 2.8: GPR profiles shown with 1 type of interpretation (light blue lines).
78

Figure 2.9: GPR profiles with interpretations of concave-up reflectors (red lines).

79

GPR Interpretation
Some of the sand-clay/muck contacts align with reflectors that indicate a
fluvial sedimentary environment. The interpreted concave-up reflectors (Figure
2.9) can be related to paleo-channels that have been scoured out by the stream
and then filled in with new material (i.e. sand). The dipping reflectors in the EastWest GPR profiles of Figure 2.9 could result from a migrating point bar, with sand
deposition, as the Jordan River has changed positions over time. The dipping
reflectors in the South-North GPR profiles in Figure 2.8 could be formed by the
deposition of sand at higher Lake Charlevoix water levels. The dipping
structures representing foresets created in delta environments when sediment
(i.e. sand) is being transported in the Jordan River is deposited when its
velocities drop at the entrance to Lake Charlevoix.
Current GPR data and borehole data are suggestive of fluvial and deltaic
environments but due to certain limitations the results are not yet conclusive.
Attenuated GPR signals made data interpretation very difficult, so only 5 lines
were used to describe the stratigraphy of the study area. No GPR velocity
measurements were collected, reducing the accuracy of reflector depths. The
snow pack thickness was not accounted for during GPR data processing. Also,
error could have also been introduced when comparing borehole data and GPR
data as the GPR was collected on top of the snow and the boreholes were
augured from the ground surface. During auguring, undisturbed sediment could
not be recovered which introduced significant uncertainty in the elevation link
between borehole data and GPR data.

80

CONCLUSIONS

GPR was effective for visualizing sediment distributions within the study
area. Borehole data showed a gradational boundary between clay and sand,
while GPR reflectors provided insights into the 2D geometry of sediment.
Borehole data aided in the interpretation of the GPR data as either deltaic or
fluvial, but no conclusion could be made about the detailed sedimentary
environment.
Potential future work could include conducting GPR surveys on the other
side of the river where sand was observed in boreholes. GPR data could be
collected at the southern end the study area to see if prominent reflectors are
observed. Marine seismic and sonar bathymetry for offshore geophysics,
combined with GPR on land, would provide more accurate estimates of the total
volume of sand that has been deposited, although dreaded volumes would have
to be included.

81

BIBLIOGRAPHY

82

Baker, P. “Response of Ground-penetrating Radar to Bounding Surfaces and
Lithofacies Variations in Sand Barrier Sequences.” Exploration Geophysics 22
(1991): 19-22.
Baker, G. “Seismic imaging shallower than three meters” Diss. U of Kansas.
1999.
Beres, M., Haeni, F. “Application of ground-penetrating-radar methods in
hydrogeologic studies.” Groundwater 29 (2000): 375-386.
Beres, M., Huggenberger, P., Green, A., Horstmeyer, H. “Using two- and threedimensional georadar methods to characterize glaciofluvial architecture.”
Sedimentary Geology 129 (1999) 1–24.

Bersezio, R., Giudici, M., Mele, M. “Combining sedimentological and
geophysical data for high-resolution 3-D mapping of fluvial architectural elements
in the Quaternary Po plain (Italy).” Sedimentary Geology 202 (2007) 230–248.

Davis, J., Annan, A. “Ground-penetrating radar high-resolution mapping of soil
and rock stratigraphy.” Geophysical Prospecting 37 (1988) 531–551.
Dearing, J., Jones, R. “Coupling temporal and spatial dimensions of global
sediment flux through lake and marine sediment records.” Global and Planetary
Change 39 (2003) 147–168.

Diodato, N. “Modelling net erosion responses to enviroclimatic changes
recorded upon multisecular timescales.” Geomorphology 80 (2006) 164–177.
Foster, G., Dearing, J., Jones, R., Crook, D., Siddle, D., Harvey, A., James, P.,
Appleby, P., Thompson, R., Nicholson, J., Loizeau, J. “Meteorological and land
use controls on past and present hydro-geomorphic processes in the pre-alpine
environment: an integrated lake–catchment study at the Petit Lac d’Annecy,
France.” Hydrological Processes 17 (2003), 3287–3305.
Great Lakes Environmental Research Laboratory. National Oceanic and
Atmospheric Administration. 2010. January 2008
<http://www.glerl.noaa.gov/data/now/wlevels/levels.html>
Jol, H., Smith, D. “Ground Penetrating Radar of Northern Lacustine Deltas.”
Canadian Journal of Earth Science 28 (1991) 1939-1947.

Kasai, M., Brierley, G., Page, M., Marutani, Tomomi, M., Trustrum, N.
“Impacts of land use change on patterns of sediment flux in Weraamaia
catchment, New Zealand.” Catena 64 (2005) 27–60.

83

Knight, R. “Ground Penetrating Radar for Environmental Applications.” Annu.
Rev. Earth Planet. Sci. 29 (2001) 229–55.
Lach, J., Wyzga, B. “Channel Incision and Flow Increase of the Upper Wisloka
River, Southern Poland, Subsequent to the Reafforestation of its Catchment.”
Earth Surface Processes and Landforms 27 (2002) 445-462.
The Lake Charlevoix Association. 2008
<http://www.lakecharlevoixassociation.org/lakelevels.html>.
Liebault, F., Gomez, B., Page, M., Marden, M., Peacook, D., Richard, D.,
“Trotter, C.M. Land-use Change, Sediment Production and Channel Response in
Upland Regions.” River Research and Applications 21 (2005) 739-756.
Michigan Geogrpahic Data Library. State of Michigan. 2010. 1
December 2008 < http://www.mcgi.state.mi.us/MIGDL/>
Nichol, S., Lian, O., Horrocks, M., Goff, J. “Holocene record of gradual,
catastrophic, and human-influenced sedimentation from a backbarrier wetland,
northern New Zealand.” Journal of Coastal Research 23 (2007) 605-617.
Noel, H., Garbolino, E., Brauer, A., Lallier-Verges, E., Beaulieu, J., Disnar, J.
“Human impact and soil erosion during the last 5000 yrs as recorded in lacustrine

sedimentary organic matter at Lac d’Annecy, the French Alps. Journal of
Paleolimnology 25 (2001) 229-244.
Pelpoa, C., Hickin, E. “Long-term bed load transport rate based on aerial-photo
and
ground penetrating radar surveys of fan-delta growth, Coast Mountains, British
Columbia.” Geomorphology 57 (2004) 169–181.
Prothero, D., Schwab, F. Sedimentary Geology, W.H. Freeman and Company,
2004.
Syvitski, J., Milliman, J. “Geology, Geography, and Humans Battle for Dominance
over the Delivery of Fluvial Sediment to the Coastal Ocean.” The Journal of
Geology, 115 (2007) 1–19.
Wallinga, J., Davids, F., Dijkmans. “Luminescence dating of Netherlands’
sediments.” Netherlands Journal of Geosciences 86 (2007) 179-196.

84

CHAPTER 3: EVALUATING THE INFLUENCE OF LAND COVER
ON NUTRIENT TRANSPORT AND IN-STREAM NUTRIENT
PROCESSING

85

INTRODUCTION
Different land uses can leach nutrients into groundwater which then flows
out to streams. Degraded water quality can lead to eutrophication of surface
waters and algal blooms (Heisler et al., 2008), and can cause stress for fish and
other aquatic species (Klingeman et al., 1998). Biogeochemistry is affected as
humans alter the natural landscape, mainly through urbanization and agricultural
development (Fitzpatrick et al, 2007).
It is important to understand how nutrient inputs from land cover are
transported to streams through groundwater fluxes. Nutrient loading into streams
been extensively studied over the last few decades (e.g., Reckhow and Chapra,
1999; and Boutt et al., 2001). Cabana and Ramussen (1996) found an increase
in the delta 15 nitrogen isotope levels of primary consumers as human density
increased; though no specific relationships between land use practices and
nitrogen isotopes in biota could be established. However, positive correlations
have been found between nitrogen isotope levels and increases in percentages
of suburban, and agricultural land cover (Vander Zanden et al., 2005).
Although correlations can help provide information where nutrient sources
come from, the transport mechanisms for nutrients between sources and sinks
need to be understood. Reckhow and Chapra (1999) points out that advances in
models can provide understanding in the transport mechanisms for nutrients,
which need to be combined with observations. Groundwater models and

86

geographic information systems (GIS) have provided additional tools to help
identify sources of nutrient loading into streams (Boutt et al., 2001).
Nitrate is a common groundwater pollutant with fertilizer and septic tanks
being common sources to near-surface aquifers (McMahon et al., 2008). Nitrate
inputs from urban landscapes are commonly believed to be smaller than from
agricultural areas (Nolan and Stoner, 2000). However, outputs from sewage
treatment plants and fertilizers from urban land cover have higher concentrations
than other land covers (Fitzpatrick et al., 2007; Close, 1989). Puckett et al.
(2008) found that short groundwater residence times can be linked to increased
in-stream nitrate concentrations.
Phosphorous transport through groundwater can be an important source
of nutrient loading to streams. The major source of phosphorous in streams is
generally considered to be associated with overland flow from the adjacent land
surface (Carpenter et al., 1998). McCobb et al. (2003) and Chardon et al. (2007)
found TDP concentrations greater than 10 mg/L in sandy soils leached from
manure piles and cow pastures.
More information is needed about loading of phosphorous to streams from
groundwater to reduce water quality degradation. Waschbusch et al. (1995)
conducted a water quality study in two suburban areas and found TDP
concentrations of 0.61 mg/L due to leaching from lawns. Observed values of
TDP leached from manure piles have ranged from 0.80-10 mg/L (McDowel and
Sharpley, 2001). Drainage from septic systems can also provide a significant

87

source of TDP to groundwater as discussed in Arnscheidt (2007) and Fitzpatrick
et al. (2007).
The transport of nitrate and TDP through the groundwater can be
significantly different. In fact, previous research generally considered surface
runoff to be the only contribution of total phosphorous (TP) to streams (Sims et
al., 1998). Tesoriero et al. (2009), across 5 drainage basins in the continental
U.S. (all with poorly drained surface soils),found that nitrate loading into streams
increased proportionally to the base flow index (BFI), while there was no
significant increase in orthophosphate with high values of BFI. However, various
forms of phosphorous have been shown to be transported into shallow
groundwater, especially the forms dissolved phosphorus which is combined
together into TDP in this paper (Sims et al., 1998). Adriano et al. (1975) found
phosphorus levels greater than 0.92 mg/L at 6.6 meters depth in sandy soils.
Numerous studies have reported groundwater concentrations of TP greater than
1.0 mg/L, which exceeded most surface water quality standards (Breesuswma et
al., 1995). Another factor contributing to the transport of TDP in the groundwater
is sediment characteristics. Vanek (1993) concluded that sand higher in calcium
retained more TDP than sand with lower or calcium content. Deep-Siliceous
sands have been found to be a significant source of TDP transported in
groundwater to streams (Sims et al., 1998).
MODFLOW (Harbaugh et al., 2000) and MT3D (A Modular 3-D MultiSpecies Transport Model, Zheng, 2006) have been used to simulate nitrate
concentrations through groundwater with simple non-reactive transport (Almasri

88

and Kaluarchchi, 2007). Lerner and Papatolios (1993) used a simpler approach
to predict groundwater nitrogen concentrations, assuming uniform recharge, to
calculate input concentrations.
The objective of this study is to illustrate how multiple land covers can be
superimposed to identify potential sources of nutrient loading for nitrate and TDP.
This study incorporates the use of groundwater and nutrient transport models
combined with a 1-D description for nutrient uptake for the Sycamore Creek
Watershed in southern Michigan.

Site Description
The different land cover types that Sycamore and Mud Creeks flow
through make the SCW a good study site for testing nutrient inputs from different
land cover types. Current land cover in the Sycamore Creek Watershed (SCW)
varies from forest to agriculture and urban areas. In the SCW, the winter-snow
and fall storms sustain perennial streams via groundwater flow. Coupled with
anthropogenic land use practices and perennial streams fed by groundwater,
nutrient loading from groundwater into streams is an important concern to
investigate. Modeling groundwater transport of nutrients to streams can provide a
powerful tool to improve the understanding of impacts of different land covers to
surface water quality.
Modeling groundwater transport requires knowledge of groundwater flow
based on sediment properties, climate inputs, and sources and sinks of an area.
The Grand River, along with some of its tributaries, makes up the southern and

89

western edges of the model boundary, and the Red Cedar River and its
tributaries form the eastern and northern edges of the model boundary (Figure
3.1). It was important to expand the model boundaries to appropriately define
the area of groundwater capture to the stream system, called the
groundwatershed, as described by Boutt (2001). The groundwatershed (280
km2) and the surface watershed (275 km2) for Sycamore and Mud Creeks are
similar, thus only the watershed will be discussed from this point on.
Lansing

Red Cedar River

Mud Creek

Grand River

Sycamore Creek
Mason

Figure 3.1: The model area along with Sycamore Creek Watershed (SCW,
highlighted in pink) is located in South Central Michigan (from MIGDL).

90

Source and Sinks
Agricultural and urban areas are the two dominant land use types within
the SCW, thus these types represent the two most likely sources for nitrate and
TDP. Agricultural sources of nitrate and TDP include leachate from application of
manure and fertilizer to crops and fields, manure piles from pastures (defined as
Herbaceous Rangeland in the Michigan Geographic Data Library (MIGDL)), and
CFOs (Confined Feeding Operations). Urban sources of nitrate and TDP include
application of fertilizers to lawns and inputs from septic tanks and sewers.
The headwaters of Sycamore and Mud Creeks both start in agricultural
areas and then merge just after passing the city of Mason, which is a small rural
community (Figure 3.2). At the lower end of the drainage basin, Sycamore Creek
flows through highly residential and industrial zones just before flowing into the
Red Cedar River.

91

Figure 3.2: Map of Current Land cover for the Sycamore Creek model area
(MIGDL).

92

METHODS

GIS Data
Modeling groundwater flow and transport across drainage basin scales
requires large data sets that describe the physical nature of the study area such
as topography, aquifer geometry, and geology. A 9-m DEM (Digital Elevation
Model) was used for topography (Figure 3.3) from NED (the National Elevation
Data set), which was created from aerial and topographic maps that represent
bare earth. The DEM provides detailed information about the top of the near
surface aquifer in which nitrate and TDP transport is simulated. The most
important use of the DEM in this study was to assign stream elevations to
Sycamore and Mud Creeks along the entire drainage system.

93

Figure 3.3: A 9 meter resolution DEM plotted with the model area (bold outline)
and the SCW (light outline) (SEAMLESS).
A GIS based Quaternary Geology map developed by Farrand and Bell
(1982) was used as a framework for hydraulic conductivity zonation in the
MODFLOW model (Figure 3.4). This glacial drift aquifer mainly consists of
medium-textured till and end moraines with some glacial outwash deposits.

94

End moraines of
medium-textured till
Glacial outwash sand and gravel
and postglacial alluvium
Medium-textured till

0

km
2.5 5 10

Figure 3.4: Surficial Quaternary Geologic Map used for geologic facies
description for MODFLOW and MT3D (Farrand and Bell, 1982).

Observed groundwater levels, aquifer geometry, and septic tank locations
were created from a single water well coverage (Figure 3.5). Insufficient oil and
gas wells were available in the model area to provide sufficient information. The
well data set is digitized from residential wells in Ingham and Eaton Counties
from MIGDL. All the wells were first filtered on the basis of locations and
elevations. In ArcGIS, the DEM elevation was extracted for each well and then
compared with the values recorded by well drillers. If the absolute difference
between these two values is greater than 1 meter, the well is excluded from
further analysis. Wells that were clearly mislocated (i.e. wells located outside of
Ingham and Eaton Counties despite being labeled as in the county) were also
removed. Wells that had major discrepancies such as negative bedrock
elevations were also removed, as this is out of the range of known elevations.
95

Additional filtering for observed heads was based on casing depth, type of
aquifer, and recorded static water level (SWL). Observations were not used from
wells if the casing depth (bottom of the well) was below the elevation of the
aquifer bottom or if the recorded aquifer material was consolidated rock.
Because most of the wells in the model area are located in the Saginaw
Formation (a major aquifer for Michigan), only 42 out of 300 wells could be used
for water level comparisons with the MODFLOW simulation.

Wells used for simulating
Septic inputs into ground
water and for creating the
bedrock surface
Observatation Wells

0

2.5

5

km
10

Figure 3.5: Wells used to construct the map of the aquifer bottom, for water level
observations, and to provide septic well inputs. Wells used inside the model
boundary (yellow line) were used to simulate nutrient inputs from septic systems
(MIGDL).

96

Most of the residential wells in the model area pump their water from the
Saginaw Aquifer, a major confined aquifer across much of Michigan. The casing
depths of these wells make them suitable to build the bedrock elevation map to
represent the bottom of the glacial drift aquifer that is the focus of this study.
Bedrock wells inside and outside the model boundary were selected, and
bedrock depths for these wells were chosen where the well intersects the top of a
confining formation (i.e. the Marshall Formation or Jurassic red beds that are
shale units). Clear outliers in the data were removed from the data set if the
bedrock elevation varied more than 40-60 meters from surrounding bedrock
elevations. Using this methodology, 41 outliers were removed. The bedrock
elevations were then interpolated using ordinary kriging in ARCGIS using a
circular variogram model (Figure 3.6A) after removing 1st order trend. Variogram
modeling parameters included 164.14 nugget, 100 m lag size, 12 lags, and a
partial sill of 100.46. After removing outliers, the interpolated surface was
converted to a grid with a 20 m cells (Figure 3.6B).

97

Semi-variance*10^-2

3.1
2.48
1.86
1.24
0.62
0

0

150 300 450 600 750 900 1050 1200
Distance, h (m)

A

275 m

200 m

0 2.5 5

km
10

B
Figure 3.6: Variogram model (using a Circular model, 255 m range, 164.14
nugget, 100 sill, 100 m lag size, 12 lags, and no anisotropy) (A) incorporated into
the interpolation of the bedrock elevation map (B).

98

Steady-state Model
The boundary conditions for the groundwater model were assigned from
the DEM, bottom geometry map, stream coverage, and model boundary. The
top elevations for the groundwater model were assigned from the DEM. The
bottom elevations for the ground water model were assigned from the bedrock
elevation map presented in Figure 3.6. Starting heads were initially estimated
using the DEM elevations. After running a preliminary model, the head solution
was exported and used as the starting heads for future models. The streams
were generated from the DEM in Arc by creating flow direction and flow
accumulation grids, which were then used to calculate flow length and stream
order metrics. Once the artificial stream network was created, the elevations
from the DEM were used for stream elevations. Along with drain elevations,
conductance values (L2/T) were assigned to each arc. The conductance was
calculated as the stream width multiplied by the hydraulic conductivity of stream
bottom sediments times the stream segment length, divided by the thickness of
the stream bottom sediments. GMS automatically calculated the stream lengths
using the information mapped over from the GIS streams coverage. An average
width of 10 m was estimated from field data. A representative hydraulic
conductivity of 2 m/d was used for the sediment hydraulic conductivity for
Sycamore and Mud Creeks. The thickness of the bottom sediments was
assumed to be 1 meter for all drain cells.
Another important step in the development of a groundwater model is
assignment of recharge across the model area. An initial recharge value of 22

99

cm/yr was estimated based on stream flow data from a USGS gage station at
Holt road, which was removed in 1997. The minimum recharge estimate was
calculated by dividing base flow measured in August 2007 by the area of the
groundwatershed for the USGS gage. Leakage between the near surface glacial
aquifer and the underlying Saginaw aquifer was not calculated because most of
the model region has a significant confining bed thickness between the aquifers.
Since no hydraulic conductivity data are available from pump tests and
slug tests within the model area, a single value of hydraulic conductivity (5 m/d)
was used as a starting value for all of the geologic zones (Freeze and Cherry,
1979; Schwartz and Zhang, 2003). The hydraulic conductivities for the main
geologic zones were then optimized using the PEST parameter estimation code.
The measured SWLs from the GIS well coverage were used as observed heads
to be compared to simulated heads from MODFLOW (Figure 3.5). One hydraulic
conductivity value for the model area was optimized to obtain the minimum
residual error between the simulated and observed heads. Parameter
optimizations for two or three geologic zones did not result in a significant
reduction in residuals, thus only one value hydraulic conductivity, optimized to
6.44 m/d, was used for all glacial deposits. Figure 3.7 shows observed vs.
simulated heads along with a spatial plot of the final head solution. Simulated
stream flow was extracted from the model results using a script and compared
against stream flows that we measured at or near base flow. Only 1 measured
flow in this study, located at Aurelius, matched well (Table 3.1). The Aurelius site

100

matched well because the recharge for this site was calculated based on the
baseflow for this station.

Table 3.1: Simulated flow and averaged measured flow at three locations

Simulated Head (m)

295

285

275

265
y = 0.946x + 14.91

255

245
245

2

R = 0.9281

255

265

275

285

295

Observed Head (m)

Figure 3.7: Simulated heads plotted relative to observed heads (A) and mapped
(B).

101

Head (m)
295
250

0 2.5 5

km
10

Figure 3.7 Continued

Groundwater transport using MT3D
Possible sources of nutrients observed during base flow are septic tank
inputs and leaching from agriculture, pasture, and urban landscapes. We
examined the potential inputs of each of these main potential sources by
simulating each land use derived source in a separate transport models then
superimposing these into one solution. This allows the influence of each land
use on water quality to be examined.

102

The benefit of simulating separate land cover types is that solutions from
the transport models can be superimposed to evaluate the effects of multiple
land uses and to estimate input concentrations, using optimization. For each
dominant land cover, unit concentrations (1 mg/L) were assigned for both Nitrate
and TDP. Residential wells within the model boundary were selected to be
septic sources (Figure 3.5). It was assumed that a septic input existed for each
residential water well because the residential water wells exist in rural areas
where city sewer systems did not exist. The unit concentration for the septic
tanks was applied to a volumetric flux of 3 m3/day, which was the estimated
water usage per household of four people (Elmwood, 2008). Polygons of
agriculture, herbaceous pasture, and urban land covers were used to assign the
starting concentrations for Nitrate and TDP. The concentrations were applied to
the recharge flux (22 cm/yr) as recharge concentrations.
When simulating groundwater transport on any scale, values for
dispersion and porosity are usually unavailable, unless derived from empirical
measures, thus representative values were chosen. The ratios of transverse to
longitudinal dispersion were selected as 500, 1E-4 m2/d for the molecular
diffusion coefficient, and 5 m was used as for the longitudinal dispersitivty
(Freeze and Cherry, 1979). A representative value for porosity of 0.3 was
assigned for all geologic zones.
Nitrate and TDP in this study are assumed to be mobile, thus the nutrient
transport is mainly by advection. The advection part of the transport equation
used by MT3D is solved using the Hybrid Method of Characteristics (HMOC)

103

which combines the Modified Method of Characteristics (a backward particle
tracking scheme) and the Method of Characteristics (a scheme that tracks
multiple particles forward in time). It is important to note that assuming mobile
species is a significant simplification for simulating groundwater fluxes to
streams, especially for the transport of TDP which can bind to iron oxides on
sediments.
Observed mass fluxes and total nutrient concentrations for Sycamore and
Mud Creeks were calculated using measured stream flow and in-stream
concentrations (Figure 3.8) during base flow conditions in the summer of 2005,
as described in Welty (2005). After starting concentrations and values for the
dispersion and advection packages were assigned, the four transport simulations
for septic inputs, agriculture, pasture, and urban landscapes were run for 150
years to account for the main period of anthropogenic settlement (Dickmann and
Leefers, 2004), based on the steady-state flows from MODFLOW. Water fluxes
from MODFLOW were multiplied by concentrations from MT3D to calculate
nutrient mass fluxes on a cell by cell basis. Measured and simulated
concentrations (mg/L) and stream flows (m3/s) for each site were multiplied
together to calculate the mass fluxes (mg/s).

104

2000.0

2000.0
1800.0
1600.0
1400.0
1200.0
1000.0
800.0
600.0
400.0
200.0
0.0

1500.0
1000.0
500.0
0.0

Nitrate Differential
Flux (mg/s)

Nitrate Total Flux
(mg/s)

Sycamore Creek

-500.0
0.592 3.287 11.701
17.235 26.34 29.226 33.15

Distance Downstream (km)
Total Nitrate Flux

Differential Nitrate Flux

Mud Creek
80.0

100.0

60.0

80.0

40.0

60.0

20.0

40.0

0.0

20.0
0.0

-20.0
0.836

A

6.34

10.238

14.512

19.946

Distance Downstream (km)

Figure 3.8: A) Nitrate observations, B) TP observations (Welty, 2005).

105

Nitrate Differential
Flux (mg/s)

100.0

120.0

Nitrate Total Flux
(mg/s)

140.0

Sycamore Creek
40.0

80.0

30.0

60.0
20.0

50.0
40.0

10.0

30.0

0.0

20.0
-10.0

10.0
0.0

TP Differential
Flux (mg/s)

TP Total Flux
(mg/s)

70.0

-20.0
0.592 3.287 11.701 17.235 26.34 29.226 33.15

Distance Downstream (km)
Total TP Flux

Differential TP Flux

Mud Creek
7.0

8.0

6.0

7.0

5.0

6.0

4.0

5.0

3.0

4.0

2.0

3.0
2.0

1.0

1.0

0.0

0.0

-1.0
0.836

B

6.34

10.238

14.512

19.946

Distance Downstream (km)

Figure 3.8 Continued

106

TP Differential
Flux (mg/s)

8.0

9.0

TP Total Flux
(mg/s)

10.0

Optimization using a simple 1-D In-stream Processing Description
This study assumes that nutrient loading into streams from groundwater can
be linearly superimposed based on separate land cover based subsurface inputs.
Equation 3.1 was used to minimize the sum of absolute residuals between
measured and simulated fluxes by adjusting scaling coefficients for each land
use type.
MIN{SUM(ABS((aAg+bU+cS+dR-rL+MWTP)-Obs))}

(3.1)

where
a, b, c, d = coefficients representing starting concentrations to be estimated
(mg/L).
Ag, U, S, and R = simulated unit fluxes (total or differential) from Agriculture,
Urban, Septic, and Pasture (mg/s).
r = processing rate of nutrients by length of stream along palustrine wetlands
(mg/s/m)
L = length of palustrine wetlands connected to the stream (m).
MWTP = Calibrated mass flux for the Mason Waste Water Treatment Plant
(mg/s)

This superposition approach could be used with multiple land use types
due to linearity of the concentrations for a stream that effectively accumulates
flow and associated solutes from multiple land use sources. Excel Solver was
used for the optimizations using equation 3.1. A description of in-stream
processing is included to account for the observed negative differential nutrient
concentrations (Figure 3.8), which is assumed to be due to wetland processing,
107

including biological and physical processing (Angier, 2002). The in-stream
processing was described using the measured stream length that flows through
palustrine wetlands upstream of each point where in-stream chemistry was
measured. The stream lengths (in meters) were then multiplied by a processing
rates (mg/s/m) which then results in a 1-D in-stream nutrient processing model to
account for the observed negative differential solute fluxes between downstream
observation locations.
Optimizations for the final nutrient concentrations included a calibrated
value for the Mason Waste Water Treatment Plant (MWTP). Using a direct mass
flux was used to help improve the optimizations because a large bias within data
existed due to large values of nutrient mass flux downstream of Jefferson. The
location of the MWTP is between Jefferson and Harper observation sites, which
explain the large increase for both TP and nitrate.

108

RESULTS AND DISCUSSION

Optimization Results
The optimization of simulated fluxes from MT3D along with in-stream
processing rates, helped explain the negative fluxes observed in the field nutrient
data. Along with optimizing input concentrations for land cover, in-stream
processing rates were also optimized for nitrate and TDP. Adding the in-stream
processing rate to the TDP optimization did not improve matches between
simulated and observed TP fluxes. However, optimizing the in-stream
processing rate significantly improved the matches between simulated and
observed nitrate fluxes.
Using values for TP and Nitrate inputs from the MWTP improved the
optimization results. After numerous optimizations of nitrate and TDP simulated
fluxes while just using in-stream processing, nutrient loading from the MWTP was
added for both nutrients. Nutrient flux, for both nitrate and TDP, was first
estimated for MWTP. After additional optimizations of input concentrations, instream processing rates, and nutrient loading rates from the MWTP for both
nitrate and TDP, the summed error between the simulated fluxes and observed
fluxes was reduced significantly (Figure 3.9). From the optimizations it is
estimated that the MWTP inputs a total flux of 1950 mg/s for nitrate and a total
flux of 98 mg/s for TP on average.
During the optimizations of nitrate and TDP simulated fluxes, four
optimization scenarios were run using 1) Observed and Simulated Total
109

Concentration; 2) Observed and Simulated Differential Concentration; 3)
Observed and Simulated Total Flux; and 4) Observed and Simulated Differential
Flux. Overall, the optimizations of differential and total nutrient flux provided
smaller summed error than using either total or differential nutrient
concentrations. Thus, only results from differential and total nutrient flux are
presented in Figure 3.9. Optimization of total nutrient flux had better R-squared
and thus the optimized input concentrations, using total nutrient flux, are used in
the final transport model.

Final MT3D Model Results
Simulations of Nitrate and TDP provided a quantitative ability to assess
loadings from each land use type for Sycamore and Mud Creeks. Using the final
concentrations from the optimizations with the simulated groundwater fluxes,
mass flux from each land use in Sycamore and Mud Creeks were calculated for
each nutrient. Mass fluxes from groundwater were calculated by entering the
optimized concentrations into MT3D and running the transport model to simulate
the final in-stream concentrations and mass fluxes for Sycamore and Mud
Creeks. These simulations indicate that urban areas are likely the primary
source for TP loading while Rangeland is the main source of nitrate into
Sycamore and Mud Creeks (Figures 3.10 to 3.14). Agriculture appears to be the
secondary source of loading for nitrate and TP in the SCW, however, there is
significant uncertainty associated with this flux as agriculture was insensitive in
the optimizations.

110

Simulated NOx (mg/s)

Nitrate Total Flux
2500
2000
1500
1000

y = 1.0099x

500

R = 0.8711

2

0

A

0

500
1000
1500
Observed NOx (mg/s)

2000

Simulated NOx (mg/s)

Nitrate Differential Flux

B

2000
1500
1000
y = 1.0132x
2
R = 0.8485

500
0
-500
-500

0

500

1000

1500

2000

Observed NOx (mg/s)
Figure 3.9: Optimization results for Simulated and Observed total and differential
fluxes for Nitrate (A, B) and TDP (C, D) with in-stream processing coefficients.

111

Simulated TDP (mg/s)

TDP Total Flux
120
100

y = 1.0828x
2
R = 0.4512

80
60
40
20
0
0

C

20

40

60

80

Observed TP (mg/s)

TP Differential Flux
Simulated TDP
(mg/s)

120
100
80
60
40
20
0
-20
-40

D

y = 1.0796x
2
R = 0.2472

-20

-10

0

10

20

Observed TP (mg/s)
Figure 3.9 Continued

112

30

40

Figure 3.10: Simulations based on the final optimized nutrient concentrations
from Agriculture, Pasture, Septic Tanks, and urban areas.
113

Sycamore Creek
80
60
50
40
30
20

Total Flux (mg/s)

70

10
0.6

3.3

11.7

17.2

26.3

29.2

0
33.2

Distance downstream (km)

Urban

Agriculture

Pasture

Septic

Observed Total Flux

Mud Creek
12

8
6
4

Total Flux (mg/s)

10

2

0.8

6.3

10.2

14.5

0
19.9

Distance downstream (km)

Figure 3.11: Simulated and Observed TDP Flux from different land cover for
Sycamore and Mud Creeks.

114

0.6

3.3

11.7

17.2

26.3

29.2

2000
1800
1600
1400
1200
1000
800
600
400
200
0
33.2

Total Flux (mg/s)

Sycamore Creek

Distance downstream (km)

Urban

Agriculture

Pasture

Septic

Observed Total Flux

Mud Creek
600
500
400
300
200

Total Flux (mg/s)

700

100
0.8

6.3

10.2

14.5

0
19.9

Distance downstream (km)

Figure 3.12: Simulated and Observed Nitrate Flux from different land cover for
Sycamore and Mud Creeks.

115

0.6

3.3

11.7

17.2

26.3

29.2

0.1
0.09
0.08
0.07
0.06
0.05
0.04
0.03
0.02
0.01
0
33.2

Concentration (mg/L)

Sycamore Creek

Distance downstream (km)

Agriculture

Pasture

Septic

Observed Concentration

Mud Creek

0.06
0.05
0.04
0.03
0.02
0.01

0.8

6.3

10.2

14.5

Concentration (mg/L)

Urban

0
19.9

Distance downstream (km)

Figure 3.13: Simulated and Observed TDP Concentration from different land
covers for Sycamore and Mud Creeks.

116

Sycamore Creek
4
3
2.5
2
1.5
1
0.5
0.6

3.3

11.7

17.2

26.3

29.2

Concentration (mg/L)

3.5

0
33.2

Distance downstream (km)

Urban

Agriculture Pasture

Observed Concentration

Septic

Mud Creek
2.5
2
1.5
1
0.5

0.8

6.3

10.2

14.5

Concentration (mg/L)

3

0
19.9

Distance downstream (km)

Figure 3.14: Simulated and Observed Nitrate Concentrations from different land
cover for Sycamore and Mud Creeks.

117

CONCLUSIONS

The results from this approach show that both nitrate and TDP fluxes can be
reasonably simulated using a simple method with linear superposition of inputs
from different land use types and septic inputs, with an added linear processing
term in wetlands along the streams.
Simulated fluxes from MT3D compared reasonably with in-stream nutrient
measurements at low base flow. The optimized total nitrate fluxes suggest that
pasture accounts for a large percent of the nitrate loading into streams at low
base flow. It also appears that TDP from urban areas are significant sources of
TP at low base flow.
Unit concentrations, 1 unit concentration (mg) per unit volume (L) or 1
mg/L, assigned in MT3D, combined with the superposition (summed simulated
concentrations) of all the land covers, enabled a simple optimization model that
can be applied to other eco-hydrology studies. Optimization of NOx and TDP
differential and total fluxes provided valuable insight on how multiple land covers
can impact nutrient loading to streams via groundwater transport. Describing the
in-stream processing with only stream length flowing through wetlands helped
explain negative differential fluxes observed within the measured in-stream
nutrient data for nitrate. Also, using a calibrated mass flux for the Mason Water
Treatment Plant increased R-squared values for optimizations of both nutrients
and helped explain large increases in mass fluxes for both nutrients downstream
of the city of Mason. Future methods will include constructing a transient
118

transport model for nitrate and TDP to assess how the loading of these nutrients
has changed along with the anthropogenic footprint on the landscape.

Acknowledgments
We would like to thank National Science Foundation and the Environmental
Protection Agency for providing funding for this project. We also thank Anthony
Kendall for developing the script and provided advice over the course of this
project. Thanks also go to Dr. Jan Stevenson, Seth Hunt, and Nick Welty for
their contributions to the data collection and sample analysis.

119

BIBLIOGRAPHY

120

Adriano, D., Novak, L., Erickson, A. Wolcott, A., Ellis, B. “Effect of long term land
disposal by spray irrigation of food processing wastes on some chemical
properties of the soil and subsurface water. Journal of Environmental Quality 4
(1975). 242-248.
Almasri, M., Kaluarachchi, J. “Modeling nitrate contamination of groundwater in
agricultural drainage basins.” Journal of Hydrology 343 (2007) 211-229.
Arnscheidt, J., Jordan, P., Li, S., McCormick, S., McFaul, R., McGrogan, H.,
Neal, M., Sims, J. “Defining the sources of low-flow phosphorous transfers in
complex catchments.” Science of the Total Environment 382 (2007) 1-13.
Breeuwsma, N., Reijerink, J., Schoumans, O. 1995. “Impact of manure on
accumulation and leaching of phosphate in areas of intensive livestock farming.”
In “Animal waste and the land-water interface, K. Steele (ed.)” Lewis Publishers,
CRC Press, New York, NY Pg. 239-251
Boutt, D., Hyndman, D., Pijanowski, B., Long, D. “Indentifying potential land usederived solute sources to stream base flow using groundwater models and GIS.”
Groundwater 9 (2001) 24-34.
Cabana, G., Rasmussen, J. “Comparison of aquatic food chains using nitrogen
isotopes.” Ecology 93 (1996) 10844-10847.
Carpernter, S., Caraco, N., Correll, D., Howarth, R., Sharpley, A., Smith, V.
“Nonpoint pollution of surface waters with phosphorous and nitrogen.” Ecological
Applications 8 (1998) 559-568.
Chardon, W., Aalderink, G., van der Salm, C. “Phosphorous leaching from cow
manure patches on soil columns.” Journal of Environmental Quality 36 (2007)
17-22.
Close, M. “Effects of septic tank effluent on chemical quality of alluvial gravel
aquifers.” New Zealand Journal of Marine and Freshwater Research 23 (1989)
275-286.
Duff, J., Tersoriero, A., Richardson, W. “Whole-stream Response to Nitrate
Loading in Three Streams Draining Agricultural Landscapes.” Journal of
Environmental Quality 37 (2008) 1133-1144.
Farrand, W.R., and Bell, D.L., Quaternary geology of Michigan. Michigan
Department of
Natural Resources, Geological Survey Division map, Lansing. 1982.

121

Fitzpatrick, M., Long, D., Pijanowski, B. “Exploring the effects of urban and
agricultural land use on surface water chemistry, across a regional drainage
basin, using multivariate statistics.” Applied Geochemistry 22 (2007) 1825-1840.
Freeze, R., Cherry, J. “Groundwater.” Prentice Hall, Inc. Upper Saddle River, NJ.
(1979).
Harbaugh, A.W., Banta, E.R., Hill, M.C., and McDonald, M.G. MODFLOW-2000,
the U.S. Geological Survey modular ground-water model -- User guide to
modularization concepts and the Ground-Water Flow Process: U.S. Geological
Survey Open-File Report 00-92. 2000
Heisler, J., Glibert, P., Burkholder, J., Anderson, D., Cochlan, W., Dennison, W.,
Dortch, Q., Gobler, C., Heil, C., Humphries, E., Lewitus, A., Magnien, R.,
Marshall, H., Sellner, K., Stockwell, D., Suddleson, M. “Eutrophication and
harmful algal blooms: A scientific consensus.” Harmful Algae 8 (2008) 3-13.
Jayawickreme, D. H., and D. W. Hyndman. “Evaluating the Influence of Land
Cover on
Seasonal Water Budgets using NEXRAD and Stream flow Data.” Water
Resources
Research 43 (2007).
Kent, D., Wilkie, J., Davis, J. “ Modeling the Movement of a pH Perturbation and
its Impact on Absorbed Zinc and Phosphate in a Wastewater-contaminated
Aquifer.” Water Resources Research 43 (2007).
Klingeman, Peter C., Beschta, Robert L., Komar, Paul D., Bradley, and Jeffery
B.,
Gravel-Bed Rivers in the Environment, Water Resources Publications, LLC.,
1998.
Lerner, D., Papatolios, K. “A simple analytical approach for predicting nitrate
concentrations in pumped groundwater.” Groundwater 31 (1993) 370-375.
Michigan Geographic Data Library. State of Michigan. 2010. 1
December 2008 < http://www.mcgi.state.mi.us/MIGDL/>
McCobb, T., LeBlanc, D., Walter, D., Hess, K., Kent, D., Smith, R. “Phosphorus
in a groundwater contaminant plume discharging to Ashumet Pond, Cape Cod,
Massachusetts, 1999.” Water-Resources Investigations Report 02-4306, US
Geological Survey (2003).
McDowell, R., Sharpley, A. “Soil phosphorous fractions in solution: influence of
fertilizer and manure, filtration and method of determination.” Chemosphere 45
(2001) 737-748.
122

Mulholland, P., Helton, A., Poole, G., Hall, R., Hamilton, S., Peterson, B., Tank,
J., Ashkenas, L., Cooper, L., Dahm, C., Dodds, W., Findlay, S., Gregory, S.,
Grimm, N., Johnson, S., McDowell, W., Meyer, J., Valett, H., Webster, J.,
Arango, C., Beaulieu, J., Bernot, M., Burgin, A., Crenshaw, C., Johnson, L.,
Niederlehner, B., O’Brien, J., Potter, J., Sheibley, R., Sobota, D., Thomas, S.
“Stream denitrification across biomes and its response to anthropogenic nitrate
loading.” Nature 452 (2008) 202-206.
Mullaney, J. “Summary of Water Quality Trends in the Connecticut River, 19681998.” American Fisheries Society Monograph 9 (2004) 273-286.
National Map SEAMLESS Server. United States Geological Survey. 2010. July
2007.
< http://seamless.usgs.gov/>
Nolan, B., Stoner, J. “Nutrients in groundwaters of the conterminous United
States 1992-1995.” Environmental Science and Technology 34 (2000) 11561165.
Puckett, L., Zamora, C., Essaid, H., Wilson, J., Johnson, H., Brayton, M., Vogel,
J. “Transport and fate of nitrate at the ground-water/surface-water interface.”
Journal of Environmental Quality 37 (2008) 1034-1050.
Reckhow, K., Chapra, S. “Modeling excessive nutrient loading in the
environment.” Environmental Pollution 100 (1999) 197-207.
Schwartz, F.W., Zhang, H., “Fundamentals of Groundwater.” John Wiley and
Sons, Inc. New York, NY. (2003).
Sims, J., Simard, R., Joern, B. “Phosphorous loss in agricultural drainage:
Historical perspective and current research.” Journal of Environmental Quality
27 (1998). 277-293.
Tesoriero, A., Duff, J., Wolock, D., Spahr, N., Almendinger, J. “Identifying
pathways and processes affecting nitrate and orthophosphate inputs to streams
in agricultural watersheds.” Journal of Environmental Quality 38 (2009). 18921900.
Vander Zanden, M., Vadeboncouer, Y., Diebel, M., Jeppesen, E. “Primary
consumer stable nitrogen isotopes as indicators of nutrient source.”
Environmental Science Technology 39 (2005). 7509-7515.
Vanek, V. “Transport of Groundwater-borne Phosphorous to Lake Bysjon, South
Sweden.” Hydrobiology 251 (1993) 211-216.

123

Village of Elmwood Park: Water Usage Facts, 2008.
http://www.elmwoodpark.org/water/Facts.htm
Waschbusch, R., Selbig, W., Bannerman, R. “Sources of phosphorus in
stormwater and street dirt from two urban residential basins in Madison,
Wisconsin, 1994-95.” U.S. Geological Survey Water Resources Investigations
Rep. 99-4021. USGS, Washington, DC
Welty, N. Exploring relationships between land use and ecohydrology using
multivariate statistics and process-based models. Diss. Michigan State
University. 2005.
Zheng, Chunmiao. MT3DMS v5.2 Supplemental User's Guide, Technical Report
to the U.S. Army Engineer Research and Development Center. 2006.

124