CIRCULATION AND EXCHANGE IN THE SAGINAW BAY-LAKE HURON SYSTEM:
OBSERVATIONS AND NUMERICAL MODELING
By
Tuan Duc Nguyen

A DISSERTATION
Submitted to
Michigan State University
in partial fulfillment of the requirements
for the degree of
Civil Engineering - Doctor of Philosophy
2014

ABSTRACT
CIRCULATION AND EXCHANGE IN THE SAGINAW BAY-LAKE HURON SYSTEM:
OBSERVATIONS AND NUMERICAL MODELING
By
Tuan Duc Nguyen
Knowledge of lake circulation is essential for addressing many issues ranging from water quality
to human and ecosystem health. Lake Huron, the third largest of the Great Lakes by volume, has
been significantly affected by natural and anthropogenic activities. Since Lake Huron is a
connecting waterway between the upper and lower Great Lakes, understanding Lake Huron
circulation and thermal structure is also important for questions involving the lower lakes.
In this study, we use a three-dimensional, unstructured grid hydrodynamic model to examine
circulation, thermal structure, ice cover extent, and exchange in the Saginaw Bay - Lake Huron
system during summer months for 3 consecutive years (2009-2011) and winter months for 2
years (2010 and 2013). The model was tested against ADCP observations of currents, data from
a Lagrangian drifter experiment in the Saginaw Bay, observations of temperature from
thermistor chains, and temperature data from the National Data Buoy Center stations. Mean
circulation was predominantly cyclonic in the main basin of Lake Huron with current speeds in
the surface layer being highest in August in summer and in January in winter. Circulation in the
Saginaw Bay was characterized by the presence of an anti-cyclonic gyre at the mouth of the
outer bay and two recirculating cells within the inner bay for both seasons. The ice cover data
extracted from Moderate Resolution Imaging Spectroradiometer (MODIS) with relatively high
spatial resolution and from the Great Lakes Ice Atlas were used to test against results obtained
from an ice model. The results show that the ice model was able to simulate lake circulation and
ice cover extent in winter season reasonably well. The percent coverage of ice reached a

maximum of 38.3% and 38.7% in mid-February in 2010 and beginning of March in 2013
respectively. New estimates are provided for the mean flushing times (computed as the volume
of the bay divided by the rate of inflow) and residence times (computed as e-folding flushing
times treating the bay as a continuously stirred tank reactor) for Saginaw Bay for summer and
winter seasons. The average flushing time (over the three months of summer and for all three
years) was 23.0 days for the inner bay and 9.9 days for the entire bay. The corresponding values
for the winter season are 43.2 days and 15.6 days respectively. The mean e-folding flushing time
was 62 days for summer and 64.7 days for winter for the inner bay and 115 days for the summer
and 114.2 days for the winter conditions for the entire bay. Empirical relations between the mean
residence time and river discharge were proposed. To characterize the behavior of river plumes
in the inner Saginaw Bay, the absolute diffusivity values in the along-shore and cross-shore
directions were calculated using data from GPS-enabled Lagrangian drifters and simulation
results based on particle transport models.

To my parents, my wife, and my lovely children

iv

ACKNOWLEDGEMENTS

First of all, I would like to thank my advisor Dr. Mantha S. Phanikumar for his guidance and
support throughout my study and this dissertation. Without his invaluable advice, patience,
encouragement and the time he spent the dissertation would not have been possible.
A special note of thanks to my committee members, Dr. Roger B. Wallace, Dr. Shu-Guang Li,
and Dr. Farhad Jaberi for serving on my the dissertation committee and for their assistance and
valuable comments.
I am also thankful to Dr. Pramod Thupaki for his suggestions throughout and numerous
discussions during my research. I would like to extend my appreciation to Drs. David J. Schwab
and Eric J. Anderson from GLERL for their support and for giving me the opportunity to work in
their laboratory during Summer 2011. Without them I would have been challenged in collecting
observational data. I would like to thank my friend and colleague Ms. Guoting Kang for her help
in processing the ice data.
I am also thankful for 4 years of scholarship from the project “Agricultural Scientific
Technology” funded by the Vietnam Ministry of Agricultural and Rural Development (MARD).
My gratitude also extends to my Vietnamese friends who have contributed to making my time at
MSU memorable.
I owe a lot to my parents for the love and constant support they provided during the entire
process. Last but not least, I would also like to thank my wife, Oanh Tran, my daughter, Ngan
Giang Nguyen, and my son, Tuan Anh Nguyen, for being wonderfully supportive and
encouraging throughout this process.

v

TABLE OF CONTENTS
LIST OF TABLES .................................................................................................................... viii
LIST OF FIGURES ..................................................................................................................... x
Chapter 1 ...................................................................................................................................... 1
Introduction .................................................................................................................................. 1
1.1 Problem description ........................................................................................................... 1
1.2 Summer Circulation and Thermal Structure ...................................................................... 5
1.3 Winter Circulation and Ice Cover Extents ......................................................................... 8
1.4 Residence and Flushing Times .......................................................................................... 9
1.5 Goals of the Thesis........................................................................................................... 11
1.6 Structure of the Thesis ..................................................................................................... 11
Chapter 2 ....................................................................................................................................
Field Observations and Analysis ..............................................................................................
2.1 Observational Data...........................................................................................................
2.1.1
Currents ..................................................................................................................
2.1.2
Temperature ...........................................................................................................
2.1.3
Lagrangian Drifters ................................................................................................
2.2 Numerical Model .............................................................................................................
2.3 Uncertainty of measurements ..........................................................................................

13
13
13
14
15
16
20
23

Chapter 3 .................................................................................................................................... 24
Summer Circulation and Thermal Structure ......................................................................... 24
3.1 Introduction ...................................................................................................................... 24
3.2 Numerical model .............................................................................................................. 25
3.2.1
Model configuration............................................................................................... 25
3.2.2
Meteorological Data............................................................................................... 31
3.3 Results and Discussions ................................................................................................... 43
3.3.1
Results .................................................................................................................... 43
3.3.1.1
Currents .......................................................................................................... 43
3.3.1.2
Temperature .................................................................................................... 49
3.3.2
Discussion .............................................................................................................. 56
3.4.3
Particle Tracking Model ........................................................................................ 63
3.4.4
Dispersion coefficient in Saginaw Bay .................................................................. 66
3.4.5
Lagrangian time-scale and length-scale ................................................................. 71
3.4 Conclusions ...................................................................................................................... 76
Chapter 4 ....................................................................................................................................
Ice Cover, Winter Circulation and Thermal Structure .........................................................
4.1 Introduction ......................................................................................................................
4.2 Observations and Numerical Methods .............................................................................
4.2.1
Numerical Model ...................................................................................................
vi

77
77
77
80
80

4.2.2
Meteorological and observational data .................................................................. 84
4.3 Results and Discussion .................................................................................................... 86
4.3.1
Results .................................................................................................................... 86
4.3.1.1
Currents .......................................................................................................... 86
4.3.1.2
Temperature .................................................................................................... 88
4.3.1.3 Ice Cover, Ice Concentration, and Ice Thickness ................................................ 97
4.3.2
Discussion ............................................................................................................ 108
4.4 Conclusions .................................................................................................................... 120
Chapter 5 ..................................................................................................................................
Residence Time in the Saginaw Bay and Lake Huron System ............................................
5.1 Introduction ....................................................................................................................
5.2 Methodology ..................................................................................................................
5.2.1
First-Order Transport Time Scales ......................................................................
5.2.2
Eulerian Method...................................................................................................
5.3 Results and Discussion ..................................................................................................
5.3.1
Results ..................................................................................................................
5.3.1.1
First-Order Transport Time Scales ...............................................................
5.3.1.2
The Eulerian Approach .................................................................................
5.3.2
Discussion ............................................................................................................
5.4 Conclusions ....................................................................................................................

122
122
122
124
124
125
126
126
126
136
146
154

Chapter 6 .................................................................................................................................. 155
Conclusions ............................................................................................................................... 155
BIBLIOGRAPHY .................................................................................................................... 158

vii

LIST OF TABLES
Table 1: Details of ADCP deployments in Lake Huron............................................................... 15
Table 2 : Details of Thermistor deployments in Lake Huron ...................................................... 16
Table 3: Details of Lagrangian drifter deployments in Saginaw Bay .......................................... 18
Table 4: Accuracy of the measurements ...................................................................................... 23
Table 5: NDBC Meteorological station locations ........................................................................ 34
Table 6: GLOS Meteorological station locations ......................................................................... 35
Table 7: NCDC Meteorological station locations ........................................................................ 38
Table 8: RMSE and standard deviations between ADCP observed current values and results
from the Saginaw Bay model (SBM)............................................................................................ 45
Table 9: Maximum values of absolute dispersion rates (in m2/s) for Lagrangian drifters and
simulated particle paths in Saginaw Bay ...................................................................................... 69
Table 10: Lagrangian Statistics for Drifters in Saginaw Bay....................................................... 73
Table 11: Maximum Ice Cover Area over Lake Huron (unit: %) ............................................... 98
Table 12: Monthly average ice thickness (in cm) for Lake Huron and Saginaw Bay in winter
months ......................................................................................................................................... 105
Table 13: Maximum ice thickness (in cm) for Lake Huron and Saginaw Bay in winter months
..................................................................................................................................................... 105
Table 14: Maximum and mean currents velocity in Lake Huron in winter months (unit: cm/s) 110
Table 15: Mean flushing time 𝑻𝒇 (in days) for inner and entire Saginaw Bay in summer months
..................................................................................................................................................... 129
Table 16: Mean flushing time 𝑻𝒇 (in days) for inner and entire Saginaw Bay in winter months
..................................................................................................................................................... 134
Table 17: Mean residence time (e-folding flushing time) in days for Inner and Entire Saginaw
Bay in summer season ................................................................................................................ 140
viii

Table 18: Mean residence time (e-folding flushing time) in days for Inner and Entire Saginaw
Bay in winter season (Without Ice Model) ................................................................................. 143
Table 19: Statistical data for daily discharge in m3/s ................................................................. 151

ix

LIST OF FIGURES

Figure 1: Map of Lake Huron showing the important physical features: Manitoulin Basin; Bay
City Basin; NC: North Channel; Alpena- Amberley Ridge; Georgian Bay. (source:
http://www.ngdc.noaa.gov/)............................................................................................................ 4
Figure 2: Aerial view of the Saginaw Bay area showing the upstream region along the Saginaw
River, the river mouth, and the Shelter Island (Photo credit : Doc Searls) .................................... 6
Figure 3: (a) Map of Lake Huron showing the important physical features. MB: Manitoulin
Basin; BCB: Bay City Basin; NC: North Channel; AAR: Alpena- Amberley Ridge. (b) Map of
Saginaw Bay showing the ADCP locations. SI: Shelter Island. 1-1 and 2-2 in Figure 3(a) are the
inner and outer Saginaw Bay transects respectively and 3-3 and 4-4 are transects for which
vertical temperature profiles are shown. ....................................................................................... 12
Figure 4: Photograph showing (a) the Lagrangian drifters, (b) a drifter in water, and (c) the
Trackpack® GPS .......................................................................................................................... 19
Figure 5: Map of Lake Huron showing the bathymetry interpolated from data downloaded from
the NOAA National Geophysical Data Center, Geophysical Data System (GEODAS). ............. 27
Figure 6: A histogram of model grid size of Lake Huron Model (LWM) .................................. 28
Figure 7: A histogram of model grid size of the Saginaw Bay Model (SBM) ............................ 29
Figure 8: The computational mesh used for the Lake Huron model (LWM) .............................. 30
Figure 9: Bathymetry and computational mesh of the Saginaw Bay model (a and b). ............... 31
Figure 10: Comparisons between vertically-averaged currents from ADCP observations and
simulations for the year 2009 at the station 216 (a, b) and the station 232 (c, d). ........................ 46
Figure 11: Comparisons between vertically-averaged currents from ADCP observations and
simulations for years 2010 at the station 211 (a,b) and 2011 at the station 8119 (c,d). ............... 47
Figure 12: Comparisons between observed and simulated water level fluctuations at station 5035
(Essexville, Michigan) for the years 2009, 2010 and 2011. ......................................................... 48

x

Figure 13: Comparisons between simulated surface temperatures and observations at the NDBC
buoys from May to September for the years 2009, 2010, and 2011. ............................................ 49
Figure 14: Comparisons between simulated surface temperatures and observations from ADCPs
for the year 2011. .......................................................................................................................... 50
Figure 15: Surface temperature contours for Lake Huron for the summer months July, Aug, and
September in 2009, 2010 and 2011............................................................................................... 53
Figure 16: Mean monthly temperature profile along the transect 3-3 in Figure 3. ...................... 54
Figure 17: Mean monthly temperature profile along the transect 4-4 in Figure 3. ...................... 55
Figure 18: Vertically-averaged currents in Lake Huron averaged over the three year period
(2009 – 2011) for each summer month (July, August and September). ....................................... 57
Figure 19: Vertically-averaged currents in Saginaw Bay averaged over the three year period
(2009 – 2011) for each summer month (July, August and September). ....................................... 58
Figure 20: Mean surface temperature contours in Saginaw Bay for the months of July, August,
and September in 2009, 2010 and 2011 ........................................................................................ 59
Figure 21: Mean vertical velocity profiles at the inlet of Georgian Bay for different layers
confirming the return flow suggested by [Beletsky, et al., 1999]. ................................................ 62
Figure 22: Comparisons between observed positions of GPS-enabled Lagrangian drifters and
simulated particle trajectories based on the Saginaw Bay hydrodynamic model. For each drifter,
the solid lines (observed) and dashed lines (simulated) are shown in the same color. ................. 65
Figure 23: Comparisons between observed and simulated values of single-particle (absolute)
horizontal diffusivities in the Saginaw Bay based on data for different drifters. ......................... 70
Figure 24: Comparisons of the Crossshore and Alongshore Lagrangian Time-Scale. ................ 75
Figure 25: Comparisons of the Crossshore and Alongshore Lagrangian Length-Scale. ............. 75
Figure 26: Schematic of Ice Model. ci is the ice concentration; a is wind stress; w is water
stresss; uw is water current speeds; ui is the ice velocity; ua is wind speed; hi is the ice thickness.
....................................................................................................................................................... 81

xi

Figure 27: Great Lakes Annual Maximum Ice Coverage (1973-2013)
(http://www.glerl.noaa.gov/data/) ................................................................................................. 85
Figure 28: Vertically-Integrated velocity comparison at location HM12 with the model. .......... 86
Figure 29: Vertically-Integrated velocity comparison at location SB32 (above) and SB33
(below) with the model. ................................................................................................................ 87
Figure 30: Time series of observed versus simulated temperature in water column at SB32,
SB33 station .................................................................................................................................. 89
Figure 31: Comparison between lake-wide averaged simulated surface temperatures and
observationals data processed by GLSEA for (a) winter 2010 and (b) winter 2013. (The GLSEA
data extracted from http://coastwatch.glerl.noaa.gov/). ................................................................ 90
Figure 32: Surface temperature contours for Lake Huron for the winter months 2010............... 92
Figure 33: Surface temperature contours for Lake Huron for the winter months 2013............... 93
Figure 34: Surface temperature contours for Saginaw Bay for the winter months 2010............. 94
Figure 35: Surface temperature contours for Saginaw Bay for the winter months 2013............. 95
Figure 36: Lake-wide averaged of simulated water temperature at surface, mid-depth, and
bottom during winter 2010 (above) and 2013 (below). ................................................................ 96
Figure 37: The time series of ice cover in percentage in the winter 2010, 2013 ......................... 99
Figure 38: Map of Lake Huron showing observed and simulated ice cover area for winter 2010.
..................................................................................................................................................... 100
Figure 39: Map of Lake Huron showing observed and simulated ice cover area for winter 2013.
..................................................................................................................................................... 103
Figure 40: Modeled Ice thickness contours for Lake Huron for the winter of 2010. ................ 106
Figure 41: Modeled Ice thickness contours for Lake Huron for the winter of 2013. ................ 107
Figure 42: Mean vertical velocity profiles at the inlet of Georgian Bay for different layers in
winter months.............................................................................................................................. 109
xii

Figure 43: Wind rose plots for winter months over Lake Huron during in the years 2010, 2013.
..................................................................................................................................................... 110
Figure 44: Monthly vertically-averaged currents in Saginaw Bay in winter months 2010 ....... 111
Figure 45: Monthly vertically-averaged currents in Saginaw Bay in winter months 2013 ....... 112
Figure 46: Monthly vertically-averaged currents in Lake Huron in winter months 2010. ........ 113
Figure 47: Monthly vertically-averaged currents in Lake Huron in winter months 2013. ........ 114
Figure 48: Monthly average Ice velocity in Lake Huron in 2010. ............................................. 116
Figure 49: Monthly average Ice velocity in Lake Huron in 2013. ............................................. 117
Figure 50: Monthly average Ice velocity in Saginaw Bay in 2010. ........................................... 118
Figure 51: Monthly average Ice velocity in Saginaw Bay in 2013. ........................................... 119
Figure 52: Time series of water exchange rates for the inner transect (1-1 in Figure 3) and the
outer transect (2-2 in Figure 3) in Saginaw Bay in summer season. .......................................... 128
Figure 53: Wind rose plots for summer months during the years 2009-2011. .......................... 130
Figure 54: Example of a diagram depicting the water exchanges in a structured grid. ............. 131
Figure 55: Map of Saginaw Bay showing the contour of flushing time in hours within the Bay
calculated based on monthly average velocity............................................................................ 133
Figure 56: Time series of water exchange rates for the inner transect (1-1 in Figure 3) and the
outer transect (2-2 in Figure 3) in Saginaw Bay in winter season. ............................................. 135
Figure 57: Contour plots of residence time for Saginaw Bay based on dye concentration
modeling. Figures (a) and (b) are for inner bay dye releases with (a) no river inflow and (b) with
river inflow. Figures (c) and (d) are results based on dye releases for the entire bay with (c) no
river flow and (d) with river flow. (Summer season). ................................................................ 138
Figure 58: Average vertically-integrated concentration in the Saginaw Bay as a function of time.
The red solid lines represent the model results while the blue dashed lines represent an

xiii

exponential decay model fit to the data. The concentration value corresponding to a value of
(1/e) is also marked. (Summer season). ...................................................................................... 139
Figure 59: Contour plots of residence time for Saginaw Bay based on dye concentration
modeling. Figures (a) and (b) are for inner bay dye releases with (a) no river inflow and (b) with
river inflow. Figures (c) and (d) are results based on dye releases for the entire bay with (c) no
river flow and (d) with river flow. (Winter season without Ice Model). .................................... 141
Figure 60: Average vertically-integrated concentration in the Saginaw Bay as a function of time.
The red solid lines represent the model results while the blue dashed lines represent an
exponential decay model fit to the data. The concentration value corresponding to a value of
(1/e) is also marked. (Winter season without Ice Model). .......................................................... 142
Figure 61: Contour plots of residence time for Saginaw Bay based on dye concentration
modeling. Figures (a) and (b) are for inner bay dye releases with (a) no river inflow and (b) with
river inflow. Figures (c) and (d) are results based on dye releases for the entire bay with (c) no
river flow and (d) with river flow. (Winter season with Ice Model). ......................................... 144
Figure 62: Average vertically-integrated concentration in the Saginaw Bay as a function of time.
The red solid lines represent the model results while the blue dashed lines represent an
exponential decay model fit to the data. The concentration value corresponding to a value of
(1/e) is also marked. (Winter season with Ice Model). ............................................................... 145
Figure 63: Mean residence time as a function of river discharge for (a) Inner Bay and (b) Entire
Bay. The blue solid dots represent the model results while the red line represents an exponential
equation fit to the data. The residence times were calculated based on the dye concentration
model under summer meteorological conditions. ....................................................................... 153

xiv

Chapter 1
Introduction
1.1

Problem description

Knowledge of circulation patterns in lakes and oceans is essential to addressing a number of
issues ranging from water quality to ecosystem modeling. Knowledge of lake circulation has
proven to be indispensable in prediction of pollution from point sources or industrial discharge
into large bodies of water. In addition, currents are a primary factor in understanding many
environmental issues such as harmful algal blooms, dead zones, and so on.
Lake Huron is a large aquatic system and the third largest of the Great Lakes by volume.
Lake Huron encompasses three distinct basins (Gerogian Bay, North Channel, and the main
basin which includes Saginaw Bay) as shown in Figure 1, with several islands (St. Joseph Island,
Bois Blanc Island, Drummond Island, Cockburn Island, Clapperton Island, Manitoulin Island,
Fitzwilliam Island, Charity Island, etc.) and bays (Georgian Bay, Saginaw Bay, Thunder Bay).
The prominent bathymetric feature in Lake Huron is the mid-lake ridge that extends
southeastward from Thunder Bay at Alpena on the west side across to Point Clark on the east
side. The ridge (Alpena-Amberley Ridge or AAR) separates the lake into the southwestern and

1

northeastern basins. The shallower southwestern basin has a maximum depth of about 90 m
while the deeper northeastern basin has a maximum depth of about 200 m. In addition to the
complex topography, the wind strongly influences circulation in Lake Huron.
As in the other Great Lakes, Lake Huron has been significantly affected by natural and
anthropogenic activities. The Lake is affected by water quality degradation, invasion of zebra
mussels [Fahnenstiel et al., 1995; Heath et al., 1995; Bierman et al., 2005; Fishman et al.,
2009], toxic contaminants and eutrophication in some localized areas [Nalepa et al., 2007].
Water quality problems such as persistance of Escherichia coli [Alm et al., 2006], and mercury
contamination [Marvin et al., 2004] are also pertinent issues. In addition, recent awareness of the
water level declining in the Lake Michigan-Huron system has made it necessary to study and
compare the circulation and exchanges with the past. Shallow areas and embayments such as
Saginaw Bay are believed to have increased benthic algal growth from invasive mussel [Hecky et
al., 2004]. Lake Huron is a connecting water way between the upper and lower lakes in the Great
Lakes. Therefore, all of these contaminants will impact not only Lake Huron ecosystem but also
the ecosystems of the lower lakes. Knowledge of the water movement and temperature
distribution are of critical importance.
Saginaw Bay is an important part of Lake Huron that receives flow from the Saginaw
River. The Saginaw Bay - Lake Huron system is significantly impacted by human activities and
figures prominently in the list of Areas of Concern [Great Lakes Water Quality Agreement,
1978]. The Saginaw River, which drains an area of over 16,000 km2, is the major source to the
Bay with an average discharge of about 100 m3/s [ Danek and Saylor, 1977]. The Bay is 87 km
long, and the width varies between 20 km and 46 km. The Bay is often divided into the inner
and outer regions based on topographic features and water circulation patterns (Figure 3a). The

2

bay has extensive shallows averaging 4.0m and 12.0m in depth in the inner and outer bays,
respectively.
The Saginaw Bay is of particular interest because it is the region that anthropogenic
activities interact most strongly with the Bay environment. A significant portion of the total
loading from industrial activity is released along the Saginaw River. Persistent contamination of
the water column and bottom sediment has resulted in degraded water quality. It is well-known
that embayed regions tend to accumulate contaminants due to their long residence times [Nixon
et al., 1996]. Offshore transport and mixing with lake waters are the dominant process for
reduction in concentrations of contaminants in embayed regions attached to large lakes [Ge, et
al., 2012; Thupaki, et al., 2010]. While mixing is dominated by small-scale hydrodynamics and
shear [Lawrence et al., 1995; Spydell et al., 2007], advective exchange between Saginaw Bay
and Lake Huron is strongly influenced by large-scale circulation patterns.

3

Figure 1: Map of Lake Huron showing the important physical features: Manitoulin Basin; Bay
City Basin; NC: North Channel; Alpena- Amberley Ridge; Georgian Bay. (source:
http://www.ngdc.noaa.gov/).

4

1.2

Summer Circulation and Thermal Structure

The large scale circulation in Lake Huron was first described by Harrington [1894], who used
data from drift bottle studies to conclude that circulation pattern is generally counterclockwise in
the main basin of Lake Huron. While drift bottles have higher error than modern Lagrangian
drifter designs due to stokes-drift, wind forces, and lack of real-time GPS tracking, many of the
conclusions Harrington arrived at have been supported by later observations. Sloss and Saylor
[1976] used observational datasets from extensive current mooring measurements to confirm that
while the northern two-thirds of the lake is dominated by counterclockwise circulation during
summer, the shallower southern part shows a much more complex circulation pattern. They also
found that stratification during summer allowed the inertial component to dominate open lake
circulation.
Although, Lake Huron and Michigan are considered one ecosystem due to their hydraulic
connectivity, Lake Huron has been poorly studied when compared to Lake Michigan. A number
of studies have used numerical models to examine circulation in Lake Huron. Sheng and Rao
[2006] used a three-dimensional hydrodynamic model based on the primitive equation, z-level
CANDIE ocean model to simulate large-scale circulation and temperature distribution in Lake
Huron. Using a nested-grid methodology, they were able to resolve small-scale hydrodynamics
in the Georgian Bay region. More recently, Anderson and Schwab [2013] used the unstructured
grid Finite Volume Coastal Ocean Model (FVCOM) to develop a Lake Michigan-Huron model
to simulate currents in the Straits of Mackinac. Bai et al., [2013] also implemented the sigmacoordinate, unstructured grid ocean model (GL-FVCOM) in all five Great Lakes with 21 terrain
following layers in the vertical to simulate circulation and thermal structure for 1993-2008.
Physical features such as shallow depths and complex bathymetry in addition to the presence of a
5

sill near the entrance of the Bay City Basin (BCB) produce a complex circulation pattern within
the Saginaw Bay. As a result, hydrodynamic models using coarse grids will not be able to
adequately resolve circulation within the Saginaw Bay.

Shelter Island

Saginaw River

Figure 2: Aerial view of the Saginaw Bay area showing the upstream region along the Saginaw
River, the river mouth, and the Shelter Island (Photo credit : Doc Searls)

Circulation in Lake Huron shows more spatial variation than other large lakes of comparable size
in North America (e.g., Lake Superior and Lake Michigan) due to the presence of several
islands, bays and other complex physical features. Studies using numerical models to examine
circulation in Lake Huron have used multi-resolution (nested) structured grids [Sheng and Rao,
2006] or coarse-resolution unstructured grids [Bai et al., 2013]. In order to reduce the numerical
6

error in simulating the thermocline and horizontal pressure gradient terms, Sheng and Rao [2006]
used the z-level ocean model CANDIE [Sheng et al., 1998]. Accurately representing the physical
features of the lake basin in the computational domain is necessary in order to simulate currents,
thermal structure as well as large-scale hydrodynamics.
Using time-series analyses and numerical simulations, Allender [1975] and Allender and
Green [1976] determined that circulation within the bay is coupled to large-scale circulation in
Lake Huron. Using a combination of Lagrangian drifter observations from current meters
deployed in Saginaw Bay, Danek and Saylor [1977] and Saylor and Danek [1976] determined
the circulation patterns inside Saginaw Bay for different wind directions. They found that when
wind direction was parallel to axis of the bay, exchange rate between the inner and outer bay was
3700m3/s. More recently, in their analysis of the seasonal circulation in the Great Lakes, Beletsky
et al.,[1999] observed that mean currents during winter are higher than during summer months.
Based on the surface current measurements, they also inferred the presence of a return flow out
of Georgian Bay at deeper depths. Schertzer et al.,[2008] have compiled the mean circulation,
inter-lake exchange and climatology for Lake Huron.
Most of the previous studies described the circulation in Lake Huron based on
observational data that is limited in space and time. Description of the annual pattern of
circulation in Lake Huron is impossible because of a lack of sufficient observations [Schertzer et
al.,2008] and overall the system is still poorly understood. Therefore, using a numerical model to
obtain a more detailed knowledge of both spatial and temporal dynamics of water movement in
the lake shows great promise.

7

1.3

Winter Circulation and Ice Cover Extents

The limited literature on the winter circulation of Lake Huron includes ice extents while the
annual cycle of ice formation and loss affects physical processes within the lake and the
ecosystem of the lake [Assel, 2005]. The winter season is characterized by nearly homogeneous
water in Lake Huron as in the other Great Lakes. The circulation during cold weather has
received very little attention not only in Lake Huron, but in all of the Great Lakes as well. The
main reason is due to lack of observational data during the winter season.
One of the first attempts to describe winter circulation in Lake Huron was carried out by
Saylor and Miller [1979]. In their study they used twenty-one current moorings deployed in
Lake Huron during the winter of 1974-1975 for approximately 6 months. They concluded that
there is a strong cyclonic flow throughout the winter with current speeds that ranged from 1 to
7cm/s with a mean of 3cm/s. Because of isothermal conditions in the water column, the winter
currents are barotropic [Beletsky et al., 1999]. It has been known that the Great Lakes have
considerable influence on the regional climate [Bates and Giorgi, 1993; Scott and Huff, 1996].
Rao and Murty [1970] used a homogeneous, vertically integrated numerical model to conclude
that winter circulation in Lake Ontario was driven by wind.
The hydrology of the Great Lakes can affect local meteorological conditions in the long
term. For example, Scott and Huff [1996] estimated that lake-effect precipitation increased 90%
downwind of Lake Huron in the winter time. In addition, understanding the major effects of
colder weather on the lake are crucial in predicting water movement patterns and water
temperature structure. More recently, Fujisaki et al.,[2012] used the Princeton Ocean Model
(POM) which included ice processes to investigate the seasonal variation of ice cover on Lake
Erie. The study showed that the packed ice cover slowed down the surface water velocity. To our
8

knowledge, simulation of lake processes beneath the ice cover in three dimensions has not been
attempted. Therefore, there are large gaps in our understanding of lake circulation, thermal
structure, ice cover extent, and the water exchange rates between Saginaw Bay and the main lake
in the winter.

1.4

Residence and Flushing Times

Residence and flushing times of an estuary are important variables used to estimate the rate of
removal of pollutants from the bay. As concluded by Monsen et al., [2002], first-order transport
time scales such as residence times and flushing times within the bay are useful for
understanding and interpreting the fate and transport of contaminants, nutrient budgets, the
occurrence of harmful algal blooms and for comparing processes and rates across different
ecosystems.
There have been a few studies of the mean residence time in Saginaw Bay and Lake
Huron system. Quinn [1992] estimated the residence time defined as the time it takes to replace
the water volume of Lake Huron as 22 years. Saylor and Danek, [1976] calculated an exchange
rate between the inner bay and outer bay of 3700m3/s and a mean residence time of 26.5 days for
the inner bay. The residence time of any bay may vary with several factors such as wind
direction, and river inflow. Since the Saginaw bay is shallow, the residence time is strongly
influenced by wind. As Saylor and Danek, [1976] pointed out, the residence time of water in the
bay is much longer with the winds perpendicular to the axis of the bay. In an earlier study Dolan
[1975] calculated a residence time of 110 days for the inner bay based on a time-variable
chloride model.
9

Until now, there have been no observational studies on residence time in Saginaw Bay
since the early work cited above. A 3D lake-wide numerical model with additional water quality
models provides a great opportunity to study residence time and flushing time for the bay in
detail for both the summer and winter seasons. This research will estimate the residence time
based on different methods that are presented in Chapter 5.
Significant progress has been made in recent years in characterizing the transport time
scales for surface water bodies using hydrodynamic and transport models [Andutta et al., 2013;
Camacho and Martin, 2013; Phelps et al., 2013; Hsu et al., 2013; Wan et al., 2013; Liu et al.,
2012; Andrado´ttir et al., 2012; Jouon et al., 2006]. Although some estimates of residence times
are available for the Saginaw Bay [Dolan, 1975; Saylor and Danek, 1976], they are nearly 40
years old and much has changed in terms of our ability to observe and model natural systems
since then. Recent concerns associated with environmental issues in the Saginaw Bay highlight
the need for accurate estimates of transport timescales in the bay and one of the objectives of this
work is to fill this gap. Furthermore, identifying the spatial patterns of residence time throughout
Saginaw Bay for the forecasting water quality is critical. Terms such as residence time and
flushing time, however, are defined and used in many different ways in the literature. Following
the suggestions of Monsen et al. [2002] and Bolin and Rodhe [1973] who recommend that the
terms be defined precisely and used with care to avoid confusion, we define these transport time
scales in the next chapter and describe how they are calculated.

10

1.5

Goals of the Thesis

The aim of this study was to characterize summer and winter circulation, thermal
structure and ice cover extent in the Lake Huron - Saginaw Bay system and to quantify residence
and flushing times within the bay. This was performed using a combination of field and satellitebased observations and high-resolution, three-dimensional, unstructured grid numerical model.

1.6

Structure of the Thesis

Chapter 1: Introduction
Chapter 2: Field Observations and Analysis
Chapter 3: Summer Circulation and Thermal Structure
Chapter 4: Ice Cover, Winter Circulation and Thermal Structure
Chapter 5: Residence Time in the Saginaw Bay and Lake Huron System
Chapter 6: Conclusions

11

Figure 3: (a) Map of Lake Huron showing the important physical features. MB: Manitoulin
Basin; BCB: Bay City Basin; NC: North Channel; AAR: Alpena- Amberley Ridge. (b) Map of
Saginaw Bay showing the ADCP locations. SI: Shelter Island. 1-1 and 2-2 in Figure 3(a) are the
inner and outer Saginaw Bay transects respectively and 3-3 and 4-4 are transects for which
vertical temperature profiles are shown.

12

Chapter 2
Field Observations and Analysis
2.1 Observational Data

Observational datasets have a very important role in terms of understanding natural processes
and to validate the numerical model. The hydrodynamic model could not be tested without
adequate data to serve as input and also for calibration purposes. The areas of interest for this
study are Lake Huron and Saginaw Bay. Field observations were undertaken in the vicinity of
Saginaw River mouth and Saginaw Bay mouth as shown in the Figure 3.
The field work was conducted to collect current velocity, and water temperature
throughout the water column. It has been known that the primary factors that affect circulation in
lakes are wind, temperature, bathymetry, and tributary flows. Therefore, in order to obtain a
better understanding of the fundamental hydrodynamics and river plume, several drifters were
released adjacent to the Saginaw River mouth. The observational data outlined above are used to
compare the results from hydrodynamic model and estimate diffusivity using Lagrangian drifters
and particle transport model. Here we present in detail the field data in the next section.

13

2.1.1

Currents
The field work was conducted to collect current velocity data in both summer and winter

seasons. The observational investigations in summer were carried out between July and
September of 2009, 2010, and 2011 using Nortek (www.nortekusa.com) Aquadopp current
profilers (2000 kHz). The ADCPs were set up in a bottom-resting, upward-looking
configuration, at different locations within the Bay and programmed to record current data with a
bin-size of 0.25m and ensemble length of 900s or 1800s to collect an extensive observational
dataset in the near-shore region of Saginaw River mouth (Figure 3). This region is characterized
by very shallow water depths (average depth 2.5m).
There are two observational investigations which contribute to winter datasets for the
year 2010 and 2013. The first field work was conducted in the winter 2010 in which the RD
Instruments ADCPs were deployed in the Saginaw Bay mouth as shown in Figure 3a as SB32
and SB33 stations. The second observational investigation was collected in the winter 2013. A
Teledyne – RD Instruments (http://rdinstruments.com) Workhorse Sentinel ADCP (600 kHz
frequency) was deployed in the area of Hammond Bay with a water depth of about 27m. The
observational dataset provided a time series of water current and temperature in the water
column. Deployment locations and additional details are listed in Table 1.

14

Table 1: Details of ADCP deployments in Lake Huron
Deployment Date

Location

Longitude

Latitude

Depth(m)

07/22/09 - 08/04/09

203

-83.90235

43.673233

3.0

08/04/09-08/18/09

216

-83.897016

43.680066

2.0

08/20/09-09/02/09

232

-83.897016

43.680066

2.9

09/02/09-09/16/09

245

-83.91055

43.699066

3.3

10/14/09-05/26/10

SB32

-83.2609

44.2699

27.0

10/20/09-05/26/10

SB33

-83.2069

44.1434

29.0

07/13/10-07/27/10

194

-83.917683

43.702916

2.4

07/30/10-08/12/10

211

-83.916666

43.702916

2.5

07/27/11-08/28/11

8119

-83.8525

43.688333

3.6

07/27/11-08/28/11

8126

-83.789166

43.670000

4.0

11/09/12-05/19/13

HM12

-84.05447

45.5264

27.0

2.1.2

Temperature

One of the primary factors affecting lake circulation is water temperature in the lake (via the
density gradient term in the momentum equation). As with other large lakes, Lake Huron is
stratified during summer time and has nearly homogeneous temperature during winter time.
Description of the thermal distribution and evolution in Lake Huron is one of the goals of the
study. Thus, the observational investigation which measures water temperature throughout the

15

water column is very useful to validate the results from the hydrodynamic model. In this study,
we used the temperature data measured by GLERL scientists using thermistor chains. Time
series measurements of water temperature were collected at 2 stations (SB32, SB33 as shown in
Figure 3) located along the boundary of outer Saginaw Bay during the summer of 2009 and the
winter of 2010. Details of the measurements at the individual stations are presented below in
which the sensor positions are in meters above the bottom.

Table 2 : Details of Thermistor deployments in Lake Huron
Deployment Date

Location

Longitude

Latitude

The sensor positions (m)

SB32

-83.2609

44.2699

[1.0;5.0;9.0;13.0;17.0;21.0;23.0]

SB33

-83.2069

44.1434

07/16/2009 05/25/2010
07/16/2009 -

[1.0;5.0;9.0;13.0;17.0;21.0;23.0;

05/25/2010

2.1.3

25.0;27.0;29.0]

Lagrangian Drifters

Lagrangian drifters have been widely used in the studies of oceans and large lakes. The
observational data was particularly valuable to estimate the flow field spatially and to validate
the hydrodynamic model. Based on the Lagrangian data the diffusion coefficient can be
estimated providing a better understanding of both fundamental hydrodynamics and
environmental transport issues.
Multiple drifters were deployed between the Saginaw River mouth and Shelter Island (SI)
(Figure 3b) for three days in end of July 2011. Deployment dates and start locations are given in
16

Table 3. The drifters were originally designed by Michael McCormick of the Great Lakes
Environmental Research Laboratory (GLERL), Ann Arbor, Michigan. The drifters used in the
present study were constructed following the original design using a fiberglass and plywood
cross frame with vinyl drogues (area of20′′ × 35′′ ) to reduce latency as shown in Figure 4.
Location of drifters were tracked in real-time using Trackpack ® GPS transmitters programmed
to report their location every 30 min. The transmitters were encased in a water resistant housing
and the manufacturer reported uncertainty in GPS location to be approximately 5m.

17

Table 3: Details of Lagrangian drifter deployments in Saginaw Bay
Drifter ID

Deployment locations

Deployment

Picking up time

time

Duration of
travel

Longitude

Latitude

(hour)

-83.846734

43.655956

07/27/2011 18:44

07/28/2011 18:44

24

-83.842742

43.661277

07/27/2011 18:21

07/28/2011 18:51

23

-83.839867

43.661309

07/27/2011 18:23

07/27/2011 23:53

5

-83.840640

43.661288

07/28/2011 20:21

07/29/2011 18:21

22

-83.839266

43.661374

07/28/2011 20:23

07/29/2011 17:53

21

-83.843365

43.658445

07/29/2011 18:44

08/01/2011 07:44

61

-83.837743

43.662919

07/29/2011 19:23

07/30/2011 03:53

8

Day 1Drifter 1
Day 1Drifter 2
Day 1Drifter 3
Day 2Drifter 1
Day 2Drifter 2
Day 3Drifter 1
Day 3Drifter 2

18

(c)

Figure 4: Photograph showing (a) the Lagrangian drifters, (b) a drifter in water, and (c) the
Trackpack® GPS

19

2.2

Numerical Model

The three-dimensional unstructured grid numerical model (FVCOM) is used to describe
circulation and thermal structure in Lake Huron and Saginaw Bay. Since Lake Huron has
complex topographic features, the use of a variable resolution mesh seems to be a suitable
approach. With the flexibility of unstructured meshes the model is able to represent complex
coastlines to a high degree of accuracy. The model solves the hydrodynamic primitive equations
using the hydrostatic assumption in the vertical direction with the Boussinesq simplification for
convective flows. The continuity, momentum, and temperature equations are shown in Equation
1 to Equation 5.
𝝏𝒖
𝝏𝒙

𝝏𝒖
𝝏𝒕

𝝏𝒚

+

𝝏𝒖

+𝒖

𝝏𝒗

𝝏𝒕

𝝏𝑷
𝝏𝒛

𝝏𝒕

𝝏𝒗

+𝒖

𝝏𝒗

𝝏𝑻

+

𝝏𝒙

𝝏𝒙

𝝏𝒘
𝝏𝒛

=𝟎

+𝒗

𝝏𝒖

+𝒗

𝝏𝒗

𝝏𝒚

𝝏𝒚

(1)

+𝒘

𝝏𝒖

+𝒘

𝝏𝒗

𝝏𝒛

𝝏𝒛

− 𝒇𝒗 = −

+ 𝒇𝒖 = −

𝟏 𝝏𝑷
𝝆𝟎 𝝏𝒙

𝟏 𝝏𝑷
𝝆𝟎 𝝏𝒚

+

+

𝝏
𝝏𝒛

𝝏
𝝏𝒛

= −𝝆𝒈

+𝒖

𝝏𝑻
𝝏𝒙

+𝒗

𝑲𝑴

𝝏𝒖

𝑲𝑴

𝝏𝒗

𝝏𝒛

𝝏𝒛

+ 𝑭𝒖

(2)

+ 𝑭𝒗

(3)

(4)

𝝏𝑻
𝝏𝒚

+𝒘

𝝏𝑻
𝝏𝒛

=

𝝏
𝝏𝒛

𝑲𝑯

𝝏𝑻
𝝏𝒛

+ 𝑭𝑻

(5)

Here x, y, and z are the Cartesian coordinates; 𝑢, 𝑣, 𝑤 are velocity components in horizontal
and vertical directions, respectively; 𝜌 is the density; 𝜌0 is the reference density; 𝑇 is the
20

temperature; 𝑃 is the pressure; 𝑓 is the Coriollis parameter; 𝑔 is the acceleration due to gravity.
The horizontal momentum diffusion terms 𝐹𝑢 , 𝐹𝑣 , have the form:

𝑭𝒖 =

𝑭𝒗 =

𝑭𝑻 =

𝝏
𝝏𝒙

𝝏
𝝏𝒚

𝝏
𝝏𝒙

𝟐𝑨𝒎

𝝏𝒖

𝟐𝑨𝒎

𝝏𝒗

𝑨𝒉

𝝏𝒙

𝝏𝒚

𝝏𝑻

+

𝝏𝒙

+

+

𝝏
𝝏𝒚

𝝏
𝝏𝒚

𝝏
𝝏𝒙

𝑨𝒎

𝝏𝒖

𝑨𝒎

𝝏𝒖

𝑨𝒉

𝝏𝒚

𝝏𝒚

+

𝝏𝒗

+

𝝏𝒗

(6)

𝝏𝒙

(7)

𝝏𝒙

𝝏𝑻

(8)

𝝏𝒚

Vertical eddy viscosity and diffusivity 𝐾𝑀 , 𝐾𝐻 are modeled using the Mellor-Yamada 2.5 level
turbulence closure scheme [Mellor and Yamada, 1982; Galperin et al., 1988]. The horizontal
diffusion coefficients 𝐴𝑚 are calculated using the Smagorinsky turbulence closure model
[Smagorinsky, 1963].

𝑨𝒎 = 𝑪∆𝒙∆𝒚

𝟏

𝝏𝒖 𝟐

𝟐

𝝏𝒙

+

𝟏 𝝏𝒖
𝟐 𝝏𝒚

+

𝝏𝒗 𝟐
𝝏𝒙

+

𝝏𝒗 𝟐
𝝏𝒚

(9)

The horizontal eddy diffusivity 𝐴𝑕 is related to the eddy viscosity 𝐴𝑚 via the turbulent Prandtl
number. 𝑃𝑟 =

𝐴𝑚

𝐴𝑕 = 1.0 was used in this study after carefully testing the effect of this

number on model results. Details of the governing equations and boundary conditions can be
found elsewhere (e.g., Chen et al., [2003]).

21

The surface boundary condition for temperature is given as follows:

𝝏𝑻
𝝏𝒛

=

𝟏

𝑸𝒏 𝒙, 𝒚, 𝒕 − 𝑺𝑾(𝒙, 𝒚, 𝒕)

𝝆𝒄𝒑 𝑲𝒉

(10)

The bottom boundary condition for temperature is given as follows:
𝝏𝑻
𝝏𝒛

=−

𝑨𝑯 𝐭𝐚𝐧 𝛂 𝝏𝑻
𝑲𝒉

(11)

𝝏𝒏

where 𝑄𝑛 𝑥, 𝑦, 𝑡 is the surface net heat flux, 𝑆𝑊(𝑥, 𝑦, 𝑡) is the shortwave radiation flux incident
at the water surface, 𝑐𝑝 is the specific heat of water, 𝐴𝐻 is the horizontal thermal diffusion
coefficient, 𝛼 is the slope of the bottom.
The surface and bottom boundary conditions for velocity u and v are:

𝑲𝒎

𝝏𝒖 𝝏𝒗

,

=

𝑲𝒎

𝝏𝒖 𝝏𝒗

=

𝝏𝒛 𝝏𝒛

,

𝝏𝒛 𝝏𝒛

𝟏
𝝆𝟎

𝟏
𝝆𝟎

𝝉𝒔𝒙 , 𝝉𝒔𝒚

(12)

𝝉𝒃𝒙 , 𝝉𝒃𝒚

(13)

where 𝜏𝑠𝑥 , 𝜏𝑠𝑦 and 𝜏𝑏𝑥 , 𝜏𝑏𝑦 are the x and y components of surface wind and bottom stresses.
The horizontal computational domain consists of non-overlapping unstructured triangular
meshes. The calculated horizontal velocity components are placed at the centroids of triangle.
While all scalar components are calculated at the nodes of triangle. More details about the model
can be found in the study of Chen et al.,[2003].

22

2.3

Uncertainty of measurements

Observational datasets were used in this study for driving and testing the numerical model.
Therefore, the accuracy of data plays a very important role for developing accurate numerical
models. The quality of meteorological data such as air temperature, wind speed, wind direction,
water surface temperature, and dew point is controlled by the National Data Buoy Center
(NDBC) and the National Climatic Data Center (NCDC). Here we present in detail the accuracy
of observational data used in this study in the Table 4.

Table 4: Accuracy of the measurements
No

Parameters

Unit

Accuracy

1

Air Temperature

Degree

+/- 1.0 C

2

Wind Speed

m/s

+/- 1.0

3

Wind Direction

Degree

+/- 10

4

Dew Point

Degree

+/- 1.0 C

5

Water Surface Temperature

Degree

+/- 1.0 C

6

ADCP velocity

cm/s

+/- 0.3% or +/- 0.3 cm/s

7

Drifter Tracks

m

5.0

23

Chapter 3
Summer Circulation and Thermal Structure

3.1

Introduction

Concerns associated with climate change impacts and the increased incidences of anthropogenic
stressors have necessitated an accurate characterization of circulation and thermal structure in
Lake Huron to obtain a better understanding of circulation, thermal structure and residence
times.
Although there have been many studies focusing on circulation in Lake Huron and
Saginaw Bay using observational datasets or numerical models (Sloss and Saylor [1976];
Allender and Green [1976] ; Danek and Saylor [1977]; Beletsky et al., [1999] ; Sheng and Rao
[2006]; Schertzer et al., [2008] ; Bai et al., [2013]) the circulation and temperature distribution in
Lake Huron are still far from being understood. For example, without information about
circulation and results from hydrodynamic models, it is impossible to solve the issues related to
water quality. Furthermore, most of the previous studies described the circulation in Lake Huron
based on observational data that is limited in space and time.
In this chapter, we use a three-dimensional, unstructured grid hydrodynamic model with a
nesting model to examine circulation and thermal structures in the Saginaw Bay - Lake Huron
24

system during the summer months for three consecutive years (2009-2011). The model was
tested against ADCP observations of currents, data from a Lagrangian drifter experiment in the
Saginaw Bay and temperature data from the National Data Buoy Center stations. Using data
from a Lagrangian drifter release conducted near the mouth of Saginaw River, we also calculate
absolute particle diffusion rates within the inner Saginaw Bay near the mouth of the Saginaw
River. The results from the numerical model are then used to calculate volumetric exchange rates
for both the inner and outer Saginaw Bay and to estimate first-order transport time scales such as
the residence time/flushing time and their inter-annual variability. The new estimated flushing
time and residence time will be presented in the Chapter 5.

3.2

Numerical model

3.2.1

Model configuration

The three-dimensional unstructured grid numerical model (FVCOM, Chen et al., [2003]) is used
to describe circulation and thermal structure in Lake Huron and Saginaw Bay. The model solves
the hydrodynamic primitive equations (1-5) presented in the previous chapter.
Two separate hydrodynamic models: (a) the Lake-wide Model (LWM) and (b) Saginaw
Bay Model (SBM) (including the Saginaw River mouth) were setup to resolve large-scale, lakewide circulation and circulation within the Saginaw Bay respectively.
The unstructured meshes for LWM and SBM were created using the Surface Water Modeling
System (SMS, www.aquaveo.com) and bathymetry at node locations was interpolated from the 6
arc-second bathymetry data downloaded from the NOAA National Geophysical Data Center,
Geophysical Data System (GEODAS) website (http://www.ngdc.noaa.gov/mgg/gdas/). In order

25

to obtain the bathymetry at node locations we used a MATLAB (The Mathworks Inc., Natick,
MA, 2013) program based on the TriScatteredInterp function to interpolate from the 6
arc-seconds bathymetry data to node locations using the Natural neighbor interpolation method
(Figure 5). The method allows us to get more accurate interpolated bathymetry data. The LWM
computational domain was treated as a closed boundary and inflow through the Straits of
Mackinac, St.Mary’s River or other tributaries and outflow to Lake St.Clair via the St.Clair River
were not included in the numerical model. The computational meshes used to resolve the
hydrodynamics and evolution of thermal structure in Lake Huron and Saginaw Bay are shown in
Figure 8 and Figure 9b, respectively.

26

Figure 5: Map of Lake Huron showing the bathymetry interpolated from data downloaded from
the NOAA National Geophysical Data Center, Geophysical Data System (GEODAS).

The LWM uses an unstructured grid with 9611 nodes and 17619 triangular elements in
the horizontal (average element size =2.5 km). In general, the LWM’s horizontal model
resolution increases from the center of the main lake, where it is as coarse as 5 km, toward the
coastlines with higher resolution near important features of about 2.0 km. With the horizontal
resolution as shown in the Figure 6, most of the cells are into the range of 2-3 km.

27

Figure 6: A histogram of model grid size of Lake Huron Model (LWM)

Circulation in Saginaw Bay is influenced by large-scale, lake-wide circulation. The SBM
uses an unstructured grid with 8236 nodes and 15575 triangular elements (average element size
approximately 200m) in the horizontal (Figure 7). The minimum element size of 45m is used
near the Saginaw River mouth. Coupling between the LWM and SBM is one-way. Simulated
hourly water levels and velocities at the boundary of the SBM, taken from the LWM, were used
to drive the SBM. Discharge from the Saginaw River, though considered negligible to LWM
circulation and hydrodynamics, was used to provide a flux boundary condition in the SBM.
Daily discharge rates measured at the USGS gaging station on the Saginaw River (#04157000)
(http://waterdata.usgs.gov/) were used to provide the boundary condition at the Saginaw River
mouth in the SBM model.

28

Figure 7: A histogram of model grid size of the Saginaw Bay Model (SBM)

The two models (LWM and SBM) for Lake Huron and Saginaw Bay have 21 vertical sigma
levels. The centers of the sigma levels are located at ([-0.0158; -0.0520; -0.0952; -0.1422; 0.1921; -0.2443; -0.2984; -0.3542; -0.4116; -0.4703; -0.5297; -0.5884; -0.6458; -0.7016; 0.7557; -0.8079; -0.8578; -0.9048; -0.9480; -0.9842]). This model setup will allow accurate
representation of the vertical temperature stratification. Based on the well-known CourantFriedrichs -Levy (CFL) stability criterion, the external and internal time steps used are 4s and
10s for the LWM and 1s and 10s for the SBM, respectively.

29

Figure 8: The computational mesh used for the Lake Huron model (LWM)

30

Figure 9: Bathymetry and computational mesh of the Saginaw Bay model (a and b).

3.2.2

Meteorological Data

The model was driven by the surface winds, heat fluxes, and river discharge. The major
difficulty with meteorological data over Lake Huron is that most of the National Data Buoy
Center Stations (NDBC) are located along the U.S shorelines in the west and the southwest and
very few are located along the Canadian shorelines in the northeast and eastern regions. In this
study we used meteorological data from the NDBC Stations and from the Great Lakes
Observation System (GLOS) (http://data.glos.us/obs/) to drive the hydrodynamic model.
The Hourly standard meteorological data such as wind speed, wind direction, and air temperature
for the model simulation periods in 2009, 2010, and 2011 were downloaded from 14 National
Data Buoy Center Stations (NDBC) (http://www.ndbc.noaa.gov/) around Lake Huron. In
addition, hourly meteorological data from 27 stations were also obtained from the Great Lakes
Observation System (GLOS). The quality of meteorological data is controlled by the National
Data

Buoy

Center

(NDBC)

and

the

National
31

Climatic

Data

Center

(NCDC)

http://gis.ncdc.noaa.gov/map/viewer). The GLOS, NCDC, and NDBC stations used in this study
are shown in Figure 3. Meteorological station locations and parameters in which were measured
at each station are presented in the Table 5, 6, and 7. In the Tables, WSPD is the wind speed,
WDIR is the wind direction, TA is the air temperature, WT is the surface water temperature,
DEWP is the dew point, and CC is the cloud cover. In order to get overwater fields data from
meteorological datasets, three main steps are required: (1) height adjustments to a common 10m
height over the water surface, (2) overland/overwater adjustments, and (3) interpolation. The
methods and procedures are described more thoroughly by Schwab [1999]. The observed
meteorological datasets such as wind speed, air temperature, cloud cover and relative humidity
were interpolated to the computational domain using a smoothened nearest neighbor method
[Schwab, 1999]. Surface heat flux, HFLX, is calculated as

𝑯𝑭𝑳𝑿 = 𝑺𝑾 + 𝑳𝑾 + 𝑯𝒔𝒆𝒏𝒔𝒊𝒃𝒍𝒆 + 𝑯𝒍𝒂𝒕𝒆𝒏𝒕

(14)

The incoming shortwave radiation was calculated as described by Bunker [1976]:

𝑺𝑾 = 𝑺𝑾𝒄𝒍𝒓𝒔𝒌𝒚 𝟏 − 𝒌𝒄

(15)

where 𝑆𝑊𝑐𝑙𝑟𝑠𝑘𝑦 is the clear-sky shortwave radiation, calculated using methods described by
Annear and Wells, [2007]. The cloud cover 𝑐 is converted from NCDC description and 𝑘 is an
empirical parameter. A value of 0.7 was used for k in this study based on [Schwab, 1999]. The
incoming long-wave radiation was calculated following [Parkinson and Washington, 1979]:

𝑳𝑾 = 𝝈𝑻𝟒𝒂 𝟏 − 𝟎. 𝟐𝟔𝟏𝒆

−𝟕.𝟕𝟕×𝟏𝟎−𝟒 𝟐𝟕𝟑−𝑻𝒂

32

𝟐

𝟏 + 𝟎. 𝟐𝟕𝟓𝒄

(16)

Where  is the Stefan-Boltzmann constant (5.67x10-8W m-2 K-4); Ta is the air temperature; c is
cloud cover. Hourly wind and heat flux values were used to force the hydrodynamic models
(LWM and SBM) used in this study. The sensible (𝐻𝑠𝑒𝑛𝑠𝑖𝑏𝑙𝑒 ) and latent (𝐻𝑙𝑎𝑡𝑒𝑛𝑡 ) heat transfers
are calculated at each grid point based on a bulk aerodynamic formulation COARE 2.6
developed by Fairall et al., [1996]. Lake-wide summer circulation for three years (2009, 2010,
and 2011) was simulated using the large-scale hydrodynamic model (LWM). In order to avoid
problems with initialization of the stratified temperature, the model simulations were started on
May 1 with a uniform vertical temperature profile determined from observations at 4 NDBC
buoys (45003, 45008, 45137, and 45143).

33

Table 5: NDBC Meteorological station locations
Station
ID
45003

Location

Latitude

Longitude

NDBC Buoy

45.351

-82.840

45008

NDBC Buoy

44.283

-82.416

Alpena
Harbor
Light, MI
DTLM4
De Tour
Village, MI
FTGM4 Fort Gratiot,
MI
GSLM4
Gravelly
Shoals Light
MI
HRBM4
Harbor
Beach, MI
LTRM4 Little Rapids,
MI
MACM4
Mackinaw
City, MI
PTIM4
Point
Iroquois, MI
RCKM4
Rock Cut,
MI
SBLM4 Saginaw Bay
Light #1, MI
SWPM4
S.W. Pier,
MI
WNEM4
West
Neebish, MI

45.060

-83.424

Period of data
availability
05-06-2009 to 11-30-2009
01-01-2010 to 12-31-2010
04-30-2011 to 12-07-2011
04-05-2012 to 12-31-2012
05-07-2009 to 11-30-2009
01-01-2010 to 12-31-2010
04-30-2011 to 12-14-2011
03-31-2012 to 12-21-2012
01-01-2009 to 06-20-2013

45.992

-83.897

01-01-2009 to 06-20-2013

43.007

-82.423

01-01-2009 to 06-20-2013

44.018

-83.537

01-02-2009 to 06-30-2013

43.847

-82.643

01-01-2009 to 06-30-2013

46.485

-84.300

01-01-2009 to 06-30-2013

45.778

-84.725

03-18-2009 to 06-30-2013

46.485

-84.632

01-01-2009 to 06-30-2013

46.265

-84.192

01-01-2009 to 06-30-2013

43.806

-83.719

01-01-2009 to 06-30-2013

46.502

-84.373

01-01-2009 to 06-30-2013

46.283

-84.205

01-01-2009 to 12-31-2012

APNM4

34

Parameters
WSPD;WDIR;
TA;WT;
WSPD;WDIR;
TA;WT;
WSPD;WDIR;
TA;
WSPD;WDIR;
TA;
WSPD;WDIR;
TA;
WSPD;WDIR;
TA;
WSPD;WDIR;
TA;
WSPD;WDIR;
WSPD;WDIR;
TA;
WSPD;WDIR;
TA;
WSPD;WDIR;
WSPD;WDIR;
TA;
WSPD;WDIR;
TA;
WSPD;WDIR;
TA;

Table 6: GLOS Meteorological station locations
Station

Latitude

Longitude

45137

Station
ID
45137

45.54

-81.01

45143

45143

44.94

-80.63

45149

45149

43.54

-82.07

45151

45151

-79.37

44.50

45152

45152

-79.72

46.23

45154

45154

-82.64

46.05

Sault Ste.
Marie

ANJ

46.48

-84.36

Alpena
County
Regional
Airport
Huron
County
Memorial
Airport
Cheboygan

APN

45.07

-83.56

BAX

43.78

-82.99

CYGM4

45.65

-84.47

HYX

43.43

-83.86

MBS

43.54

-84.08

Saginaw
County
H.W.
Browne
Airport
MBS
International
Airport

35

Period of data
availability
07-15-2009 to 11-31-2009
05-01-2010 to 11-19-2010
05-01-2011 to 11-30-2011
07-15-2009 to 11-31-2009
05-01-2010 to 11-22-2010
05-01-2011 to 11-30-2011
07-15-2009 to 11-31-2009
05-04-2010 to 11-17-2010
07-15-2009 to 10-06-2009
05-17-2010 to 11-02-2010
07-15-2009 to 10-15-2009
05-08-2010 to 10-13-2010
05-15-2011 to 10-28-2011
07-15-2009 to 10-16-2009
05-07-2010 to 11-28-2010
05-12-2011 to 11-14-2011
07-15-2009 to 12-31-2009
01-01-2010 to 11-30-2010
05-01-2011 to 11-30-2011
06-01-2012 to 06-20-2013
07-15-2009 to 12-31-2009
01-01-2010 to 11-30-2010
05-01-2011 to 11-30-2011
06-01-2012 to 06-20-2013
07-15-2009 to 12-31-2009
01-01-2010 to 11-30-2010
05-01-2011 to 11-30-2011
06-01-2012 to 06-20-2013
07-15-2009 to 12-31-2009
01-01-2010 to 11-30-2010
05-01-2011 to 11-30-2011
06-01-2012 to 06-20-2013
07-15-2009 to 12-31-2009
01-01-2010 to 11-30-2010
05-01-2011 to 11-30-2011
06-01-2012 to 06-20-2013
07-15-2009 to 12-31-2009
01-01-2010 to 11-30-2010
05-01-2011 to 11-30-2011
06-01-2012 to 06-20-2013

Parameters
WSPD;WDIR;
TA;WT;
WSPD;WDIR;
TA;WT;
WSPD;WDIR;
TA;WT;
WSPD;WDIR;
TA;WT;
WSPD;WDIR;
TA;WT;
WSPD;WDIR;
TA;WT;
WSPD;WDIR;
TA;DEWP;
WSPD;WDIR;
TA;DEWP;
WSPD;WDIR;
TA;DEWP;
WSPD;WDIR;

WSPD;WDIR;
TA;DEWP;

WSPD;WDIR;
TA;DEWP;

Table 6 (cont’d)
Station

Station
ID
MCD

Latitude

Longitude

45.87

-84.64

OscodaWurtsmith
Airport

OSC

44.45

-83.40

Port Hope

P58

44.02

-82.80

St. Clair
County
International

PHN

42.92

-82.53

Pellston
Regional
Airport of
Emmet
County
Airport
Presque Isle
Light

PLN

45.56

-84.79

PRIM4

45.36

-83.49

Sturgeon
Point

SPTM4

44.71

-83.27

Collingwoo
d Automatic
Weather
Reporting
System
Goderich
Automatic
Weather
Reporting
System

WCO

44.50

-80.22

WGD

43.77

-81.72

Mackinac
Island
Airport

36

Period of data
availability
07-15-2009 to 12-31-2009
01-01-2010 to 11-30-2010
05-01-2011 to 11-30-2011
06-01-2012 to 06-20-2013
07-15-2009 to 12-31-2009
01-01-2010 to 11-30-2010
05-01-2011 to 11-30-2011
06-01-2012 to 06-20-2013
07-15-2009 to 12-31-2009
01-01-2010 to 11-30-2010
05-01-2011 to 11-30-2011
06-01-2012 to 06-20-2013
07-15-2009 to 12-31-2009
01-01-2010 to 11-30-2010
05-01-2011 to 11-30-2011
06-01-2012 to 06-20-2013
07-15-2009 to 12-31-2009
01-01-2010 to 11-30-2010
05-01-2011 to 11-30-2011
06-01-2012 to 06-20-2013

Parameters
WSPD;WDIR;
TA;DEWP;
WSPD;WDIR;
TA;DEWP;
WSPD;WDIR;
TA;DEWP;
WSPD;WDIR;
TA;DEWP;
WSPD;WDIR;
TA;DEWP;

07-15-2009 to 12-31-2009
01-01-2010 to 11-30-2010
05-01-2011 to 11-30-2011
06-01-2012 to 06-20-2013
01-01-2010 to 11-30-2010
05-01-2011 to 11-30-2011
06-01-2012 to 06-20-2013
07-15-2009 to 12-31-2009
01-01-2010 to 11-30-2010
05-01-2011 to 11-30-2011
06-01-2012 to 06-20-2013

WSPD;WDIR;

07-15-2009 to 12-31-2009
01-01-2010 to 11-30-2010
05-01-2011 to 11-30-2011
06-01-2012 to 06-20-2013

WSPD;WDIR;
TA;DEWP;

WSPD;WDIR;
TA;DEWP;CC
WSPD;WDIR;
TA;DEWP;

Table 6 (cont’d)
Station
Ont

Station
ID
YAM

Latitude

Longitude

46.48

-84.50

Elliot Lake
Supplement
ary Aviation
Weather
Reporting
Ont

YEL

46.33

-82.56

YQA

44.97

-79.30

Ont

YVV

44.75

-81.10

Ont

YZE

45.88

-82.57

Sarnia
Airport

YZR

43.00

-82.32

37

Period of data
availability
01-01-2010 to 11-30-2010
05-01-2011 to 11-30-2011
06-01-2012 to 06-20-2013
07-15-2009 to 12-31-2009
01-02-2010 to 11-30-2010
04-01-2011 to 11-30-2011
07-15-2009 to 12-31-2009
01-01-2010 to 11-30-2010
05-01-2011 to 11-30-2011
06-01-2012 to 06-20-2013
07-15-2009 to 12-31-2009
01-01-2010 to 11-30-2010
05-01-2011 to 11-30-2011
06-01-2012 to 06-20-2013
07-15-2009 to 12-31-2009
01-01-2010 to 11-30-2010
05-01-2011 to 11-30-2011
06-01-2012 to 06-20-2013
07-15-2009 to 12-31-2009
01-01-2010 to 11-30-2010
05-01-2011 to 11-30-2011
06-01-2012 to 06-20-2013

Parameters
WSPD;WDIR;
TA;DEWP;CC
WSPD;WDIR;
TA;DEWP;

WSPD;WDIR;
TA;DEWP;
WSPD;WDIR;
TA;DEWP;
WSPD;WDIR;
TA;DEWP;
WSPD;WDIR;
TA;DEWP;

Table 7: NCDC Meteorological station locations
Station
ID
711720
712600
712610
712680
712700
712960
713140
714390
714600
714620
715320
715340
716300
716310
716330
716940
717040

Location
Parry
SoundCanada
Sault Ste
Marie Canada
Goderich Canada
Elliott Lake
- Canada
Colling
Wood Canada
Egbert –
Canada
Barrie-Oro Canada
Cove Island
- Canada
Killarney Canada
Great Duck
Island Canada
Muskoka
Awos Canada
Borden Canada
Muskoka Canada
Mount
Forest Canada
Wiarton
Arpt Canada
Beatrice
Climate Canada
Sarnia Canada

Latitude

Longitude
-80.033

Period of data
availability
01-01-2009 to 06-20-2013

45.333

DEWP;CC

46.483

-84.517

01-01-2009 to 06-20-2013

DEWP;CC

43.767

-81.717

01-01-2009 to 06-20-2013

DEWP;CC

46.350

-82.567

01-01-2009 to 06-20-2013

DEWP;CC

44.500

-80.217

01-01-2009 to 06-20-2013

DEWP;CC

44.233

-79.783

01-01-2009 to 06-20-2013

DEWP;CC

44.483

-79.550

01-01-2009 to 06-20-2013

DEWP;CC

45.333

-81.733

01-01-2009 to 06-20-2013

DEWP;CC

45.967

-81.483

01-01-2009 to 06-20-2013

DEWP;CC

45.633

-82.967

01-01-2009 to 06-20-2013

DEWP;CC

44.967

-79.300

01-01-2009 to 06-20-2013

DEWP;CC

44.250

-79.917

01-01-2009 to 06-20-2013

DEWP;CC

44.967

-79.300

01-01-2009 to 06-20-2013

DEWP;CC

43.983

-80.750

01-01-2009 to 06-20-2013

DEWP;CC

44.750

-81.100

01-01-2009 to 06-20-2013

DEWP;CC

45.133

-79.400

01-01-2009 to 06-20-2013

DEWP;CC

43.000

-82.317

01-01-2009 to 06-20-2013

DEWP;CC

38

Parameters

Table 7 (cont’d)
Station
ID
717300
717310
717330
717334
717460
717670
719560
720321
720375

720629

722007

722125

725384

Location
Sudbury Canada
North Bay Canada
Gore Bay
Airport Canada
Elliot Lake
Muni Canada
Sarnia
Climate Canada
Tobermory
Rcs Canada
Gore Bay
Climate Canada
Cheboygan
Co – United
States
Drummond
Island United
States
Jack
Barstow United
States
Tuscola
Area United
States
Saginaw Co
HW Brown
- United
States
St Clair Co
Intl - United
States

Latitude

Longitude
-80.800

Period of data
availability
01-01-2009 to 06-20-2013

46.617

DEWP;CC

46.350

-79.433

01-01-2009 to 06-20-2013

DEWP;CC

45.883

-82.567

01-01-2009 to 06-20-2013

DEWP;CC

46.350

-82.567

01-01-2009 to 06-20-2013

DEWP;CC

43.000

-82.300

01-01-2009 to 06-20-2013

DEWP;CC

45.233

-081.633

01-01-2009 to 06-20-2013

DEWP;CC

45.883

-82.567

01-01-2009 to 06-20-2013

DEWP;CC

45.654

-84.519

01-01-2009 to 06-20-2013

DEWP;CC

46.007

-83.743

01-01-2009 to 06-20-2013

DEWP;CC

43.663

-84.261

01-01-2009 to 06-20-2013

DEWP;CC

43.459

-83.446

01-01-2009 to 06-20-2013

DEWP;CC

43.433

-83.867

01-01-2009 to 06-20-2013

DEWP;CC

42.911

-82.529

01-01-2009 to 06-20-2013

DEWP;CC

39

Parameters

Table 7 (cont’d)
Station
ID
725386
725406

726379
726390

726395

727340

727344

727347

727417
727435

994063

Location
Port Hope United
States
Huron Co
Mem United
States
MBS Intl United
States
Alpena/Phel
ps Colli United
States
Oscoda
Wurtsmith United
States
Sault Ste
Marie United
States
Chippewa
Co Intl United
States
Pellston
Rgnl Arpt United
States
Presque Isle
Co - United
States
Mackinack
Island United
States
Presque Isle
Light United
States

Latitude

Longitude
-82.793

Period of data
availability
01-01-2009 to 06-20-2013

44.022

DEWP;CC

43.780

-82.986

01-01-2009 to 06-20-2013

DEWP;CC

43.533

-84.080

01-01-2009 to 06-20-2013

DEWP;CC

45.072

-83.564

01-01-2009 to 06-20-2013

DEWP;CC

44.450

-83.40

01-01-2009 to 06-20-2013

DEWP;CC

46.479

-84.357

01-01-2009 to 06-20-2013

DEWP;CC

46.250

-84.467

01-01-2009 to 06-20-2013

DEWP;CC

45.564

-84.793

01-01-2009 to 06-20-2013

DEWP;CC

45.407

-83.813

01-01-2009 to 06-20-2013

DEWP;CC

45.865

-84.637

01-01-2009 to 06-20-2013

DEWP;CC

45.356

-83.492

01-01-2009 to 06-20-2013

DEWP;CC

40

Parameters

Table 7 (cont’d)
Station
ID
994064

994065

994979

995980
997257

997258
997260

997265

997266
997268
997359

Location
Sturgeon Pt
Light GL United
States
Tawas City
Glos Weat United
States
Gravelly
Shoals United
States
Georgian
Bay Canada
De Tour
Village United
States
Ft Gratoit United
States
Harbor
Beach United
States
Point
Iroquois United
States
Rock Cut United
States
SW Pier United
States
Alpena United
States

Latitude

Longitude
-83.273

Period of data
availability
01-01-2009 to 06-20-2013

44.713

DEWP;CC

44.256

-83.443

01-01-2009 to 06-20-2013

DEWP;CC

44.017

-83.533

01-01-2009 to 06-20-2013

DEWP;CC

45.533

-81.017

01-01-2009 to 06-20-2013

DEWP;CC

45.983

-83.900

01-01-2009 to 06-20-2013

DEWP;CC

43.017

-82.417

01-01-2009 to 06-20-2013

DEWP;CC

43.850

-82.633

01-01-2009 to 06-20-2013

DEWP;CC

46.483

-84.633

01-01-2009 to 06-20-2013

DEWP;CC

46.267

-84.183

01-01-2009 to 06-20-2013

DEWP;CC

46.500

-84.367

01-01-2009 to 06-20-2013

DEWP;CC

45.067

-83.417

01-01-2009 to 06-20-2013

DEWP;CC

41

Parameters

Table 7 (cont’d)
Station
ID
997697

998203

998242

998255

998339

Location
Little
Rapids United
States
West
Neebish United
States
Saginaw
Bay Light 1
- United
States
Mackinaw
City United
States
Mouth of
Black River
- United
States

Latitude

Longitude
-84.300

Period of data
availability
01-01-2009 to 06-20-2013

46.483

DEWP;CC

46.283

-84.20

01-01-2009 to 06-20-2013

DEWP;CC

43.800

-83.717

01-01-2009 to 06-20-2013

DEWP;CC

45.783

-84.717

01-01-2009 to 06-20-2013

DEWP;CC

42.975

-82.419

01-01-2009 to 06-20-2013

DEWP;CC

42

Parameters

3.3

Results and Discussions

Lake-wide summer circulation for three years (2009, 2010, and 2011) was simulated using the
large scale hydrodynamic model (LWM). In order to avoid problems with initialization of the
stratified temperature, the model simulations were started in May with zero velocity and uniform
vertical temperature profile determined from observations at 4 NDBC buoys (45003, 45008,
45137, and 45143). As noted by Beletsky and Schwab [2001], the spin-up time for the
hydrodynamic model was very short and the effect of the initial conditions on velocity was
negligible after about a week of simulation time.
In order to quantify model performance and the agreement between data and the model,
we used the root-mean-square error (RMSE) index:

𝑹𝑴𝑺𝑬 =

𝒊=𝒏
𝒊=𝟏

𝒙𝒊 −𝑿𝒊

𝟐

𝒏

(17)

where xi denotes the observed data and Xi is the simulation results; and n is the number of data
points.
3.3.1

Results

3.3.1.1 Currents
Figures 10 and 11 show typical comparisons between the observed and simulated verticallyaveraged velocities (east-west velocity u and north-south velocity v) for all three years based on
the Saginaw Bay Model. The RMSE and standard deviation values for the hydrodynamic model

43

comparisons are summarized in Table 8. The comparisons indicate that model performance was
relatively poor for the year 2009. The main reason is due to the missing data for the GLOS
stations from May to July in the year 2009 while we started the model at the beginning of May in
order to avoid problems with initialization of the stratified temperature. Therefore, the model
used only data from the NDBC stations between May and July. This was not an issue for the
remaining two years (2010 and 2011). In all cases, considering the fact that the ADCPs were
deployed in very shallow areas of the bay (2 to 3 m depth) where the effects from waves and
boat traffic might be significant, we conclude that the model was able to describe the verticallyaveraged velocities reasonably well.

44

Table 8: RMSE and standard deviations between ADCP observed current values and results
from the Saginaw Bay model (SBM)

Deployment Date

RMSE (cm/s)

Standard Deviation

(Eastward,

(Eastward, Northward velocity)

Location
Northward
velocity)
Observed

Modeled

08/04/09-08/18/09

216

2.53, 2.68

1.00,2.00

3.4,4.7

08/20/09-09/02/09

232

3.82,3.9

1.49,2.33

3.51,4.83

07/13/10-07/27/10

194

3.6, 5.7

1.0,0.90

2.78,5.72

07/30/10-08/12/10

211

2.67,3.09

2.07, 2.89

3.8,5.2

07/27/11-08/28/11

8119

3.24, 2.23

2.43,2.36

3.10,2.47

07/27/11-08/28/11

8126

4.6,8.9

3.31,6.36

5.89,4.17

45

Figure 10: Comparisons between vertically-averaged currents from ADCP observations and
simulations for the year 2009 at the station 216 (a, b) and the station 232 (c, d).
46

Figure 11: Comparisons between vertically-averaged currents from ADCP observations and
simulations for years 2010 at the station 211 (a,b) and 2011 at the station 8119 (c,d).
47

In addition to the vertically-averaged currents, the numerical model was able to simulate the
water levels. The comparisons showed very good agreement between the observed and simulated
water level fluctuation at station 5035 (Figure 12).

Figure 12: Comparisons between observed and simulated water level fluctuations at station 5035
(Essexville, Michigan) for the years 2009, 2010 and 2011.
48

3.3.1.2 Temperature
The numerical model was also able to simulate the evolution of surface temperatures accurately
(RMSE values less than 1.7oC), as shown by the comparisons between model results and
observations at NDBC buoys and by ADCP presented in Figure 13 and 14, respectively.

Figure 13: Comparisons between simulated surface temperatures and observations at the NDBC
buoys from May to September for the years 2009, 2010, and 2011.

49

Figure 14: Comparisons between simulated surface temperatures and observations from ADCPs
for the year 2011.

While comparisons with observations made at a point can be used to quantify model
accuracy, dynamic lake-wide processes can be understood by examining the spatial variability in
temperature and currents are more clearly shown in the contour plots presented in Figure 15, and
18. Vertically-integrated currents based on the original unstructured grid, averaged over each
summer month for the three year simulation period, were interpolated to a uniform Cartesian grid
with a step size of 5 km and shown as vector plots in Figure 18. The dominant current patterns
are marked using red lines with arrows showing the direction of flow. The lake-wide circulation
became progressively more pronounced and complex with more gyres appearing as time
progressed from July to September while circulation in the Saginaw Bay was strong during the
month of August in which the direction of maximum wind velocity (Figure 53) is along the axis
of the bay. Details of summer circulation in the Saginaw Bay are shown in Figure 19 by
interpolating the vertically-integrated currents from the unstructured grid (averaged over the 3
50

year period) to a uniform Cartesian grid with a resolution of 2 km. We notice the presence of an
anticyclonic (that is, clockwise-rotating) gyre at the mouth of the outer Saginaw Bay for all three
summer months. As noted by others [Sheng and Rao, 2006; Beletsky et al., 1999], mean
circulation in Lake Huron is predominantly cyclonic (counterclockwise) during the summer
months (Figure 18). Generally, current speeds are higher closer to the shore and dominated by
the alongshore component. Currents in the offshore regions of the lake are smaller and organized
into gyres strongly influenced by topographical features (e.g., the Manitoulin Basin (MB) and the
Alpena-Amberley Ridge (AAR) in Figure 1). As expected, vertically-averaged currents in the
shallower parts of the lake (south and south-west) are higher than values in the deeper parts of
the lake (north and north-east). Exchange between Lake Huron and the two main bays (Saginaw
Bay and Georgian Bay) is complex and influenced by a number of factors, including topography,
stratification as well as orientation of the shoreline. For example, Figure 21 shows that mean
surface flow is into the Georgian Bay which is associated with a return flow in the deeper layers
out of Georgian Bay.
As in the other Great Lakes, Lake Huron shows a significant seasonal variation in the
thermal structure. The contour plots of surface temperature presented in Figure15 show the
gradual warming of the surface layers as the summer progresses. As expected, the shallower
parts of the lake warm more quickly and reach the peak of about 25oC in August. However, the
shallower Saginaw Bay reaches this value early by July. The dramatic inter-annual variability in
surface temperature is also clear from the results presented in Figure 15, which show the
temperature contours for all three summer months and for all three years considered (20092011). The vertical variability in temperature is shown more clearly by the temperature profile of
the water column along the transects 3-3 and 4-4 in Figures 16 and 17, respectively. The thermal

51

structure of Lake Huron is affected by the central Alpena-Amberley Ridge (AAR) in the main
lake as well as the presence of the numerous islands. The vertical thermal structure is similar and
comparable to previous observations [Saylor and Miller, 1979; Boyce et al., 1989], and results
from the numerical simulations in the Great Lakes [Beletsky and Schwab, 2001; Sheng and Rao,
2006]. We can also see the presence of Ekman flow-driven coastal upwelling in the northeast of
the lake [Plattner et al., 2006] as well as topographically driven upwelling regions associated
with sharp changes in bathymetry (e.g. AAR). The impact of these upwelling regions can also be
observed in Figure 15 in the form of significantly lower surface water temperatures located
north-east of the AAR. These simulated features are also observed in the GLSEA (Great
Lakes Surface Environmental Analysis) daily composite temperature images produced from the
thermal AVHRR channels of NOAA’s polar orbiting weather satellites [Plattner et al., 2006;
Leshkevich et al., 1992].

52

Figure 15: Surface temperature contours for Lake Huron for the summer months July, Aug, and
September in 2009, 2010 and 2011.
53

Figure 16: Mean monthly temperature profile along the transect 3-3 in Figure 3.
54

Figure 17: Mean monthly temperature profile along the transect 4-4 in Figure 3.

55

3.3.2

Discussion

Results from the high-resolution Saginaw Bay hydrodynamic model (SBM) are presented in
Figures 10 and 11. Summary statistics for the comparisons between the ADCP observations and
model results are provided as RMSE values in Table 8.
Large-scale circulation in Lake Huron is influenced its key physical features (several large
islands, exchange with Georgian Bay and Saginaw Bay, and the mid-lake AAR ridge).
Field observations and results from numerical models have shown that circulation in the main
lake is generally counterclockwise (cyclonic) but not as organized as in Lake Superior due to the
presence of the mid-lake ridge. The lake-wide numerical model was able to simulate this feature
during the summer months. The mean current speed was around 3 cm/s while currents as high as
20cm/s were found in Eastern Georgian Bay. Exchange with Georgian Bay was uniform over the
summer months and the mean surface current speed in the channel connecting Georgian Bay
with Lake Huron was about 4cm/s. Surface currents into Georgian Bay during the summer
months are associated with a return flow at greater depths as indicated by Beletsky et al., [1999].
This feature was confirmed by our numerical model as shown by the vertical profile of mean
current velocity at the mouth of Georgian Bay presented in Figure 21. The mean velocity in the
deeper layers is about 1 cm/s and generally flows from Georgian Bay into Lake Huron. Here the
surface layer coincides with the epilimnion and deeper layer coincides with the hypolimnion.
Therefore, it is likely that the direction of flow at a particular depth is influenced by the strength
of the stratification and changes from year to year. Inter-annual variability of currents was
evaluated using mean current speeds for different months. We found that current speeds were
highest during the year 2010, with a mean current speed of about 4 cm/s.

56

Figure 18: Vertically-averaged currents in Lake Huron averaged over the three year period
(2009 – 2011) for each summer month (July, August and September).
57

Figure 19: Vertically-averaged currents in Saginaw Bay averaged over the three year period
(2009 – 2011) for each summer month (July, August and September).

58

Figure 20: Mean surface temperature contours in Saginaw Bay for the months of July, August,
and September in 2009, 2010 and 2011

59

As a result of differences in area, volume and size of entrance channel, the circulation
patterns within Saginaw Bay and Georgian Bay are very different. Mean summer circulation in
Georgian Bay was largely counterclockwise with a mean speed of about 3cm/s. However,
circulation in Saginaw Bay was much more complex and depended on the wind direction as
noted by Saylor and Danek [1976] and Danek and Saylor [1977]. The circulation at the entrance
of Saginaw Bay was either clockwise (for southwest wind direction) or counterclockwise (for
northeast wind direction) depending on the wind direction. Winds blowing perpendicular to the
axis of the bay did not significantly affect currents within the bay. Current speeds within
Saginaw Bay are generally higher than in the lake resulting in higher mixing due to dispersion.
Another important feature of circulation within Saginaw Bay was the difference between
hydrodynamics in the inner and outer bays and the effect on exchange between the two. We
found that currents along the coast are aligned with the prevailing wind direction. Currents in the
central (deeper) part of the bay flow in the opposite direction.
Similar to other large, deep lakes, the thermal structure in Lake Huron shows significant
seasonal variability. The hydrodynamic model simulations were initialized with constant
temperature based on the observed water temperature measured at the buoys. As shown by the
vertical temperature profiles presented in Figures16 and 17, the process of stratification is
generally completed by August. The thermocline was located at a depth of about 9m with an
inter-annual variability of less than 4 m. By including a surface wind-wave mixing scheme, Bai
et al., [2013] reproduced a reasonable thermal structure with a sharp thermocline located at a
depth of about 15 m in the Great Lakes in August. Since the hydrodynamic models use the
σ−coordinate, we expect a marginally diffused thermocline as noted in [Beletsky and Schwab,
2001]. Due to the lack of any temperature profile measurements in the vertical during our

60

simulating period, a quantitative comparison of the accuracy of the thermal model is not
presented; however, comparisons of temperatures in deeper parts of the lake measured by
ADCPs are shown in Figure 14 which show a good agreement between observations and the
model (RMSE values are less than 0.6oC). As shown in Figure 15, the southern (shallower) part
of the lake warms more quickly. However, as shown in Figure 13, surface water temperature in
the lake has a high inter-annual variability especially in the beginning of summer (June). These
variations in monthly mean water temperatures could be due to global climate dynamics. The
influence of the mid-lake ridge on thermal structure in the main lake is also clear from Figure 16.
The deepening thermocline reaches the mid-lake ridge (~20 m) by around August. At this point
exchange between the hypolimnion of the northeastern and southwestern basins is limited and
exchange between the two basins is largely in the surface (epilimnion) layers.
Several researchers have observed topography-driven upwelling in the marine
environment [e.g., Figueroa and Moffat, 2000]. However, such features have not been observed
in the Great Lakes to the best of our knowledge. Given the highly variable topography in Lake
Huron (compared to the other Laurentian Great Lakes), conditions are favorable for the
upwelling and hydraulic jumps at several locations the main lake as well as near Saginaw Bay.
Coastal upwelling possibly due to Ekman flow driven by the prevailing wind is observed
near the northeastern part of the lake. Figure 16 shows the temperature profile along the
transect 3-3 marked in Figure 3a. During the stratified summer period, we can see that
colder

waters

are forced upward near the AAR. This indicates the presence of a

topographically-driven upwelling and possible ventilated thermocline in central Lake Huron.
Another region of interest for this study is the sill at the mouth of Saginaw Bay (BCB in Figure
3). Since the depth here is about15m and the upper-mixed layer reaches the lake-bottom, no

61

upwelling regions are evident here. However, it is possible that strong currents as a result of
winds from the southwest could setup conditions favorable for a hydraulic jump.

Figure 21: Mean vertical velocity profiles at the inlet of Georgian Bay for different layers
confirming the return flow suggested by [Beletsky, et al., 1999].

62

3.4.3

Particle Tracking Model

The high-resolution Saginaw Bay hydrodynamic model was able to resolve the smallscale bathymetry features and their influence on the circulation within the bay. We used these
results to simulate the trajectory of Lagrangian drifters that were released near the mouth of the
Saginaw River.
The Lagrangian particle tracking consists of solving a nonlinear system of ordinary
differential equations as follows (Chen et al., [2003]):

𝒅𝒙
𝒅𝒕

= 𝑼 𝒙 𝒕 ,𝒕

(18)

where x is the particle position at a time t, dx/dt is the rate of change of the particle position in
time and U(x,t) is the 3-D velocity field. FVCOM model used the explicit Runge-Kutta multi
step methods to solve the equations in which the particle position at a time t is calculated from
solving the discrete integral:

𝒙 𝒕 = 𝒙 𝒕𝒏 +

𝒕
𝑼
𝒕𝒏

𝒙 𝒕 , 𝝉 𝒅𝝉

(19)

Assume that x(tn) = xn is the position of a particle at time t=tn,
x(tn+1) = xn+1 is the new position of a particle at time tn+1= tn+∆t and can be determined by the 4th
order explicit Runge-Kutta method as follows:

𝝐𝟏 = 𝒙𝒏

(20)

63

𝟏

𝝐𝟐 = 𝒙𝒏 + ∆𝒕. 𝑼 ∈𝟏

(21)

𝟐

𝟏

𝝐𝟑 = 𝒙𝒏 + ∆𝒕. 𝑼 ∈𝟐

(22)

𝝐𝟒 = 𝒙𝒏 + ∆𝒕. 𝑼 ∈𝟑

(23)

𝟐

𝒙𝒏+𝟏 = 𝒙𝒏 + ∆𝒕

𝑼(∈𝟏 )
𝟔

+

𝑼(∈𝟐 )
𝟑

+

𝑼(∈𝟑 )
𝟑

+

𝑼(∈𝟒 )
𝟔

(24)

where ∆t is the time step. Details of the particle tracking model can be found elsewhere (e.g.,
Chen et al., [2003]).
Currents from SBM model are used to simulate particles through the Bay. In order to
reduce the effects from initial conditions of current we used the “restart mode” in the model. In
this mode, the particles were released with the initial conditions based on the current velocity
calculated from the hydrodynamic model. This technique allows the particles to move
immediately after being released.
The results from the Lagrangian drifter study and comparisons with observed drifter
tracks are presented in Figure 22. As shown by the starting points (initial positions of the drifters
at time t = 0) marked in the figure, the drifters were released at locations close to the mouth of
the Saginaw River where the currents are complex due to plume dynamics and the presence of
Shelter Island (SI in Figure 3b). As shown in Figure 22, the Lagrangian drifters moved in either
an eastward or northward direction after release near the Saginaw River mouth. The general
direction of the drifters was affected by the prevailing winds and indicative of the complex
plume dynamics near the mouth of a large river. The model simulated particle trajectories
64

compare well with GPS locations of the Lagrangian drifters on all the days. It must be noted that
during the field study, it was observed that the drifters would often get in the shallower depths
and get entangled in vegetation present close to the shore in this part of the bay. As a result,
mean buoy speeds were used to determine the status of the buoy and comparisons with model are
presented only till the first grounding event. Long-term drifter trajectories might track the
gyres in the inner and outer Saginaw Bay that are predicted by the numerical model.

Figure 22: Comparisons between observed positions of GPS-enabled Lagrangian drifters and
simulated particle trajectories based on the Saginaw Bay hydrodynamic model. For each drifter,
the solid lines (observed) and dashed lines (simulated) are shown in the same color.

65

3.4.4

Dispersion coefficient in Saginaw Bay

Studies of dye (tracer) diffusion in the open lake and ocean have found that while mixing rates in
open waters follow the Okubo 4/3 power-law scaling [Okubo, 1971], mixing rates close to shore
are higher due to shear effects resulting from lateral and bottom boundaries [Lawrence et al.,
1995; Peeters et al., 1996; Stoket and Imberger, 2003; Ojo, et al., 2006]. Based on data from a
continuous Rhodamine-WT dye release study conducted during summer 2008, [Thupaki et al.,
2013] estimated a horizontal dispersion coefficient of about 5.6 m2/s close to the shore in
southern Lake Michigan. Dispersion rates in the near-shore can also be calculated using
Lagrangian drifters, which have the additional advantage of higher sampling rate [Spydell et al.,
2007; Nekouee, 2010]. Spydell et al., [2007] calculated the absolute (one-particle) and relative
(two-particle) dispersion rates using statistical analyses of drifter observations. They found that
the relative diffusion rate 𝐷 for Lagrangian particles followed the scaling 𝐷~𝑙2/3 in the surfzone, where 𝑙 is the particle separation distance. Brown et al., [2009] found that rates are higher
in the presence of rip currents with relative diffusion rates following the relation 𝐷~𝑙1/5 on
rip-channel beaches. In the presence of sheared along-shore currents Spydell et al., [2009] found
that along-shore diffusivity is greater than cross-shore diffusivity, but still lower than the
diffusion rate predicted for 2D inertial sub-range 𝐷~𝑙4/3 . Absolute diffusion rates can also be
calculated from single-point Lagrangian drifter statistics by assuming isotropic turbulence. Using
particle diffusion trajectories in the surf-zone, Spydell et al., [2007] found that diffusion rates in
the along-shore and cross-shore directions varied from 4.5 m2/s to 0.7 m2/s over a two-day
period.

66

Data from the Lagrangian drifter study were also used to calculate the absolute dispersion
rates near the mouth of Saginaw River. The recorded GPS locations based on the drifter positions
were used to calculate the one-point dispersion statistics in Saginaw Bay near the mouth of the
Saginaw River. Using the hydrodynamic results from the SBM, a Lagrangian particle-tracking
model was used to simulate the drifter trajectories. The simplest case for dispersion in turbulent
flow involves the statistical analysis of particles released from a point source in isotropic
turbulent flow. Lagrangian statistics provide useful estimates of time and space scales over
which flow fields are correlated and are based on the velocity autocorrelation function:

𝝆 𝝉 =

𝒖𝒊 𝒕 𝒖𝒊 (𝒕+𝝉)

(25)

𝒖𝒊 (𝒕)𝟐

where 𝜌(𝜏) is the velocity autocorrelation function as a function of the time lag 𝜏, 𝑢𝑖 for i=1,2
denote the fluctuations in the along-shore and cross-shore components of the velocities (𝑈𝑖 , i =
1,2) respectively defined as 𝑢𝑖 = 𝑈𝑖 − 𝑈𝑖 , where 𝑈 denotes the time-average of the velocity
over the period T, where T is the total period of integration. Taylor [1921] established that oneparticle diffusivity in x-direction 𝜅𝑥𝑥 is related to the derivative of the variance in displacement
[Pope, 2000]:

𝒌𝒙𝒙 =

𝟏 𝒅𝝈𝟐𝒙𝒙 𝒕
𝟐

𝝈𝟐𝒙𝒙 = 𝟐𝒖′𝟐

where 𝑢′ = 𝑈12

1/2

(26)

𝒅𝒕

𝒕
𝟎

𝒕 − 𝝉 𝝆 𝝉 𝒅𝝉

(27)

is the root mean squared velocity in the x-direction and t and 𝜏 denote the

time and time-lag respectively. 𝜎𝑥𝑥 is the standard deviation of particle displacements and can be
67

calculated by integrating the Lagrangian auto-correlation (auto-covariance) function 𝜌(𝜏). The
equation to calculate one-particle diffusivity in the y-direction 𝑘𝑦𝑦 is similar to Equation 26.
Diffusivity values can be calculated using the zonal and meridional components of the velocity.
However, since the Saginaw Bay shoreline is almost perpendicular to the axis of the bay near the
mouth of the Saginaw River where the drifters were released, we used the along-shore and crossshore velocities to compute the diffusivities in those directions. We used a MATLAB (The
Mathworks Inc., Natick, MA, 2013) program based on the autocorr function to calculate
diffusivity values in the alongshore (x) and cross-shore (y) directions for observed and simulated
drifter tracks.
Quality controlled data from the Lagrangian drifters were used to calculate the observed
dispersion rates in Saginaw Bay near the Saginaw River plume. Results from the drifter model
were used to calculate the model predicted dispersion coefficients and are presented in Table 9.
Comparisons between the observed and simulated alongshore and cross shore dispersion
coefficients are shown in Figure 23 for different drifters. Combined wind and wave processes
contribute to the observed drifter tracks and dispersion in the field, therefore some differences
can be expected with comparisons based on a purely hydrodynamic model. The comparisons in
Figure 23 show that the observed and simulated dispersion coefficients have similar magnitude
and trends. For example, equation (26) predicts a linear spreading rate for short times and
square-root spreading for the other extreme of very large times. Observed and simulated
spreading rates are generally in agreement in describing these limits (Figure 23).
Figure 23 shows a comparison of the model derived Lagrangian single-point (absolute)
dispersion rates in the along shore and cross-shore directions (kxx, kyy) with observed values. It is
clear from the observed values that dispersion is (a) highly anisotropic and (b) variable with
68

time. The model is able to predict some of this variability. However, with time the difference
between observed and simulated values increases.

Table 9: Maximum values of absolute dispersion rates (in m2/s) for Lagrangian drifters and
simulated particle paths in Saginaw Bay
Cross-shore diffusivity

Along shore diffusivity

Day
Observed

Modeled

Observed

Modeled

Day 1-Drifter 2

1.33

7.0

3.9

4.1

Day 2-Drifter 2

0.16

0.44

2.87

3.1

Day 3-Drifter 2

0.66

0.29

12.5

5.7

Closer to the mouth of Saginaw River, we also notice the presence of two additional, smaller
recirculating cells behind the Shelter Island (not shown in the figures). These flow features, not
resolved by the lake-wide model, play an important role in controlling the transport of
contaminants in the Saginaw River plume.

69

Figure 23: Comparisons between observed and simulated values of single-particle (absolute)
horizontal diffusivities in the Saginaw Bay based on data for different drifters.

70

3.4.5

Lagrangian time-scale and length-scale

The Lagrangian integral time and length scales are useful for estimating the decorrelation scales
of the surface water flow. The mean speed (𝑆), mean alongshore velocity (𝑈), and mean
crossshore velocity (𝑉) were calculated according to McCormick et al., [2002/2008]:

𝑺=
𝑼=
𝑽=

𝟏
𝒏
𝟏
𝒏

𝟏
𝒏

𝒏
𝒊=𝟏

𝒖𝟐𝒊 + 𝒗𝟐𝒊

𝒏
𝒊=𝟏 𝒖𝒊

(28)

(29)

𝒏
𝒊=𝟏 𝒗𝒊

(30)

where 𝑢𝑖 , 𝑣𝑖 are the East and North components of the velocity for time index i, and n is the total
number of time segments for each drifter. The velocity autocorrelation function (equation 25)
forms the basis for the computation of the integral length and time scales. The function is
evaluated separately in each direction (x and y) to compute the length and time scales in both
directions.
The Lagrangian integral time-scale (𝑳𝑻𝒊𝒎𝒆 ) is computed as:

𝑳𝑻𝒊𝒎𝒆 =

𝑻
𝝆
𝟎

𝝉 𝒅𝝉

(31)

where the upper bound of integration (T) is defined as the first zero crossing of the
autocorrelation function 𝜌(𝜏). The Lagrangian integral time-scales are associated with integral
length-scales (𝑳𝒍𝒆𝒏𝒈𝒕𝒉) as follows:

71

𝑳𝒍𝒆𝒏𝒈𝒕𝒉 = 𝒖′𝟐𝒊

𝟏

𝟐

𝑻
𝝆
𝟎

𝝉 𝒅𝝉

(32)

where 𝑢𝑖′ (i=1,2) denote the fluctuations in the along-shore and cross-shore components of the
velocities defined as 𝑢𝑖′ = 𝑢𝑖 − 𝑢𝑖 , where 𝑢𝑖 denotes the time-average of the velocity over the
period T.
Details of the drifter deployments in Saginaw Bay including drifters ID, deployment locations
and duration of travel for each drifter are shown in the Table 3. Here we present the Lagrangian
statistics for drifters in the Table 10.

72

Table 10: Lagrangian Statistics for Drifters in Saginaw Bay
Drifter

Mean

Mean

Mean

Lagrangian

Lagrangian

ID

Speed

Alongshore

Cross-

integral time-scale

integral length-

(𝑺)

velocity

shore

(hour)

scale

(cm/s)

(𝑼)

velocity

(cm/s)

(𝑽)

Along-

Cross-

Along-

Cross-

(cm/s)

shore

shore

shore

shore

(m)

Day1-

3.15

-1.4

1.88

2.02

1.402

115.48

123.8

Drifter1

2.54

0

2.02

2.68

0.92

144.15

44.08

Day1-

3.07

-1.69

1.76

3.38

2.57

138.14

194.99

Drifter2

2.83

-0.15

1.82

2.52

0.59

188.59

39.30

Day2-

2.68

2.48

0.3

0.927

1.365

31.11

47.39

Drifter1

2.49

2.16

-0.24

1.47

0.557

110.99

18.70

Day2-

2.90

2.69

0.12

1.02

1.15

36.39

42.51

Drifter2

3.38

3.23

-0.11

0.448

0.538

41.19

17.58

Day3-

4.81

4.15

1.14

2.37

5.80

442.13

351.79

Drifter1

4.04

3.77

0.81

0.065

0.63

16.33

46.76

Day3-

4.00

3.90

-0.39

0.958

0.746

38.19

25.82

Drifter2

4.64

4.46

-0.32

1.088

0.575

136.53

36.35

MEAN

3.44

1.69

0.8

1.77

2.17

133.5

131.1

3.32

2.24

0.66

1.38

0.64

106.3

33.8

(Note: The simulation and observation data are presented in black and red color respectively).

73

The Lagrangian statistics for Saginaw Bay are shown in the Table 10 and Figures 24, 25.
Examining Table 10, the overall mean observed and simulated speed for all drifters are 3.44 cm/s
and 3.32 cm/s, respectively where the means for the velocity components are 1.69 cm/s (alongshore velocity) and 0.8 cm/s (cross-shore velocity) obtained from simulation. The corresponding
numbers are 2.24 cm/s and 0.66 from observation data.

In theory, the Lagrangian integral time and length scales will describe correctly the behavior of a
water parcel if the drifter perfectly tags the movement of the water parcel (that is, there is no
slippage due to wind and wave effects). It should be noted there are no drifters can meet this
condition. The decorrelation time calculated from observation datasets ranged from 0.06h to
2.68h in the alongshore direction, and from 0.53h to 0.92h in the cross-shore direction. The
corresponding numbers calculated from simulation data are 0.92h to 3.38h in the alongshore
direction, and from 0.74h to 5.8h in the cross-shore direction. The Lagrangian integral lengthscales were calculated following equation 32. The results of the Lagrangian length-scales are
shown in the Table 10 and Figure 25. The results show that the integral length scales based on
observations are anisotropic. The anisotropy of the Lagrangian length scales are also observed by
McCormick et al.,[2002] for Lake Michigan. In comparison, the length scales based on
simulation are nearly isotropic. The difference can be attributed to the lack of a wave model that
accounts for the combined effects of wind and wave processes on the observed drifter tracks.
Future modeling efforts including wave processes may be needed in order to improve the
simulation of the Lagrangian drifters.

74

Crossshore Lagrangian Time-Scale (hour)

7
6
5
4

Modeled
Observed

3
2
1
0
0

1

2

3

4

Alongshore Lagrangian Time-Scale (hour)

Crossshore Lagrangian Length-Scale (m)

Figure 24: Comparisons of the Crossshore and Alongshore Lagrangian Time-Scale.

400
350
300
250

Modeled
Observed

200
150
100
50
0
0

100

200

300

400

500

Alongshore Lagrangian Length-Scale (m)

Figure 25: Comparisons of the Crossshore and Alongshore Lagrangian Length-Scale.

75

3.4

Conclusions

Summer circulation in Lake Huron and Saginaw bay was examined using an unstructured
grid hydrodynamic model. The model was tested using ADCP observations in Saginaw Bay,
surface temperature measurements at NDBC buoys, and Lagrangian drifter trajectories. Using
the unstructured grid model we were able to simulate the small-scale circulation features in
Saginaw Bay and Georgian Bay as well as the lake-wide circulation in Lake Huron. The
hydrodynamic model was able to reproduce the anti-cyclonic (counterclockwise) lake-wide
circulation patterns that have been observed in earlier numerical and field studies on Lake
Huron. We found that mean current speeds close to the shore can be as high as 20 cm/s
compared with the lake-wide average of 3 cm/s. The complex topography coupled with
stratification and rapid warming of the southern part of the lake results in high inter-annual
variability in surface temperature values. Several potential upwelling regions have been
identified where topographically-driven upwelling brings cold water from the deeper layers to
the surface. One such region is in the strait connecting Lake Huron with Georgian Bay. In
addition to upwelling, our simulations show that in the surface layer (epilimnion) mean flow
during summer is from Lake Huron into Georgian Bay, which sets up a return flow in the deeper
layers (hypolimnion) from the bay into the lake. Since the hypolimnion is 3-4 times thicker than
the epilimnion in this part of the lake, surface velocities are several times larger in magnitude
than the velocities in the deeper layers.
Results from the Lagrangian drifter study in Saginaw Bay showed the highly anisotropic and
time-variable nature of dispersion near the mouth of Saginaw River. The numerical model was
able to predict the drifter trajectories and dispersion coefficients reasonably well. However,
differences between observed and predicted values increased with time.

76

Chapter 4
Ice Cover, Winter Circulation and Thermal
Structure
4.1

Introduction

When simulating hydrodynamics in lakes over the seasonal or annual cycle, most studies avoid
simulating for winter time due to the lack of observational datasets. The limitations of
instruments (e.g., the ability to measure currents accurately when the water surface is covered
with ice) plus adverse weather conditions are the primary reasons for the lack of many studies on
winter circulation in the Great Lakes in the past. As a result, most of the evaluation about
environmental consequences related to human activities was based on the hydrodynamic results
for ice-free conditions. Meanwhile, the hydrodynamic conditions with ice cover are different
from those of an ice-free condition. The ice cover significantly influences lake hydrodynamics
and water temperature distribution by inhibiting wind stress and heat flux at the surface. Besides
their importance in understanding lake hydrodynamics under winter conditions, a knowledge of
and the ability to accurately predict ice cover, ice thickness, ice concentration, and the duration
of ice in Lake Huron is very useful for icebreaking operations and for managing the lake
77

ecosystem during winter periods. The winter season is characterized by nearly homogeneous
water in Lake Huron as in the other Great Lakes. In a normal year, the Great Lakes have an ice
cover partially from December to April.
One of the first attempts to describe winter circulation in Lake Huron was carried out by
Saylor and Miller [1979]. In the study they used twenty-one current moorings deployed in Lake
Huron during the winter of 1974-1975 for approximately 6 months. They concluded that there is
a strong cyclonic flow persisting throughout the winter with current speeds that ranged from 1 to
7cm/s and a mean of 3cm/s. The winter currents are barotropic because of isothermal conditions
in the water column [Beletsky et al., 1999]. It has been known that the Great Lakes have
considerable influence on the regional climate [Bates and Giorgi, 1993; Scott and Huff, 1996].
In addition, ice cover is also considered an important indicator of regional climate change
[Robertson et al., 1992]. More recently, Bai et al., [2012] found that both North Atlantic
Oscillation (NAO) and El Nino-Southern Oscillation (ENSO) have impacts on Great Lakes ice
cover. Rao and Murty [1970] used a homogeneous, vertically integrated numerical model to
conclude that winter circulation in Lake Ontario was driven by wind. Understanding the major
effects of colder weather on the lake is crucial in predicting water movement patterns and the
thermal structure. In a study related to the trends of ice cover in the Great Lakes, Assel et al.,
[2003] pointed out that the distributions of ice cover among the five Laurentian Great Lakes
were influenced by the spatial pattern of water depth on each lake and regional differences in air
temperature. More recently, Fujisaki A.F et al., [2012] used the Princeton Ocean Model (POM)
which included ice processes to investigate the seasonal variation of ice cover on Lake Erie.
Their study showed that the packed ice cover slowed down the surface water velocity. To our
knowledge, simulation of lake process beneath the ice cover in three dimensions has not been

78

attempted. Therefore, there are large gaps in our understanding of the lake circulation,
temperature structure, the role of ice cover in Lake Huron during the winter time.
In this chapter, we investigate the evolution of ice cover, ice concentration, circulation,
and thermal structure over Lake Huron during winter time using a hydrodynamic model that is
coupled with an ice model. We used the unstructured grid FVCOM numerical ocean model
originally developed by Chen et al., [2003] coupled with ice module to describe the winter
circulation in Lake Huron during the winters of 2009-2010 and 2012-2013. The main reason for
implementing the model for the specified years is based on the available observational data.
FVCOM has been successfully applied in many previous studies for lakes or oceans but to date
the application of the model coupled with ice model has not been seen from previous studies. In
this chapter, the comparisons between the field data during the winter seasons, satellite remote
sensing data and model simulations of currents, temperature, and ice cover extent are presented.
The observed ice concentration data are downloaded from the Great Lakes Ice Atlas
(http://www.glerl.noaa.gov/data/ice/) and used to test with the simulated ice concentration. The
winter simulation results are also used to estimate the mean flushing times and residence times
for the Saginaw Bay that will be presented in Chapter 5. These novel results are expected to aid
in our understanding of how contaminants are flushed out of the bay during winter seasons and
how lake processes are changing in this important ecosystem.

79

4.2

Observations and Numerical Methods

4.2.1

Numerical Model

We use an unstructured-grid, three-dimensional hydrodynamic model of Lake Huron coupled
with a two dimensional ice model to understand the nature of circulation and the extent of ice
cover over the Lake Huron during winter seasons and to quantify exchange rates between the
Saginaw Bay and Lake Huron. The unstructured-grid ICE model was designed in the Cartesian
coordinate fully coupled to The Los Alamos Community Ice Code (CICE) [Hunke et al., 2010].
The computational domain to investigate the winter circulation and ice cover is the same as the
one that was used in the Chapter 3 which has 9611 nodes and 17619 triangular elements in
horizontal (Figure 4). The model has 21 vertical sigma levels. Based on the Courant- Friedrichs Levy (CFL) numerical stability criteria, the time steps are 4s and 10s for external and internal
scheme. The time steps for the ice model are 40s and 20s respectively.
As stated above, the hydrodynamic model is fully coupled with a two-dimensional ice
model. The model solves the hydrodynamic primitive equations (Equation 1 to Equation 5) using
the hydrostatic assumption in the vertical direction with the Boussinesq simplification for
convective flows. Vertical eddy viscosity and diffusivity are modeled using the Mellor-Yamada
2.5 level turbulence closure scheme [Mellor and Yamada, 1982; Galperin et al., 1988]. The
horizontal diffusion coefficients are calculated using the Smagorinsky turbulence closure model
[Smagorinsky, 1963].

80

In the presence of ice, the surface boundary conditions for momentum equations are given as
follows:

𝝏
𝝏𝒛

𝝏
𝝏𝒛

𝑲𝒎

𝝏𝒖𝒘

𝑲𝒎

𝝏𝒗𝒘

𝝏𝒛

𝝏𝒛

=

=

𝟏
𝝆𝒘

𝟏
𝝆𝒘

𝒄𝝉𝒘𝒙 + 𝟏 − 𝒄 𝝉𝒂𝒙

(33)

𝒄𝝉𝒘𝒚 + 𝟏 − 𝒄 𝝉𝒂𝒚

(34)

where c is the ice concentration (unit:dimensionless) , 𝜏𝑤𝑥 , 𝜏𝑤𝑦 are the water stresses in x and y
components, respectively. The movement of ice is driven by a combination of water flows and
wind as pointed out by Leppäranta, M. [2005].

Figure 26: Schematic of Ice Model. ci is the ice concentration; a is wind stress; w is water
stresss; uw is water current speeds; ui is the ice velocity; ua is wind speed; hi is the ice thickness.

81

Ice Dynamic Processes Equations
The movements of ice have an important role on lake circulation and thermal process. The
drifting of ice will lead to a redistribution of water temperature and it is influenced by a
combination of water movement and wind stress. The full nonlinear governing equations of the
ice dynamics processes are given as:

𝒎

𝝏𝒖𝒊

𝒎

𝝏𝒗𝒊

𝝏𝒕

𝝏𝒕

+ 𝒎 𝒖𝒊

𝝏𝒖𝒊

+ 𝒎 𝒖𝒊

𝝏𝒗𝒊

𝝏𝒙

𝝏𝒙

+ 𝒗𝒊

𝝏𝒖𝒊

+ 𝒗𝒊

𝝏𝒗𝒊

𝝏𝒚

𝝏𝒚

− 𝒎𝒇𝒗𝒊 =

+ 𝒎𝒇𝒖𝒊 =

𝝏𝝈𝟏𝒋
𝝏𝒙𝒋

𝝏𝝈𝟐𝒋
𝝏𝒙𝒋

− 𝒎𝒈

𝝏𝑯𝟎

− 𝒎𝒈

𝝏𝑯𝟎

𝝏𝒙

𝝏𝒚

+ 𝝉𝒂𝒙 + 𝝉𝒘𝒙

(35)

+ 𝝉𝒂𝒚 + 𝝉𝒘𝒚

(36)

where m is the mass of ice per unit area, ui, vi are the x and y components of the ice velocity,
respectively.H0 is the surface elevation 𝜎𝑖𝑗 is the internal stress tensor. 𝜏𝑎𝑥 , 𝜏𝑎𝑦 are the surface
wind stresses in x and y components, respectively.𝜏𝑤𝑥 , 𝜏𝑤𝑦 are the water stress in x and y
components, respectively.
The surface wind and water stresses are formulated as:

𝝉𝒂 = 𝒄𝒊 𝑪𝒂 𝝆𝒂 𝒖𝒂 𝒖𝒂 𝐜𝐨𝐬 𝝓 + 𝒌𝒖𝒂 𝐬𝐢𝐧 𝝓

𝝉𝒘 = 𝒄𝒊 𝑪𝒘 𝝆𝒘 𝒖𝒘 − 𝒖𝒊

𝒖𝒘 − 𝒖𝒊 𝐜𝐨𝐬 𝜽 + 𝒌 𝒖𝒘 − 𝒖𝒊 𝐬𝐢𝐧 𝜽

(37)

(38)

where ci is the ice concentration; Ca, Cw are the drag coefficient of air and water, respectively.
𝑢𝑖 is the ice velocity. 𝜌𝑎 , 𝜌𝑤 are the density of air and water, respectively. 𝜙, 𝜃 are the air and
water turning angles, respectively. k is unit vector.
82

The temperature condition at the water surface is specified as follow:

𝝏𝑻
𝝏𝒛

=

𝟏
𝝆𝒘 𝒄𝒑 𝑲𝒉

𝟏 − 𝒄𝒊 𝑭𝒂𝒊𝒓 − 𝟏 − 𝒄𝒊 𝑭𝑺𝑾𝒂 − 𝒄𝒊 𝑭𝒔𝒘𝒊𝒄𝒆 + 𝒄𝒊 𝑭𝑶𝑰 + (𝟏 − 𝒄𝒊 )𝑭𝒇𝒓𝒆𝒆𝒛𝒆

(39)

Where 𝐹𝑎𝑖𝑟 is the heat flux from the atmosphere to the lake; 𝐹𝑆𝑊𝑎 is the shortwave irradiance of
the open area; 𝐹𝑠𝑤𝑖𝑐𝑒 is the shortwave irradiance penetrated through ice; 𝐹𝑂𝐼 is the net heat flux
between the lake and ice; 𝐹𝑓𝑟𝑒𝑒𝑧𝑒 is the heat gained when ice grows over open water; ci is ice
concentration. Further details of the governing equations, boundary conditions, and
thermodynamics can be found in Chen et al. [2013]).
The ice transport equation is to describe the evolution of the ice thickness distribution:
𝝏𝒈 𝒙,𝒉,𝒕
𝝏𝒕

= −𝛁. 𝒈𝒖𝒊

(40)

where g(x,h,t) is the thickness distribution function containing the ice area, ice volume. ui is the
horizontal ice velocity. ∇= (

𝜕

,

𝜕

𝜕𝑥 𝜕𝑦

).

The ice model the distribution of ice thickness is sensitive to the simulation results [Hunke E.C,
2010]. In this study, we assumed there are five categories of ice thickness for each grid node
which is followed as the linear equation:

𝑯𝒊 = 𝒊 𝒅𝟏 + 𝒅𝟐 𝒊 − 𝟏

(41)

𝒅𝟏 = 𝟏. 𝟎 𝒏𝒄𝒂𝒕

(42)

𝒅𝟐 = 𝟎. 𝟓 𝒏𝒄𝒂𝒕

(43)

83

where Hi is ice thickness of ith category of ice; ncat is the number of category of ice thickness.
The thicknesses of ice in each category, therefore, are 20cm, 60cm, 120cm, 200cm, and 300cm
respectively. With the absence of ice cover, the water temperature of the lake is controlled by
exchange of heat flux between the surface water and atmospheric forcing, while with the
presence of ice cover, the heat flux is recalculated based on the ice cover conditions.
Winter seasons usually start on 23rd December in previous year and last until 22nd March
in the next year. In order to avoid confusion in the winter time definition, we use the next year as
the year of winter. The hydrodynamic model simulations were started in May 2009 and 2012
with zero velocity and uniform vertical temperature profile in order to avoid problems with
initialization of the stratified temperature. The initial temperature conditions were set up at 2.5
0

C and 3.5 0C for 2009 and 2012 respectively. The default setting in the model is used to

calculate wind stress and net heat flux which were then used in the Ice Model.

4.2.2

Meteorological and observational data

The same meteorological stations and data sources were used to download the
atmospheric data as described in the Chapter 2. The model used hourly atmospheric forcing data
such as wind speed, wind direction, surface air temperature, relative humidity, and cloud cover
fraction for the years 2009-2010, and 2012-2013. It must be noted that due to severe weather
conditions the surface water temperature was not observed at National Data Buoy Center
(NDBC) stations during winter. Therefore, we used the analysis of surface water temperature
from The Great Lakes Surface Environmental Analysis (GLSEA) to evaluate the model results.
The GLSEA data is a daily water temperature derived from NOAA polar-orbiting satellite
imagery obtained through the Great Lakes Coast Watch program. There are two thermistor
84

chains (SB32, SB33) deployed in the area of the mouth of Saginaw Bay (Figure 3a). The time
series of water temperature was observed through the water column. We also deployed three
bottom-mounted, up-looking acoustic Doppler current profilers (ADCPs) during the winter
season and collected hydrodynamic data to test the numerical models. The locations of the
ADCPs (SB32, SB33, and HM12) are presented on Figure 3. Note that the two years we
simulated are the mild winters in Great Lakes area compare to other years within 30 years
(Figure 27) where annual maximum ice coverage over Great Lakes are about 30% and 38.4% for
2010 and 2013 respectively. The long-term average is 51.4% as shown in Figure 27.

Figure 27: Great Lakes Annual Maximum Ice Coverage (1973-2013)
(http://www.glerl.noaa.gov/data/)

85

4.3

Results and Discussion

4.3.1

Results

4.3.1.1 Currents
The long time series of current measurements obtained in winters of 2010 and 2013 were
used to test the model results. Results from the hydrodynamic model coupled with two
dimensional ice model show that the lake-wide model with an average element size of 2.5 km is
able to resolve the circulation in the lake. Figures 28 and 29 show the comparisons between the
observed and simulated velocities (eastward and northward velocity) for the years 2010 and 2013
during winter time. The model results agree reasonably well with the observed data as shown in
Figures 28 and 29.

Figure 28: Vertically-Integrated velocity comparison at location HM12 with the model.

86

Figure 29: Vertically-Integrated velocity comparison at location SB32 (above) and SB33
(below) with the model.
87

It should be noted that ADCPs were deployed in difference areas in Lake Huron. For the
winter of 2010, ADCPs were deployed in the area near the mouth of Saginaw Bay in the
southern part of lake while for the winter of 2013 the ADCP was deployed in the area of
Hammond Bay in the northern part of lake. We concluded that the model was able to describe
large-scale circulation reasonably well. However, the simulated results for winter 2010 did not
capture some peaks of the water velocities at SB32 and SB33 stations during January as shown
in Figure 29 but for the long-term simulation the model was able to capture the water velocity
relatively well.
4.3.1.2 Temperature
Observations of water temperature at the two thermistor chain stations in the area of the Saginaw
Bay mouth were used to validate the model during winter 2010. Each temperature mooring
measured at every 4 meters vertically from July 2009 to the end of May 2010. In this chapter, we
present the comparisons between the simulation results and observations during the winter time
only (Figure 30). Accurate modeling of water temperature in vertical direction from surface to
bottom is an important test because it shows the accuracy of thermodynamics in the model. In
order to test the model, the time series of temperature distribution as a function of depth was
presented in Figure 30. The results show that simulated water temperature was in good
agreement with observed data. The model produced the evolution of temperature vertically from
stratified to homogeneous conditions relatively consistent with thermistor chain measurements at
both moorings SB32 and SB33 (Figure 30). The comparisons of vertical profiles of the
temperature clearly show that the hydrodynamic model coupled with ice model was able to
reproduce the important features of temperature in the water column in the period of transition
time from fall season to winter season.
88

Figure 30: Time series of observed versus simulated temperature in water column at SB32,
SB33 station

89

In addition to the comparisons at thermistor chain moorings, modeled water temperature were
also compared to lake-wide averaged surface temperature (GLSEA) derived from NOAA polarorbiting satellite imagery obtained through the Great Lakes Coast Watch program as shown in
Figure 31.

Figure 31: Comparison between lake-wide averaged simulated surface temperatures and
observationals data processed by GLSEA for (a) winter 2010 and (b) winter 2013. (The GLSEA
data extracted from http://coastwatch.glerl.noaa.gov/).

While comparisons with thermistor chain observations made at a point can be used to quantify
model accuracy, evolution of thermal process in Lake Huron and Saginaw Bay in winter season
can be understood by examining the spatial variability temperature which is shown in contour
plot in Figures 32, 33, 34, 35. At the end of fall (mid-November) the lake is almost homogenous
in temperature both vertically and horizontally with temperatures of 80C (Figure 30, 32, and 33).
As cooling of air continues in late fall, the water temperatures decrease starting in the shallow
areas such as the inner Saginaw Bay, some areas located in the Northern Bay and the Georgian

90

Bay. In March, Lake Huron is almost homogeneous at a temperature of 30C. The ice begins to
form when the water temperature decreases to freezing point.
The monthly average temperature map in Saginaw Bay was presented in Figures 34, 35. They
clearly show that the bay was at a homogeneous temperature in November and April. There is a
well-defined “tthermal boundary” located near the inner transect 1-1 in December. The
difference in temperature between the inner bay and outer bay is about 2.50C – 3.00C in
December. As cooling air continues to decrease in temperature the thermal boundary moves
lakeward in winter time until the bay reaches homogeneous temperature conditions.

91

Figure 32: Surface temperature contours for Lake Huron for the winter months 2010.

92

Figure 33: Surface temperature contours for Lake Huron for the winter months 2013.

93

Figure 34: Surface temperature contours for Saginaw Bay for the winter months 2010.

94

Figure 35: Surface temperature contours for Saginaw Bay for the winter months 2013.

95

The overturning of temperature in large lakes exhibits a seasonal variation. The lake-wide
average of temperatures at surface, mid-depth, and bottom through the winter seasons of 2010,
2013 are shown in Figure 36. The temperature at which the fall overturns is about 5 0C by midDecember and the spring overturns is about 4 0C by the end of the next April. The result is
consistent with the conclusions of Saylor and Miller [1979].

Figure 36: Lake-wide averaged of simulated water temperature at surface, mid-depth, and
bottom during winter 2010 (above) and 2013 (below).
96

4.3.1.3 Ice Cover, Ice Concentration, and Ice Thickness
The comparisons between modeled and observed ice cover are shown in Figures 38, 39. The
observed ice concentration data are downloaded from the Great Lakes Ice Atlas
(http://www.glerl.noaa.gov/data/ice/). The ice atlas products provide the daily ice concentration
data set. The observational ice cover data sets were used to compare with the simulated ice cover
extent. The distributions of observed ice concentration are reasonably well represented by the ice
model. Because over-lake precipitation and evaporation were not included in the ice model
differences between observed and modeled ice cover can be expected. Most of the ice cover in
Lake Huron occurred in Saginaw Bay, North Channel, and Georgian Bay during January and
February while there are large fractions of open water in the body of Lake Huron during the two
winters. The Saginaw Bay is almost entirely ice covered especially in the inner bay area with the
maximum ice concentration reaching 100%. The high concentration of ice will influence the
flushing rate between the outer and inner bays and between the (entire) bay and the open lake.
In Figure 37 we plot the time series of modeled ice cover in percentage (red and blue lines)
compared to observational data (red and blue solid dots) over Lake Huron for the winters of 2010
and 2013. The modeled percentage of ice cover is calculated as in the equation:

𝑷𝑰𝑪𝑬 =

𝒏
𝒊=𝟏 𝑨𝒊

𝑨𝒍𝒂𝒌𝒆

× 𝟏𝟎𝟎 %

(44)

where: PICE is the total ice cover expressed as a percentage; Ai is the area of cell which has ice;
Alake is the total area of Lake Huron and equal 5.88E+04 km2; n is the total number of cells which
have ice. It should be noted that it is assumed ice cover is fully filled at each cell where it has ice
formation at the centroid. This assumption may affect to the result of the total area covered by
97

ice in calculation where the grid size is large. Ice cover on Lake Huron started to form by
beginning of December for the years 2009, and 2012. The ice cover increased rapidly on Julian
day of 10 (2010) and 20 (2013) and it reached the maximum value of 38.3% and 38.6% in 2010
and 2013, respectively. Maximum values of simulated and observed ice areas are shown in Table
11. With the range of the annual maximum ice cover over Lake Huron varying from 45% to 79%
[Assel et al., 2003] the winters of 2010 and 2013 are considered mild winters. It is likely that the
model results have overestimated fractional ice cover as shown in Figure 37. The modeled ice
areas did not match with the observed data from mid-February to March while the model results
agree reasonably well with observations from January to mid-February. The difference between
simulations and observations areas covered by ice can be explained by the ice thickness
distribution used in the model. Because of the lacking of ice thickness data fields in Lake Huron,
therefore, in this study as stated before we assumed the ice thickness have 5 ice categories that
are 20cm, 60cm, 120cm, 200cm, and 300cm respectively while the distribution of ice thickness
is sensitive to the simulation results [Hunke E.C, 2010]. Furthermore, the exclusion of snow
cover may affect to the ice growth and melting processes because it is related to surface albedo.

Table 11: Maximum Ice Cover Area over Lake Huron (unit: %)
Modeled

Observed

2010

38.3

38.7

2013

38.6

45.6

98

Figure 37: The time series of ice cover in percentage in the winter 2010, 2013

99

Observed

Model

JD13-2010

JD13-2010

JD14-2010

JD14-2010

JD19-2010

JD19-2010

Figure 38: Map of Lake Huron showing observed and simulated ice cover area for winter 2010.

100

Figure 38 (cont’d)

Observed

Model

JD25-2010

JD25-2010

JD32-2010

JD32-2010

JD35-2010

JD35-2010

101

Figure 38 (cont’d)
Observed

Model

JD39-2010

JD39-2010

JD56-2010

JD56-2010

JD68-2010

JD68-2010

102

Observed

Model

JD 10 - 2013

JD 10 - 2013

JD 15 - 2013

JD 15 - 2013

JD 22 - 2013

JD 22 - 2013

Figure 39: Map of Lake Huron showing observed and simulated ice cover area for winter 2013.

103

Figure 39 (cont’d)

Observed

Model

JD 24 - 2013

JD 24 - 2013

JD 40 - 2013

JD 40 - 2013

JD 48 - 2013

JD 48 - 2013

104

In Figures 40, 41 we plot the contours of modeled ice thickness for some selected days in
January, February, and March for winters of 2010, 2013. Mean ice thickness over the three
months of winter for the two years was 3.4 cm for Lake Huron and 7.1 cm for the Saginaw Bay.
The results of ice thickness in winter months for Lake Huron and Saginaw Bay for the year 2010
and 2013 are presented in Tables 12, 13. The maximum of ice thickness was reached in February
with a value of 95.4 cm for Lake Huron. In Saginaw Bay the ice thickness reached the maximum
value of 35.4 cm in March 2013.

Table 12: Monthly average ice thickness (in cm) for Lake Huron and
Saginaw Bay in winter months
Year

Areas

Jan

Feb

Mar

Lake Huron

1.7

4.8

1.2

Saginaw Bay

4.5

11.5

2.7

Lake Huron

1.5

6.2

5.0

Saginaw Bay

3.7

11.8

8.8

2010

2013

Table 13: Maximum ice thickness (in cm) for Lake Huron and Saginaw
Bay in winter months
Year

Areas

Jan

Feb

Mar

Lake Huron

60

53.5

36.0

Saginaw Bay

21.7

22.5

22.4

Lake Huron

42.2

95.4

40.0

Saginaw Bay

30.4

23.4

35.4

2010

2013

105

JD13-2010

JD32-2010

JD19-2010

JD39-2010

JD25-2010

JD56-2010

Figure 40: Modeled Ice thickness contours for Lake Huron for the winter of 2010.
106

JD 22 - 2013

JD 63 - 2013

JD 40 - 2013

JD 67 - 2013

JD 48 - 2013

JD 76 – 2013

Figure 41: Modeled Ice thickness contours for Lake Huron for the winter of 2013.
107

4.3.2

Discussion

During the winter season, large-scale circulation in Lake Huron is influenced by a combination
of key physical features (islands, exchange with Saginaw Bay and Georgian Bay, and the midlake AAR ridge) and the movement of ice. With the presence of ice, the results from numerical
model have shown that winter circulation in the main lake is counterclockwise. The mean current
speed in Lake Huron was 1.14 cm/s over winter months for the two simulation years. The mean
vertical velocity profiles at the inlet of Georgian Bay are illustrated in Figure 42, in which the
positive and negative values indicate inflow and outflow currents to and from Georgian Bay,
respectively. The velocity profiles in winter time were unusual in February and March in
comparison with the earlier reported results by Beletsky et al., [1999]. The results show that there
is an outflow at the surface associated with the inflow at the bottom at Georgian Bay mouth. The
mean velocity of inflow into Georgian Bay in deeper layers is about 1 cm/s. The presence of
flow reversal at Georgian Bay mouth can be explained by the role of dominant wind direction
under homogeneous conditions. As shown in Figure 43 the predominant wind directions over
Lake Huron in January and February are SW and NE respectively. Furthermore, the circulation
can change from year to year or season to season. Circulation in Saginaw Bay was characterized
by the presence of a cyclonic gyre at the mouth of the Saginaw bay and two recirculating cells
within the inner bay as shown in Figures 44, 45. The role of ice cover on the current velocity was
clearly found as presented in Figures 44, 45 with low current velocities in the inner bay
especially in February. This finding is supported by the conclusion from Fujisaki et al., [2012]
who concluded that the packed ice cover slowed down the water velocity. This finding also has
important implications in that it indicates that the flushing time in winter is expected to be longer
than the flushing time in summer for the inner bay. The reduction of current velocities in winter
108

implies lower mixing rates in Saginaw Bay. Another important feature of circulation in winter
time that there is a narrow band of flow of about 2km width along the mid-lake ridge from
Thunder Bay at Alpena on the west side across to Point Clark on the east side.

Figure 42: Mean vertical velocity profiles at the inlet of Georgian Bay for different layers in
winter months.

109

Table 14: Maximum and mean currents velocity in Lake Huron in winter months (unit: cm/s)
Year

December

January

February

March

Mean

Max

Mean

Max

Mean

Max

Mean

Max

2010

1.2

12.0

1.6

9.7

1.4

15.0

1.3

8.3

2013

0.72

4.5

1.35

10.45

0.6

4.93

1.0

7.4

Figure 43: Wind rose plots for winter months over Lake Huron during in the years 2010, 2013.

110

Figure 44: Monthly vertically-averaged currents in Saginaw Bay in winter months 2010
111

Figure 45: Monthly vertically-averaged currents in Saginaw Bay in winter months 2013
112

Figure 46: Monthly vertically-averaged currents in Lake Huron in winter months 2010.

113

Figure 47: Monthly vertically-averaged currents in Lake Huron in winter months 2013.
114

As shown by the observed and simulated vertical temperature profiles presented in Figure 30
the lake is almost entirely homogeneous in temperature by mid-November. During the
homogeneous period, we can see that the thermal boundary between the inner and outer bays
in Saginaw Bay where there is a significant gradient in temperature between the two zones.
The presence of this thermal temperature can be associated with the bathymetric features in
the bay. Since the mean depth in the inner bay is about 4m, the water temperature decreases
relatively quickly than in the deeper layers of the outer bay as cooling of ambient air
continues to occur.
The monthly vertically-averaged currents in Figures 44, 45 clearly show the role of ice cover
on water movement underneath ice, especially in the inner area in the Saginaw Bay. It is
interesting to note that the inner bay almost reached a state of stagnation in February (Figure
44b). The flushing times for inner and entire bay, therefore, are expected to be higher than
the flushing times during the summer season. It is likely that the model results have
overestimated fractional ice cover over Lake Huron. The ice model needs further
improvements in future for more accurate simulation of ice cover. The monthly average ice
velocities in Lake Huron for the years 2010 and 2013 are illustrated in Figures 48, 49. In
Georgian Bay the ice packs tend to move in the Northeastern direction. It is likely that in
Georgian Bay the direction of ice movement is not strongly influenced by the wind direction
but the combination of wind stress and water movement. In Saginaw Bay the ice movement
is associated with the circulation in the Bay (Figures 50, 51). We can see that the ice
movement tends to follow the two re-circulating cells within the inner bay in March 2010
and January 2013.

115

Figure 48: Monthly average Ice velocity in Lake Huron in 2010.
116

Figure 49: Monthly average Ice velocity in Lake Huron in 2013.
117

Figure 50: Monthly average Ice velocity in Saginaw Bay in 2010.

118

Figure 51: Monthly average Ice velocity in Saginaw Bay in 2013.

119

4.4

Conclusions

The unstructured grid hydrodynamic model coupled with an ice model is used to investigate the
large-scale circulation, thermal structure, and ice cover in Lake Huron and Saginaw Bay.
Simulation is carried out for two mild winters of 2010 and 2013. The model included multiple
ice categories. The model reproduced reasonably well, the lake surface temperature, circulation,
ice cover, ice concentration, and ice thickness in Lake Huron. The model results were tested
against ADCP observations of currents, thermistor chain data, and observed ice cover data. Mean
circulation in winter time was predominantly cyclonic in the main body of Lake Huron with
current speeds being highest in January. Thermal structure in winter was characterized by a
homogeneous condition that was found by both observations and simulations in the study. The
lake is almost entirely homogeneous in temperature by mid-November. The thermal boundary
between the inner and outer bays in Saginaw Bay was found during winter time with a difference
in temperature of approximately 2.50 to 3.00C. The important features of fall and spring overturn
were also reproduced by the model. As expected, the model results show that the significant role
of ice cover on lake circulation in winter. The inner area in Saginaw Bay was almost in a
stagnation state during the ice formation period. As a result, the flushing times are expected to be
higher than the corresponding flushing times in summer season (this question will be explored in
detail in the next chapter). We found that the mean current speed in Lake Huron in winter time
was 1.14 cm/s compared to the summer time average of 3 cm/s found in an earlier chapter. The
remarkable feature we found is the unusual nature of currents at the inlet of the Georgian Bay.
There is an outflow at the surface layer associated with the inflow at the bottom at the Georgian
Bay mouth. The model results show that ice cover reached a maximum of 38.3% and 38.7% in
mid-February in 2010 and beginning of March in 2013 respectively. We found that the mean ice
120

thickness over the three months in winter for 2 years was 3.4 cm for Lake Huron and 7.1 cm for
the Saginaw Bay.

121

Chapter 5
Residence Time in the Saginaw Bay and Lake
Huron System
5.1

Introduction

Agricultural, industrial, and municipal activities, both within the Lake Huron basin and in
upstream region of the Saginaw River, have resulted in pollution by a variety contaminants and
the subsequent degradation of water quality. Persistent contamination of the water column and
bottom sediment has resulted in degraded water quality. The embayed regions tend to
accumulate contaminants due to their longer residence times as indicated by Nixon et al., [1996].
Residence time and flushing time of an estuary are important variables used to estimate
the rate of removal of pollutants from the bay. Therefore, a spatial map of the residence time and
flushing time of the bay can be useful. Understanding the mean residence time and spatial
variation in the bay is of critical importance in resolving issues such as contamination, water
quality (e.g., harmful algal blooms, bacterial contamination etc.), and for ecosystem analysis;
therefore, residence time is an important index that scientists should pay careful attention.
Monsen et al., [2002] concluded that first-order transport time scales such as residence times and

122

flushing times within the bay are useful for understanding, interpreting, and comparing processes
and rates across different ecosystems. For example, when attempting to solve issues related to
the transport of matter in lakes or a bay, we need to know how the exchange takes place. Recent
concerns associated with environmental issues in the Saginaw Bay highlight the need for
accurate estimates of transport time scales.
There have been a few studies into the Saginaw Bay and Lake Huron systems that have
provided the preliminary description of the mean residence times in these areas. Quinn [1992]
estimated that the residence time, which is defined as the time it takes to replace the water
volume of Lake Huron, is 22 years. Saylor and Danek, [1976] calculated the exchange rate
between the inner bay and outer bay to be 3700m3/s, while the mean residence time was 26.5
days for the inner bay. The residence time of any bay may vary with several factors such as wind
direction or river inflow. Since the Saginaw bay is shallow, the residence time is strongly
influenced by wind. As Saylor and Danek, [1976] pointed out, the residence time of water in the
bay is much longer when the winds are perpendicular to the axis of the bay. Earlier study by
Dolan [1975] calculated a residence time of 110 days for the inner bay based on a time-variable
chloride model. However, the river discharge was not included in this model.
Until now, there have been no observations and studies on residence times in the Saginaw
Bay since the early work cited above. High resolution hydrodynamic models with additional
water quality models can be used to estimate the residence times in more detail than with an
observational dataset alone. Therefore, the aim of this chapter is to estimate the residence time
for Saginaw Bay and to describe the residence time distribution zones using two approaches.
Methodologies for computing flushing time and residence time will be presented in the next
session. The distribution map will provide “dead-zones” within the Bay that will increase the

123

capability to solve and manage environmental problems. The contribution of Saginaw River
discharge as a source of fresh water to the residence time in the entire and the inner bays will be
quantified and the relationship between them will be established.

5.2

Methodology

5.2.1 First-Order Transport Time Scales
The terms, residence time and flushing time, were used in different studies but sometimes there
was confusion in the definition [Sheldon and Alber, 2002]. One of the simplest and most widely
used transport time scales is the flushing time (𝑇𝑓 ) defined as “the ratio of the mass of a scalar in
a reservoir to the rate of renewal of the scalar” [Geyer, 1997; Monsen et al., 2002]. Flushing time
can be calculated by dividing the volume of the bay (V) by the volumetric flow rate into (𝑄𝑖𝑛 ) or
out of (𝑄𝑜𝑢𝑡 ) the bay (assumed to be equal in this idealized case based on the assumptions of
steady-state, and constant-volume):

𝑻𝒇 = 𝑽 𝑸
𝒊𝒏

(45)

The volume V can be calculated from bathymetric data over the surface area of the bay. The
flushing time as defined here was variously referred to in the past as “residence time” [Go´mezGesteira et al., 2003; Chapra 2008], “turnover time” [Sheldon and Alber, 2006; Takeoka, 1984],
“flushing time” [Fischer et al., 1979; Monsen et al., 2002; Delhez et al., 2004], “water exchange
rate” [Kraines et al., 2001] and “water renewal time” [Andrefouet et al.,2001]. The flushing
time is an important physical index for determining the capability of the bay itself in term of

124

flushing out the pollutant. Equation 45 shows that the flushing time is directly proportional to the
value of the volume of the bay. In a water body with a long flushing time index it will take long
time to flush out the pollutant. In other words, a long flushing time may be a potential condition
to increase concentration and accumulate the pollutant within a semi-enclosed water body like
Saginaw Bay.

5.2.2 Eulerian Method
Another important transport timescale is the residence time (𝜏𝑅 ), which is a measure of “how
long a water parcel, starting at a specified location within the water body, will stay within the
boundaries of the water body before exiting” [Monsen et al., 2002]. By treating the bay as a
CSTR (continuously stirred tank reactor, [Chapra, 2008]), we can calculate the residence time
for every grid cell in the domain as the e-folding flushing time [Monsen et al., 2002] – that is, the
time at which 36.78% (𝑒 −1 ) of the initial mass still remains within the bay. As pointed out by
[Monsen et al., 2002], the flushing time (𝑇𝑓 ) is an integrative system measure while the residence
time (𝜏𝑅 ) is a local measure that varies from one grid cell to the other.
The dye concentration model was used to define the residence time (𝜏𝑅 ) as the time it takes
for the vertically integrated dye concentration at every grid cell to drop below the 1/e value
(36.78%).
The dye concentration was solved using the following equations:

𝝏𝑫𝑪
𝝏𝒕

+

𝝏𝑫𝒖𝑪
𝝏𝒙

+

𝝏𝑫𝒗𝑪
𝝏𝒚

+

𝝏𝒘𝑪
𝝏𝝈

−

𝟏 𝝏
𝑫 𝝏𝝈

𝑲𝒉

𝝏𝑪
𝝏𝝈

125

− 𝑫𝑭𝒄 = 𝑫𝑪𝟎 𝒙, 𝒚, 𝝈, 𝒕

(46)

where C is the concentration of the dye, D is the total depth, u,v, and w are the x,y, and 
components of the water velocity, Kh is the vertical diffusion coefficient, Fc is the horizontal
diffusion term, and C0 is the initial dye concentration.
The dye concentration is placed at every grid cell throughout all layers in the water column. The
residence time within the bay is then calculated as the time it takes for each element’s verticallyintegrated concentration to drop as an exponential equation from an initial condition of 100 ppm
to 36.78 ppm [Burwell et al., 2000]. A High resolution numerical model is expected to provide
accurate estimates of residence time in the bay and that is one of the objectives of this chapter.
More details will be presented in the next session.

5.3

Results and Discussion

5.3.1

Results

5.3.1.1 First-Order Transport Time Scales

Due to the recirculating nature of the flow within the bay, if we consider a transect separating the
bay from the rest of the lake (as shown in Figure 3), then flow enters the bay over a portion of
the transect and exits over the remaining portion (Figure 18, 19). The relative extent of these
inflow and outflow regions as well as the magnitude and direction of velocities along the
transect exhibit unsteadiness throughout the simulation period. These complexities introduce
difficulties in the application of the simple flushing time concept.

126

The exchanges between the inner and the outer bays and between the outer bay and the
open lake fundamentally describe the mean residence time and contaminant flushing rates in
Saginaw Bay. The complexity of flows and exchange rates can be seen in Figure 19, Figure 20
which shows the results of vertically-averaged currents from the LWM superimposed on surface
temperatures in Saginaw Bay for a three year simulation period. The presence of two large
recirculating cells in the inner bay can be clearly noticed (e.g., in July and August for all three
years) and the dominant flow features can be related to the predominant wind patterns. Shown as
wind rose plots in Figure 51 the length of the histogram denotes the frequency of winds coming
from a certain direction while the color denotes the magnitude of the wind velocity.
Equation 47 was used to calculate the mean volumetric inflow into Saginaw Bay at the
inner (1-1) and outer (2-2) transects shown in Figure 3a to facilitate the application of the
flushing time concept. Horizontal velocities from the lake-wide model were interpolated to
points on the Saginaw Bay transects and used in Equation 47 to calculate the mean inflow 𝑄𝑖𝑛
volume during the summer and winter months. The inflow volume at inner and outer transects
were used to estimate the flushing rates for the inner and outer bay [Saylor et al., 1976]. The
flushing time (𝑇𝑓 ) for the inner (and entire) bay are presented in Table 15 for summer months.

𝑸𝒊𝒏 =

𝒕𝟐
𝒕
(𝒕𝟐 −𝒕𝟏 ) 𝟏
𝟏

𝑣𝑖𝑛 =

𝑩
𝟎

𝑯
𝒗𝒊𝒏 𝒅𝒉
𝟎

𝑣. 𝑛
0

𝒅𝒃 𝒅𝒕

(47)

if 𝑣. 𝑛 > 0
if 𝑣. 𝑛 ≤ 0

Where 𝑣 is the velocity (m/s) from the hydrodynamic model, 𝐵 is the length of the inner (or
outer) transect, 𝐻 is the depth, 𝑡1 , 𝑡2 are the start and end times for the calculations, and 𝑛 is a unit
127

vector normal to the transect pointing into the bay (towards southwest). V is the volume of the
inner (or entire) bay.
Figure 52 shows the vertically averaged volumetric influx as a function of time at the inner
(and outer) transects for the summer months of 2009-2011. A significant difference between
exchange at the outer transect 2-2 and at the inner transect 1-1 was found. Another notable
feature is that the exchange at the outer transect has higher fluctuation compared to the one at the
inner transect.

Figure 52: Time series of water exchange rates for the inner transect (1-1 in Figure 3) and the
outer transect (2-2 in Figure 3) in Saginaw Bay in summer season.
128

Table 15: Mean flushing time 𝑻𝒇 (in days) for inner and entire Saginaw Bay in summer
months
Year

Areas

July

August

September

Inner Bay

23.9

26.2

19.7

Entire Bay

8.8

13.9

9.1

Inner Bay

26.3

25.0

19.5

Entire Bay

9.1

8.9

8.6

Inner Bay

25.1

24.4

16.8

Entire Bay

11.6

10.1

9.2

2009

2010

2011

129

Figure 53: Wind rose plots for summer months during the years 2009-2011.

130

As stated in chapter 3, the circulation in Saginaw Bay is complex with either clockwise or
counterclockwise circulation depending on the wind direction. The important feature of
circulation in the bay was also characterized by the difference between hydrodynamics in the
inner and outer bays and the effect on exchange between the two. The complexity of the
circulation suggests that the hydrodynamics in both the inner and outer bays plays a significant
role in the flushing time distribution. Here, the flushing time varies from one grid cell to other.
In order to find “Dead-Zones” and “Flusing-time Maps” we calculated the flushing time
based on the results from the hydrodynamic model for every grid cell within the bay. Following
this approach, the eastward and northward velocities obtained from the hydrodynamic model
were interpolated into a high resolution structured grid of 200mx200m. Interpolated velocities
were then used to calculate the mean volume inflow into each box with consideration of flux
between adjacent boxes (Figure 54).

Qin
Qout

Vi

Qout

Qin

Figure 54: Example of a diagram depicting the water exchanges in a structured grid.
The flushing time for each structured grid cell can be calculated as the volume of each cell
divided by the rate of inflow:
131

𝝉𝒊 =

𝑽𝒊

(48)

𝒒𝒊

Where: 𝜏𝑖 ,𝑉𝑖 , and 𝑞𝑖 are the flushing time, volume, and volumetric flow rate into box ith,
respectively. It should be noted that the calculation is based on the assumptions of steady-state,
and constant-volume.
The flushing times were calculated for each summer month for 3 constitutive years (20092011) based on vertically integrated currents in the Saginaw Bay. The spatial map in Figure 55
shows the flushing time is mostly closely associated with the bathymetry map. The shallower
areas have a shorter flushing time. The variability in estimates of flushing times within the bay
can have important implications for stagnation areas between two large recirculating cells within
the inner bay. It is important to recognize that the method of obtaining the total flushing time for
the bay by accumulating flushing times from individual cells cannot be relied upon.

132

July - 2009

August - 2009

September - 2009

July - 2010

August - 2010

September - 2010

July - 2011

August - 2011

September - 2011

Figure 55: Map of Saginaw Bay showing the contour of flushing time in hours within the Bay
calculated based on monthly average velocity.

133

Circulation in the Saginaw Bay in winter months was characterized by two large recirculating
cells in the inner bay just as is observed in the summer months. The same approaches and
equations were applied to calculate the water exchange and flushing time for the bay for winter
in the years 2010 and 2013. Since the vertically integrated current velocity in winter is lower
than in summer as shown in Figure 46, Figure 47, the flushing times are expected to be longer
than in summer time. The vertically averaged volumetric inflow as a function of time in winter at
the transect 1-1 and transect 2-2 are presented in Figure 56. The complexity of the exchange
rates at outer transect in winter time can be seen by the fluctuation of volumetric influx and the
difference of this index between the two winters (Figure56b).

Table 16: Mean flushing time 𝑻𝒇 (in days) for inner and entire Saginaw Bay in winter months
Year

Areas

Dec-2009/2012

Jan

Feb

Mar

Inner Bay

15.2

33.9

108.0

18.1

Entire Bay

13.9

11.7

18.5

11.2

Inner Bay

20.8

27.0

75.7

46.9

Entire Bay

18.8

19.7

16.6

14.5

2010

2013

134

Figure 56: Time series of water exchange rates for the inner transect (1-1 in Figure 3) and the
outer transect (2-2 in Figure 3) in Saginaw Bay in winter season.

135

5.3.1.2 The Eulerian Approach

In this approach, the bay was treated as a CSTR (continuously stirred tank reactor, [Chapra,
2008]). The residence time for every grid cell in the domain is then calculated as the time it takes
for only 36.78% (e-1) of the initial mass to remain within the bay. It is important to recognize
that, due to zero-dimensionality and instantaneous mixing assumptions, the ideal CSTR case do
not exactly represent real systems. However, the measures are still useful for understanding and
interpreting the behavior of contaminants within the bay.
The same hydrodynamic model used for calculating the flushing time (𝑻𝒇 ) was used to
calculate the residence time (expressed as the e-folding flushing time) assuming that the Saginaw
Bay was a CSTR. Here it is possible to use different approaches as described in [Burwell et al.,
2000]. A Lagrangian approach in which particles are released and tracked in every grid cell and
an Eulerian approach in which the residence time is calculated based on dye concentration
modeling are both options. We have used the second approach.
In this method, the dye concentration was calculated following the instantaneous dye
release at an initial concentration of C0 = 100 ppm in every grid cell in the bay area. The
residence time is defined as the time it takes for each element’s vertically-integrated
concentration to decrease below 1/e of its initial vertically-integrated concentration [Burwell et
al., 2000]. Since the residence time (𝝉𝑹 ) calculated using this definition would be different for
every grid cell, an integrative system measure (𝝉𝑹 ) can be obtained by computing the spatial
average of the vertically-integrated concentrations for all grid cells in the domain. This average
concentration is only a function of time and does not depend on space. The residence time for the
bay (𝝉𝑹 ) is then calculated as the time it takes for the average vertically-integrated concentration

136

in the bay to fall below (1/e) of its initial concentration (i.e., below an average (C/C0) value of
0.3678).
To understand the effect of the Saginaw River on the residence time, we calculate
residence times with and without river flow entering the bay. For each case, the hydrodynamic
model was run with flow (or no inflow) from the river and an initial concentration of 100 ppm in
every grid cell. For the case with river inflow, the bay was subjected to the initial concentration
of 100 ppm but the river water was assumed to be at a zero concentration; therefore, the
residence time estimated for this case represents the effect of flushing of the initial water in the
bay by the river. It is important to recognize that flushing time and residence time as calculated
here can be expected to yield different estimates of transport timescales since the flushing time
concept assumes that advective exchange is the only mechanism of transport, while the
concentration-based residence time estimate considers the processes of both advection and
dispersion in the modeling.
In order to investigate the effect of the ice cover on the residence time, we also calculated
the residence times with and without ice model. The results will provide the ranges of the
residence time for the winter with less ice coverage to the case of the severe winter with high ice
cover.
The summer contour map of residence times in the inner bay and the entire bay with and
without river inflow and the relationship between the vertically-integrated concentration
averaged over the entire and inner bay and the residence time are presented in Figures 57 and 58
respectively. Similar contour maps for winter season are shown in Figure 59 and 60 for the case
of without ice model and Figure 61 and 62 for the case of with ice model.

137

Figure 57: Contour plots of residence time for Saginaw Bay based on dye concentration
modeling. Figures (a) and (b) are for inner bay dye releases with (a) no river inflow and (b) with
river inflow. Figures (c) and (d) are results based on dye releases for the entire bay with (c) no
river flow and (d) with river flow. (Summer season).

138

Figure 58: Average vertically-integrated concentration in the Saginaw Bay as a function of time.
The red solid lines represent the model results while the blue dashed lines represent an
exponential decay model fit to the data. The concentration value corresponding to a value of
(1/e) is also marked. (Summer season).

139

The empirical relation curves between the residence time and average vertically-integrated
concentration for the inner and entire bay in summer season are as follow:

For Inner Bay:
With River Inflow: (R2 = 0.85)
𝑪
𝑪𝟎

= 𝟎. 𝟕𝟒𝟑𝒆−𝟎.𝟎𝟏𝟎𝟗𝟓𝒕

(49)

Without River Inflow: (R2 = 0.79)
𝑪
𝑪𝟎

= 𝟎. 𝟕𝟔𝟕𝒆−𝟎.𝟎𝟎𝟖𝟎𝟑𝒕

(50)

For Entire Bay:
With River Inflow: (R2 = 0.95)
𝑪
𝑪𝟎

= 𝟎. 𝟖𝟖𝟗𝒆−𝟎.𝟎𝟎𝟖𝟏𝟔𝒕

(51)

Without River Inflow: (R2 = 0.93)
𝑪
𝑪𝟎

= 𝟎. 𝟗𝟐𝟔𝒆−𝟎.𝟎𝟎𝟔𝟐𝟒𝒕

(52)

where t is the residence time in day.
Table 17: Mean residence time (e-folding flushing time) in days for Inner and Entire Saginaw
Bay in summer season
Year

2011

Areas

Case
No river inflow

With river inflow

Inner Bay

112

62

Entire Bay

125

115

140

Figure 59: Contour plots of residence time for Saginaw Bay based on dye concentration
modeling. Figures (a) and (b) are for inner bay dye releases with (a) no river inflow and (b) with
river inflow. Figures (c) and (d) are results based on dye releases for the entire bay with (c) no
river flow and (d) with river flow. (Winter season without Ice Model).

141

Figure 60: Average vertically-integrated concentration in the Saginaw Bay as a function of time.
The red solid lines represent the model results while the blue dashed lines represent an
exponential decay model fit to the data. The concentration value corresponding to a value of
(1/e) is also marked. (Winter season without Ice Model).

142

The empirical relations curves for winter season for the case of without ice model are as follow:
For Inner Bay:

With River Inflow: (R2 = 0.985)
𝑪
𝑪𝟎

= 𝟎. 𝟗𝟖𝟗𝒆−𝟎.𝟎𝟏𝟕𝟐𝟗𝒕

(53)

Without River Inflow: (R2 = 0.98)
𝑪
𝑪𝟎

= 𝟏. 𝟎𝟑𝒆−𝟎.𝟎𝟏𝟓𝟔𝟑𝒕

(54)

For Entire Bay:
With River Inflow: (R2 = 0.952)
𝑪
𝑪𝟎

= 𝟏. 𝟏𝟏𝒆−𝟎.𝟎𝟎𝟗𝟕𝒕

(55)

Without River Inflow: (R2 = 0.935)
𝑪
𝑪𝟎

= 𝟏. 𝟏𝟔𝒆−𝟎.𝟎𝟎𝟖𝟕𝟗𝒕

(56)

Table 18: Mean residence time (e-folding flushing time) in days for Inner and Entire Saginaw
Bay in winter season (Without Ice Model)
Year

2013

Areas

Case
No river inflow

With river inflow

Inner Bay

67.7

64.7

Entire Bay

131.8

114.2

143

Figure 61: Contour plots of residence time for Saginaw Bay based on dye concentration
modeling. Figures (a) and (b) are for inner bay dye releases with (a) no river inflow and (b) with
river inflow. Figures (c) and (d) are results based on dye releases for the entire bay with (c) no
river flow and (d) with river flow. (Winter season with Ice Model).
144

Figure 62: Average vertically-integrated concentration in the Saginaw Bay as a function of time.
The red solid lines represent the model results while the blue dashed lines represent an
exponential decay model fit to the data. The concentration value corresponding to a value of
(1/e) is also marked. (Winter season with Ice Model).

145

The relations between the residence time and average vertically-integrated concentration in
winter season for the inner and entire bay for the case of with ice model are shown in Figure 62.
The results are quite different from the case of without ice model. For the case of with ice model,
after about 226 days, the mean vertically-integrated concentration in the inner bay drop to 95%
and 93.4% with no river inflow and with river inflow, respectively. Corresponding number for
the entire bay are 71% and 69.8% after 156 days and 154.7 days, respectively.
Considering the residence time including mixing processes revealed the behavior of the bay
under real meteorological conditions. Considerably different residence times were found for the
cases with the river and without the river. This leads to the suggestion of the relationship
between the residence time and the river discharge that will be presented in the next session.

5.3.2

Discussion

The mean residence time for contaminants entering the Saginaw Bay from point and non-point
sources in the bay was estimated from two different methods. The first method used results from
the lake-wide numerical model. The flushing rates for water (and dissolved contaminants) in the
inner and outer bay were calculated using Equation (45) for different summer months during
2009-2011 and winter months during 2010, 2013 are presented in Table 15 and Table 16,
respectively.
The mean flushing time for all three summers (and all 3 years) was 23.0 days for the
inner bay and 9.9 days for the entire bay. In an earlier study, Saylor et al. [1976] estimated a
mean flushing time of 26.5 days for the inner bay. This difference is explained by the change in
volume of the inner bay (from 8.5 km3 in 1976 to 6.5 km3 in 2012) over the past four decades

146

due to falling lake levels. The results presented in Table 15 clearly show that flushing rate for
the outer bay is almost twice as much as for the inner bay. This could be due to the role played
by lake-wide circulation on flow in the outer bay. Currents within the Saginaw Bay are strongly
affected by the direction of wind. Wind blowing out of the southwest or northeast is associated
with higher mean current speeds and therefore, greater mixing and dilution as well as a smaller
mean residence time. No such trend is apparent for the mean flushing rate for the outer bay,
which is affected to a greater degree by the large-scale circulation. Exchange at the outer transect
2-2 (Figure 52a) is significantly higher compared to the exchange at the inner transect 1-1
(Figure 52b). We also notice that the rate of exchange between the inner and outer bays (Figure
52b) shows a significantly larger variation in summer time (varies by nearly four orders of
magnitude) compared to the exchange between the outer bay and the lake (Figure 52a). One of
the reasons for this large difference in the variability of exchange rates is the widely different
cross-sectional areas at sections 1-1 and 2-2 in Figure 3a as well as the presence of an
anticyclonic gyre in the outer bay (shown in Figure 19). The cross-sectional area for the inner
bay transect 1-1 is approximately 0.09 km2 while the area is 1.31 km2 at the outer transect 2-2.
While in summer months the flushing time in the inner bay is nearly twice as much for
the entire bay, the flushing time is almost similar in the inner and entire bay in the beginning of
winter (in December) (Table 16). As the cooling air continues and with the presence of ice in
January, February, and March there is a significant difference for the flushing time between the
inner bay and entire bay. The maximum value of flushing time for the inner bay reached to 108
days in 2010. The difference can be explained by the role of ice cover during the winter time
especially in inner bay area. The mean flushing time for two winters was 43.2 days for the inner
bay and 15.6 days for the entire bay.

147

Using the First-Order transport time scale approach we calculated the flushing time for
each structured grid cell within the bay as described in Figure 54 and Equation 48. The results
show that as the month progresses in summer, the spatially flushing time distribution does not
change except for the year 2011 (Figure 55). The areas with higher values of flushing times were
found in the inner bay. It is clearly shown that the role of two gyres within the inner bay causes
the occurrence of stagnation areas where there are longer flushing times.
Using an Eulerian dye concentration-based approach, we calculated the mean residence
time (e-folding flushing time) in the bay as described earlier for both summer and winter
seasons. The residence time results are shown in Table 17 and Table 18 for summer and winter
season respectively. For summer season, the range of values is from 62 days to 125 days
depending on whether inflow from the Saginaw River is considered or not. An earlier study by
Dolan [1975] calculated a residence time of 110 days for the inner bay and 52 days for the entire
bay based on a time-variable chloride model. Their results, based on “well-mixed” mass balance
models [Chapra, 2008], allow the bay to be divided into a few (e.g.,three) segments taking both
advective and dispersive exchange across segments into account. Since the present work is based
on fully three-dimensional hydrodynamic and transport models, our residence time estimates can
be expected to differ from those of Dolan [1975]. The results presented in Table 17 and Table 18
show that residence times calculated by the Eulerian dye concentration-based model are greater
than the mean flushing times calculated using the first method which considers only advective
exchange. This difference can be explained by the increased time due to mixing and dilution in
the entire bay.
For summer months, the low residence times occur in the middle of the bay, and
high residence times occur along the shoreline and the Saginaw river mouth areas (Figure

148

57). With dye release into the entire bay, low residence times were found in the areas near the
mouth of the Saginaw Bay (5 to 50 days) while the remaining areas have relatively high
residence times (more than 140 days, Figure 57). The figure clearly shows the important role
played by the gyres near the mouth of the outer Saginaw Bay on residence times. As can be seen
from Figure 57, river inflow has a significant effect on the residence time of the inner bay
although the effect is almost insignificant when the entire bay is considered (Figure 57d). This is
not surprising since the inner bay is fairly shallow and is readily influenced by riverine discharge
while the gyres and the larger lake-wide circulation mainly influence the outer bay residence
time. Figure 58 shows the vertically-integrated concentrations within the bay (averaged over the
spatial extent of the bay) as a function of time for the different cases considered (inner bay
and entire bay with and without river flow). The figure shows the models results and the
best exponential-fit to the models results. Generally, the declining concentrations follow an
exponential curve with R2 values ranging from 0.79 to 0.95. The empirical equations between
the residence time and average vertically-integrated concentration are shown in the Equation 49
to Equation 52 for inner and outer bay.
For winter months without ice model, the residence time patterns in the inner bay are
similar for the case of with river and no river inflow with mean residence times of 67.7 days and
64.7 days respectively. The low residence times occur in the middle of the inner bay with the
value of about 60 days. For the case of dye release into the entire bay, the longer mean residence
times were found for the case of with and no river inflow of 131.8 days and 114.2 days
respectively. The lowest residence times were found in the areas in the northwest of the mouth of
the Saginaw Bay of about 75 days. This value is relatively higher the results from summer
season (5 to 50 days). The results from the winter simulation clearly show the important role of

149

the gyres near the mouth of the Saginaw Bay on residence time. The empirical equations
between the residence time and average vertically-integrated concentration are shown in the
Equation 53 to Equation 56 for inner and outer bay.
For winter months with ice model, the minimum of average vertically-integrated
concentration is 69.8% for the entire bay. For the inner bay, the dye concentration decrease only
5% through winter months. In this case, we still see the effect of the gyres near the mouth of the
Saginaw Bay where the residence times are about 130 days in the outer bay area.
The results of residence time from the Eulerian method above indicate large variation in
residence time between the cases includes river inflows and no river inflows. It implies the
significance role of the Saginaw River discharge to the residence time in the Saginaw Bay. A
very important question concerns the residence time in the Saginaw Bay under different Saginaw
River inflows, but this is still poorly understood. The variation of residence times are expected to
be consistent with how much the river discharge flows into the bay. In other words, there should
be a relationship between the residence time of the Bay and the river discharge.
The experimental investigations were undertaken to provide additional insights into the
relationship between river discharge and residence time. No difference is made in the dye
concentration modeling except the input of river discharge in which the river discharge was kept
as a constant through the whole simulation. Although, it may be not a real condition, it is very
seldom that there is a steady river flow rate for a long time but it is still useful for understanding
and interpreting the behavior of system to the river inflow.

150

Table 19: Statistical data for daily discharge in m3/s
Min

25th percentile

Median

Mean

75th percentile

Max

15.2

49.55

81.27

114.9

140.01

523.86

We used a statistical data of the USGS station on Saginaw River listed in Table 19. The
residence time estimated from dye concentration analysis versus fresh water input rate from the
Saginaw River shown in Figure 63. Exponential least-squares regression fit to the results are also
given for each case of simulation.
The correspondence between river flow rate and the mean residence time is high with R2
of about 0.97 (Figure 63). For the entire bay, the residence time is nearly 115 days at minimum
of flow rate from river of 15.2 m3/s, while under maximum of flow rate at 523.86 m3/s the
residence time to be 93 days. For the inner bay, the residence times are 85 days and 55 days,
respectively. There is not large difference in the residence time index at the long-term mean and
at the minimum flow rate for both entire and inner bay cases. While the residence time is
expected to be strongly influenced by the wind speed and wind direction, the strong connection
between the mean residence time and river discharge could be caused by other factors such as
mixing processes and dilution as well. It should be noted that there is no significant difference in
mean residence times for cases where the value of river discharge is less than the 75th percentile,
as illustrated in Figure 63. The declining mean residence time as a function of river discharge
follows an exponential curve with relatively high value of R2 of 0.97 as shown in Equation 57
and Equation 58.

151

For Inner Bay: (R2 = 0.976)

𝝉 = 𝟖𝟖. 𝟔𝟑𝒆−𝟎.𝟎𝟎𝟎𝟖𝟐𝟖𝟖𝑸

(57)

For Entire Bay: (R2 = 0.971)

𝝉 = 𝟏𝟏𝟕. 𝟕𝒆−𝟎.𝟎𝟎𝟎𝟒𝟏𝟗𝟐𝑸

(58)

Where Q is river discharge (unit: m3/s); 𝜏 is mean residence time (unit: in day).

152

Figure 63: Mean residence time as a function of river discharge for (a) Inner Bay and (b) Entire
Bay. The blue solid dots represent the model results while the red line represents an exponential
equation fit to the data. The residence times were calculated based on the dye concentration
model under summer meteorological conditions.

153

5.4

Conclusions

Residence time, an important environmental attribute of the bay is determined for both summer
and winter seasons. The data and interpretations in the study will bring about a broader
understanding of the residence time in Saginaw Bay. The numerical model result is used to
calculate the water mass exchange rates and residence time for Saginaw Bay. The influx rates
(m3/s) at two transects (1-1 and 2-2 in Figure 3a) across Saginaw Bay were used to estimate the
mean flushing times for the inner and outer Saginaw Bay. The mean flushing time for the inner
bay (~23.0 days) was approximately twice that for the entire bay (~9.9days) in summer season.
This difference is explained by the topographical features (sill) and the gyres that limit exchange
in the inner bay along with the large-scale circulation features that increase the exchange rate for
the outer bay. Comparing our results with earlier studies [Saylor et al., 1976] on the exchange
rate for Saginaw Bay, we found that the difference is explained by the change in volume of the
inner bay in the last four decades due to falling lake levels. There is a notable feature that the
mean flushing time for the inner bay (~ 43.2 days) was a significant difference from the
corresponding number for the entire bay (~15.6 days) for winter months.
In the summer, the mean e-folding flushing time was found to be 62 days (2 months) for the
inner bay and 115 days (3.7 months) for the entire bay. While in winter, the mean e-folding
flushing time was 64.7 days and 114.2 days for the inner and entire bay respectively. The
declining mean residence time as a function of river discharge was found as an exponential curve
with high value of R2 of 0.97 for both case of inner and entire bay.

154

Chapter 6
Conclusions
The aim of this study was to describe lake circulation, thermal structure, and exchange in the
Saginaw Bay - Lake Huron system under summer and winter conditions. The summer and winter
circulations were examined using an unstructured grid hydrodynamic model based on the
FVCOM. Several field datasets based on ADCP, thermistor chain and Lagrangian drifters as well
as satellite-based datasets for ice cover extent were used to test the numerical models. The
comparisons between the simulated and observed results show that the model was able to
reproduce the circulation patterns and temperature distribution in the lake very well. The model
was run for 3 constitutive years in summer and 2 years in winter seasons.
The hydrodynamic model was able to reproduce the anti-cyclonic (counterclockwise)
lake-wide circulation patterns that have been observed in earlier numerical and field studies on
Lake Huron for both summer and winter seasons. We found that mean current speeds in Lake
Huron were 3 cm/s and 1.14 cm/s during summer and winter periods respectively. The complex
topography coupled with stratification and rapid warming of the southern part of the lake results
in high inter-annual variability in surface temperature values during summer. Several potential
upwelling regions have been identified where topographically-driven upwelling brings cold
water from the deeper layers to the surface. One such region is in the strait connecting Lake
Huron with Georgian Bay. In addition to upwelling, our simulations show that in the surface
155

layer (epilimnion) mean flow during summer is from Lake Huron into Georgian Bay, which sets
up a return flow in the deeper layers (hypolimnion) from the bay into the lake. Since the
hypolimnion is 3-4 times thicker than the epilimnion in this part of the lake, surface velocities
are several times larger in magnitude than the velocities in the deeper layers.
In winter time, the simulations show that there is reverse flow at the mouth of Georgian Bay with
the outflow from Georgian Bay to Lake Huron at the surface layer and inflow at the bottom
layer. The results show that the ice model was able to reproduce the thermal evolution, presence
of ice cover, ice concentration, and ice thickness reasonably well. The effects of ice cover on the
circulation were also demonstrated where the inner bay was almost in a stagnation state during
February.
Based on the model results, the water mass exchange rates for Saginaw Bay were calculated.
The influx rates (m3/s) at two transects (1-1 and 2-2 in Figure 1) across Saginaw Bay were used
to estimate the mean flushing times for the inner and outer Saginaw Bay. The mean flushing time
for the inner bay (~23.0 days) was approximately twice that for the entire bay (~9.9days) in
summer time. This difference is explained by the topographical features (sill) and the gyres that
limit exchange in the inner bay and the large-scale circulation features that increase the exchange
rate for the outer bay. Comparing our results with earlier studies [Saylor et al., 1976] on the
exchange rate for Saginaw Bay, we found that the difference is explained by the change in
volume of the inner bay in the last four decades due to falling lake levels. The mean flushing
time for the inner and entire bay in winter time was found is 43.2 days and 15.6 days,
respectively. The mean e-folding flushing time was 62 days for summer and 64.7 days for winter
for the inner bay and 115 days for the summer and 114.2 days for the winter conditions for the
entire bay.

156

Results from the Lagrangian drifter study in Saginaw Bay showed the highly anisotropic and
time-variable nature of dispersion near the mouth of Saginaw River. The numerical model was
able to predict the drifter trajectories and dispersion coefficients reasonably well.

157

BIBLIOGRAPHY

158

BIBLIOGRAPHY
Andréfouët, S., Pages, J., & Tartinville, B. (2001). Water renewal time for classification of atoll
lagoons in the Tuamotu Archipelago (French Polynesia).Coral reefs, 20(4), 399-408.
Andutta, F. P., P. V. Ridd, and E. Wolanski (2012), The age and the flushing time of the Great
Barrier Reef waters. Continental Shelf Research.
Allender, J. H. (1975), Numerical simulation of circulation and advection-diffusion processes in
Saginaw Bay, Michigan. (Doctoral dissertation, University of Michigan).
Allender, J. H., & A. W. Green (1976), Free mode coupling of Saginaw Bay and Lake Huron.
Journal of Great Lakes Research, 2(1), 1-6.
Alm, E, W, Janice Burke, and Erin Hagan (2006), Persistence and Potential Growth of the Fecal
Indicator Bacteria, Escherichia coli, in Shoreline Sand at Lake Huron. J. Great Lakes
Res. 32:401–405.
Anderson, E. J., D. J. Schwab (2013), Predicting the oscillating bi-directional exchange flow in
the Straits of Mackinac. Journal of Great Lakes Research 39(4):9 pp.
Andradóttir, H. Ó., F. J. Rueda, J. Armengol, and R. Marcé(2012), Characterization of
residence time variability in a managed monomictic reservoir. Water Resources
Research, 48(11).
Annear, R. L., and S. A. Wells (2007), A comparison of five models for estimating clear‐sky
solar radiation. Water Resources Research, 43(10).
Assel, R. A. (2005). Classification of Annual Great Lakes Ice Cycles: Winters of 1973–
2002. Journal of climate, 18(22).
Assel, R. A., Robertson, D. M., Hoff, M., & Selgery, J. (1995). Climatic change implications of
long-term (1823–1994) ice records for the Laurentian Great Lakes. Ann. Glaciol, 21, 383386.
Assel, R., Cronk, K., & Norton, D. (2003). Recent trends in Laurentian Great Lakes ice
cover. Climatic Change, 57(1-2), 185-204.
Bai, X., J. Wang, C. Sellinger, A. Clites, and R. Assel (2012), Interannual variability of Great
Lakes ice cover and its relationship to NAO and ENSO, J. Geophys. Res., 117, C03002.

159

Bai, X., J. Wang, D.J. Schwab, Y. Yang, L. Luo, G. A. Leshkevich, and S. Liu
(2013),Modeling 1993-2008 climatology of seasonal general circulation and thermal
structure in the Great Lakes using FVCOM. Ocean Modeling, 65:40-63.
Beletsky, D., J. H. Saylor, and D. J. Schwab (1999), Mean circulation in the Great Lakes. J.
Great Lakes Res., 25,78–93.
Beletsky, D., & Schwab, D. J. (2001). Modeling circulation and thermal structure in Lake
Michigan: Annual cycle and interannual variability. Journal of Geophysical Research:
Oceans (1978–2012), 106(C9), 19745-19771.
Bennett Jr, T. J. (1982). A coupled atmosphere-sea ice model study of the role of sea ice in
climatic predictability. Journal of the Atmospheric Sciences, 39(7), 1456-1465.
Berliand, T. C. (1960). Method of climatological estimation of global radiation.Meteorol.
Gidrol, 6, 9-12.
Bierman, V. J, Jr. Jagjit Kaur, Joseph V. DePinto, Timothy J. Feist, and David W. Dilks,
Modeling the Role of Zebra Mussels in the Proliferation of Blue-green Algae in Saginaw
Bay, Lake Huron. J. Great Lakes Res. 31:32–55
Bolin, B., and , H. Rodhe (1973), A note on the concepts of age distribution and transit time in
natural reservoirs. Tellus, 25(1), 58-62.
Brown, J., J. MacMahan, and A. Reniers (2009), Surf zone diffusivity on a rip-channeled beach.
Journal of Geophysical Research 114.C11: C11015.
Boyce, F. M., M. A. Donelan, P. F. Hamblin, C. Murthy, T. Simons (1989), Thermal structure
and circulation in the Great Lakes. Atmosphere-Ocean 27, 607–642.
Bunker, A. F. (1976), Computations of surface energy flux and annual air-sea interaction cycles
of the North Atlantic ocean.
Burwell, D., M. Vincent, M. Luther, and B. Galperin (2000), Modeling residence times: Eulerian
versus Lagrangian, in Estuarine and Coastal Modeling, eds. M.L. Spaulding and H.L.
Butler, ASCE, Reston, VA, pp. 995-1009.
Camacho, R. A., and J. L. Martin (2012), Hydrodynamic Modeling of First-Order Transport
Timescales in the St. Louis Bay Estuary, Mississippi. Journal of Environmental
Engineering, 139(3), 317-331.
Chapra, S. C. (2008), Surface Water-Quality Modeling, Waveland Press, Long Grove, IL, USA.
Chen, C., H. Liu, R. C. Beardsley (2003), An unstructured, finite-volume, three-dimensional,
primitive equation ocean model: application to coastal ocean and estuaries. J. Atmos.
Oceanic Technol., 20, 159–186.

160

Chen, C. et al, An unstructured, finite-volume, three-dimensional, primitive equation ocean
model: FVCOM user manual (2013).
Danek, L. J., and J. H. Saylor (1977), Measurements of the summer currents in Saginaw Bay,
Michigan.Journal of Great Lakes Research, 3(1), 65-71.
Delhez, E. J. M., G. Lacroix, E. Deleersnijder (2004b), The age as a diagnostic of the dynamics
of marine ecosystem models. Ocean Dynamics 54, 221–231.
Dolan, D. M. (1975), Saginaw Bay residence time. Internal Working Paper, US EPA, Large
Lakes research Station, Grosse Ile, Michigan.
Fahnenstiel, G. L., Lang, G. A., Nalepa, T. F., & Johengen, T. H. (1995). Effects of Zebra
Mussel (Dreissena polymorpha) Colonization on Water Quality Parameters in Saginaw
Bay, Lake Huron. Journal of Great Lakes Research, 21(4), 435-448.
Fairall, C. W., Bradley, E. F., Rogers, D. P., Edson, J. B., & Young, G. S. (1996). Bulk
parameterization of air‐ sea fluxes for tropical ocean‐ global atmosphere coupled‐ ocean
atmosphere response experiment. Journal of Geophysical Research: Oceans (1978–
2012), 101(C2), 3747-3764.
Figueroa, D., & Moffat, C. (2000). On the influence of topography in the induction of coastal
upwelling along the Chilean coast. Geophysical Research Letters, 27(23), 3905-3908.
Fischer, H. B. (Ed.). (1979). Mixing in inland and coastal waters.
Fishman, D. B., Adlerstein, S. A., Vanderploeg, H. A., Fahnenstiel, G. L., & Scavia, D. (2009).
Causes of phytoplankton changes in Saginaw Bay, Lake Huron, during the zebra mussel
invasion. Journal of Great Lakes Research,35(4), 482-495.
Galperin, B., Kantha, L. H., Hassid, S., & Rosati, A. (1988). A quasi-equilibrium turbulent
energy model for geophysical flows. Journal of the Atmospheric Sciences, 45(1), 55-62.
Ge, Z., R. L. Whitman, M. B. Nevers, M. S. Phanikumar, and M. N. Byappanahalli (2012b),
Nearshore
hydrodynamics
as
loading
and
forcing
factors
forEscherichiacolicontamination at an embayed beach,Limnol. Oceanogr.,57(1), 362–
381, doi:10.4319/lo.2012.57.1.0362.
Geyer, W. R. (1997), Influence of wind on dynamics and flushing of shallow
estuaries. Estuarine, Coastal and Shelf Science, 44(6), 713-722.
Go´mez-Gesteira, M., M. deCastro, and R. Prego (2003), Dependence of the water residence
time in Ria of Pontevedra (NW Spain)on the seawater inflow and the river discharge.
Estuarine Coastal and Shelf Science 58, 567–573.

161

Harrington, M. W. (1894), Currents of the Great Lakes as deduced from the movements of bottle
papers during the seasons of 1892 and 1893. U.S. Weather Bureau, Washington, D.C.
Heath, R. T., Fahnenstiel, G. L., Gardner, W. S., Cavaletto, J. F., & Hwang, S. J. (1995).
Ecosystem-Level Effects of Zebra Mussels (Dreissena polymorpha): An Enclosure
Experiment in Saginaw Bay, Lake Huron. Journal of Great Lakes Research, 21(4), 501516.
Hecky, R. E., Smith, R. E., Barton, D. R., Guildford, S. J., Taylor, W. D., Charlton, M. N., &
Howell, T. (2004). The nearshore phosphorus shunt: a consequence of ecosystem
engineering by dreissenids in the Laurentian Great Lakes. Canadian Journal of Fisheries
and Aquatic Sciences, 61(7), 1285-1293.
Hsu, K., M. T. Stacey, and R. C. Holleman (2011), Exchange between an estuary and an
intertidal marsh and slough. Estuaries and Coasts, 1-13.
Hunke, E. C., & Lipscomb, W. H. (2010). CICE: the Los Alamos sea ice model, documentation
and software user’s manual, Version 4.1.
Kraines, S. B., Isobe, M., & Komiyama, H. (2001). Seasonal variations in the exchange of water
and water-borne particles at Majuro Atoll, the Republic of the Marshall Islands. Coral
Reefs, 20(4), 330-340.
Leppäranta, M. (2005). The drift of sea ice. Springer.
Liu, Z., H. Wang, X. Guo, Q. Wang, and H. Gao (2012), The age of Yellow River water in the
Bohai Sea. Journal of Geophysical Research: Oceans (1978–2012), 117(C11).
Jouon, A., P. Douillet, S. Ouillon, and P. Fraunié (2006), Calculations of hydrodynamic time
parameters in a semi-opened coastal zone using a 3D hydrodynamic model. Continental
Shelf Research, 26(12), 1395-1415.
Marvin C, S. Painter, and R. Rossmann (2004), Spatial and temporal patterns in mercury
contamination in sediments of the Laurentian Great Lakes. Environmental Research 95
(2004) 351–362.
McCormick, Miller, G.S., Murthy, C.R., Rao, Y.R., and Saylor, J.H. 2002. Tracking coastal flow
with surface drifters during the episodic events Great Lakes experiment. Verh. Internat.
Verein. Limnol. 28:365–369.
McCormick, M. J., Manley, T. O., Beletsky, D., Foley III, A. J., & Fahnenstiel, G. L. (2008).
Tracking the surface flow in Lake Champlain. Journal of Great Lakes Research, 34(4),
721-730.
Mellor, G. L., and T. Yamada (1982), Development of a turbulence closure model for
geophysical fluid problems, Rev. Geophys., 20, 851–875,1982.
162

Miller, G. S., & Saylor, J. H. (1981). Winter temperature structure in Lake Huron. Journal of
Great Lakes Research, 7(3), 201-206.
Monsen, N. E., J. E. Cloern, L. V. Lucas, and S. G. Monismith (2002), A comment on the use of
flushing time, residence time, and age as transport time scales, Limnology and
Oceanography, 47(5), pp. 1545-1553.
Nalepa, T. F., Fanslow, D. L., Pothoven, S. A., Foley III, A. J., & Lang, G. A. (2007). Long-term
trends in benthic macroinvertebrate populations in Lake Huron over the past four
decades. Journal of Great Lakes Research, 33(2), 421-436.
Nekouee, N. (2010), Dynamics and Numerical Modeling of River Plumes in Lakes, NOAA
Technical Memorandum GLERL-151, Great Lakes Environmental Research Lab, Ann
Arbor, MI, USA.
Niemelä, S., Räisänen, P., & Savijärvi, H. (2001). Comparison of surface radiative flux
parameterizations: Part I: Longwave radiation. Atmospheric Research, 58(1), 1-18.
Nixon, S. W., J. W. Ammerman, L. P. Atkinson, V. M. Berounsky, G. Billen, W. C. Boicourt,
W. R. Boynton, T. M. Church, D. M. Ditoro, R. Elmgrens, J. H. Garber, A. E. Giblin’o,
R. A. Jahnke, N. J. P. Owens, M. E. Q. Pilson, S. P. Seitzinger (1996), The fate of
nitrogen and phosphorus at the land-sea margin
of the North
Atlantic
Ocean,Biogeochemistry, 35: 141-180.
Ojo, T. O., J. S. Bonner, and C. Page (2006), Observations of shear-augmented diffusion
processes and evaluation of effective diffusivity from current measurements in Corpus
Christi Bay, Cont. Shelf Res.,26(6), 788–803.
Okubo, A. (1971), Oceanic diffusion diagram,Deep Sea Research,18, 789–902.
Parkinson, C. L., and Washington, W. M. (1979), A large‐scale numerical model of sea
ice. Journal of Geophysical Research: Oceans (1978–2012), 84(C1), 311-337.
Peeters, F., A. Wüest, G. Piepke, and D. M. Imboden (1996), Horizontal mixing in lakes,
J.Geophys.Res., 101, 18-18.
Phelps, J. J., J. A. Polton, A. J. Souza, and L. A. Robinson (2013), Hydrodynamic timescales in
a hyper-tidal region of freshwater influence.Continental Shelf Research.
Pickett, R. L. (1980). Observed and predicted Great Lakes winter circulations. Journal of
Physical Oceanography, 10(7), 1140-1145.
Plattner, S., D. M. Mason, G. A. Leshkevich, D. J. Schwab, and E. S. Rutherford (2006),
Classifying and forecasting coastal upwellings in Lake Michigan using satellite derived
temperature images and buoy data. Journal of Great Lakes Research, 32(1), 63-76.
163

Pope, S. B. (2000), Turbulent flows. Cambridge University Press.
Rao, D. B., & Murty, T. S. (1970). Calculation of the steady state wind-driven circulations in
Lake Ontario. Archiv für Meteorologie, Geophysik und Bioklimatologie, Serie A, 19(2),
195-210.
Robertson, D. M., Ragotzkie, R. A., & Magnuson, J. J. (1992). Lake ice records used to detect
historical and future climatic changes. Climatic Change, 21(4), 407-427.
Saylor, J. H., and L. J. Danek (1976), Wind-Driven Circulation of Saginaw Bay, Coastal
Engineering, pp. 3262-3275, ASCE.
Saylor J. H and G. S. Miller (1979), Lake Huron winter circulation, J.Geophys. Res.
Oceans,84:3237–3252.
Schwab D.J (1999), Lake Michigan Mass Balance Study: Hydrodynamic Modeling Project,
NOAA Technical Memorandum ERL GLERL-108.
Schertzer, W. M., R. A. Assel, D. Beletsky, T. E. Croley II, B. M. Lofgren, J. H. Saylor, D. J.
Schwab (2008), Lake Huron climatology, inter-lake exchange and mean circulation.
Aquatic Ecosystem Health & Management, 11(2):144–152.
Sheldon, J. E., and M. Alber (2006), The calculation of estuarine turnover times using
freshwater fraction and tidal prism models: A critical evaluation. Estuaries and
Coasts, 29(1), 133-146.
Sheng, J., Wright, D. G., Greatbatch, R. J., & Dietrich, D. E. (1998). CANDIE: A new version of
the DieCAST ocean circulation model. Journal of Atmospheric & Oceanic Technology, 15(6).
Sheng, J., Y. R. Rao (2006), Circulation and thermal structure in Lake Huron and GeorgianBay:
Application of a nested-grid hydrodynamic model, Continental Shelf Research, 26 1496–
1518.
Simons, T. J., Murthy, C. R., & Campbell, J. E. (1985). Winter circulation in lake
Ontario. Journal of Great Lakes Research, 11(4), 423-433.
Sloss, P. W., and J. H. Saylor (1976), Large-scale current measurements in Lake Huron,
J.Geophys. Res.81:3069–3078.
Smagorinsky, J. (1963). General circulation experiments with the primitive equations: I. The
basic experiment*. Monthly weather review, 91(3), 99-164.
Spydell, M., F. Feddersen, R. T. Guza, and W. E. Schmidt (2007), Observing surf-zone
dispersion with drifters, Journal of Physical Oceanography, 37(12), 2920-2939.

164

Spydell, M. S., F. Feddersen, and R. T. Guza (2009),Observations of drifter dispersion in the
surfzone: The effect of sheared alongshore currents, J. Geophys. Res. Oceans,114:
C07028.
Stocker, R., and J. Imberger(2003), Horizontal Transport and Dispersion in the Surface Layer of
a Medium-Sized Lake,Limnology and Oceanography, 48(3), 971-982.
Takeoka, H. (1984), Fundamental concepts of exchange and transport time scales in a coastal
sea. Continental Shelf Research 3, 311–326.
Taylor, G. I. (1921), Diffusion by continuous movements, Proc. London Math. Soc.,20,196–212.
Thupaki, P., M. S. Phanikumar, D. Beletsky, D. J. Schwab, M. B. Nevers, and R. L. Whitman
(2010). Budget analysis of Escherichia coli at a southern Lake Michigan
beach. Environmental science & technology, 44(3), 1010-1016.
Thupaki, P., M. S. Phanikumar and R. L. Whitman (2013), Solute dispersion in the coastal
boundary layer ofsouthern Lake Michigan,J. Geophys. Res. Oceans, 118(3),pp. 16061617, doi:10.1029/2012JC008286.
Wan, Y., C. Qiu, P. Doering, M. Ashton, D. Sun, and T. Coley (2013). Modeling residence
time with a three-dimensional hydrodynamic model: Linkage with chlorophyll a in a
subtropical estuary. Ecological Modelling, 268, 93-102.

165