NOAA Technical Memorandum ERL GLERL-108

LAKE MICHIGAN MASS BALANCE STUDY:

HYDRODYNAMIC MODELING PROJECT

David J. Schwab

NOAA, Great Lakes Environmental Research Laboratory, Ann Arbor, MI

Dmitry Beletsky

Dept. of Naval Architecture and Marine Engineering, Univ. of Michigan, Ann Arbor, MI

Cooperative Institute for Limnology and Ecosystems Research/NOAA, Great Lakes Environmental

 Research Laboratory and Univ. of Michigan, Ann Arbor, MI

Great Lakes Environmental Research Laboratory

Ann Arbor, Michigan

October 1998

NAT OI NAL OCEANIC AND ATMOSPHERIC ADMINISTRAT OI N

U.S.DEPARTMENT OF COMMERCE UNITED STATES

DEPARTMENT OF COMMERCE

William Daley

Secretary

NATIONAL OCEANIC AND

ATMOSPHERIC ADMINISTRATION

D. James Baker

Under Secretary for Oceans

and Atmosphere/Administrator

Environmental Research

Laboratories

James L. Rasmussen

Director
2

NOTICE

Mention of a commercial company or product does not constitute an endorsement

by the NOAA Environmental Research Laboratories. Use of information from

this publication concerning proprietary products or the tests of such products for

publicity or advertising purposes is not authorized. This is GLERL Contribution

No. 1111.

ACKNOWLEDGEMENTS

Support provided by the USEPA MED-Duluth, Large Lakes Research Station -

Grosse Ile, MI (USEPA IAG Number DW-13936717).We would like to thank

Rick Towler and Glenn Muhr for programming support, Cathy Darnell for editing, Jim Saylor, Jerry Miller, Mike McCormick, and Nathan Hawley for providing current, temperature, and sediment data, and Ray Assel and George Leshkevich for providing ice data. The National Biological Service provided temperature survey data in 1994-95. The National Weather Service and the Illinois State

Water Survey provided meteorological observations. The National Data Buoy

Center provided meteorological and oceanographic observations at buoys. The

National Ocean Service provided water level data. The EPA Great Lakes National

Program Office (Chicago) provided CTD survey data in 1994-95. Municipal

water intake plants provided water temperature and turbidity observations at

Benton Harbor, MI, Bridgman, MI, Holland, MI, Ludington, MI, Muskegon, MI,

Muskegon Heights, MI, St. Joseph, MI, Grand haven, MI, Cudahy, WI, Green

Bay, WI, Kenosha, WI, Manitowoc, WI, Marinett, WI, Sheboygan, WI, Oak

Creek, WI, East Chicago, IN, Borman Park, IN, Ogden Dunes, IN, Michigan City,

IN, Whiting, IN, and Chicago, IL.
3

CONTENTS

PAGE

Abstract ........................................................................................................................................................... 6

1. Introduction ...................................................................................................................................................... 6

1.1 LMMBS Modeling Framework......................................................................................................... 6

1.2 Hydrodynamic Modeling Program Overview................................................................................... 7

2. Model Description ............................................................................................................................................ 8

2.1 Wave Model ....................................................................................................................................... 8

2.2 Hydrodynamic Model...................................................................................................................... 11

2.3 Bathymetric Grids ........................................................................................................................... 14

3. Forcing Functions ........................................................................................................................................... 14

4. Model Results and Comparison with Observations ....................................................................................... 19

4.1 1982 - 1983 Period .......................................................................................................................... 19

4.1.1 Temperature ........................................................................................................................... 22

4.1.2 Currents .................................................................................................................................. 22

4.1.3 Water Levels .......................................................................................................................... 27

4.1.4 Waves ..................................................................................................................................... 29

4.2 1994 - 1995 Period .......................................................................................................................... 29

4.2.1 Temperature ........................................................................................................................... 29

4.2.2 Currents .................................................................................................................................. 32

4.2.3 Water Levels .......................................................................................................................... 32

4.2.4 Waves ..................................................................................................................................... 32

5. Aggregation of Hydrodynamic Model Output for Level II Water Quality Model ......................................... 32

6. Interpretation/Discussion ................................................................................................................................ 38

7. Conclusions .................................................................................................................................................... 40

8. References ...................................................................................................................................................... 41

9. Products ......................................................................................................................................................... 44

9.1 Publications ..................................................................................................................................... 44

9.2 Presentations .................................................................................................................................... 44

9.3 Final Report CD: Graphics, Programs, Observations, and Selected Model Results

(including averages and aggregation) .............................................................................................. 46

9.3.1 Model Codes .......................................................................................................................... 47

9.3.2 Shell Scripts ........................................................................................................................... 49

9.3.3 Programs for Input/Output Analysis (IDL) ............................................................................ 50

9.4 Model Input and Output on CD - Binary (XDR) and ASCII .......................................................... 51
4

Figures

PAGE

Figure 1. LMMBS modeling framework ......................................................................................................... 7

Figure 2. Model grid.......................................................................................................................................... 9

Figure 3. Observation network for 1982-1983 ................................................................................................ 15

Figure 4. Observation network for 1994-1995 ................................................................................................ 16

Figure 5. 30-day average wind field, April 1982 ............................................................................................ 20

Figure 6. 30-day average total surface heat flux, April 1982 .......................................................................... 21

Figure 7. Simulated mean temperature profile for 1982-1983 ........................................................................ 23

Figure 8. Time-series of simulated water temperature versus observed at 45007 and

45002 for 1982-1983 ....................................................................................................................... 24

Figure 9. Stream-function for summer and winter circulation for 1982-1983 ................................................ 26

Figure 10. Progressive vector diagrams of simulated currents versus observed at stations

1,2,3, and 23 for 1982-1983 ............................................................................................................. 28

Figure 11. Time series of simulated water levels versus observed for October-November 1983.

Green Bay, WI (a); Milwaukee, WI (b); Calumet Harbor, IL (c) .................................................... 30

Figure 12. 30 day average wind wave height, April 1982................................................................................. 31

Figure 13. Time series of simulated wave height versus observed for April-May 1982,

NDBC buoy 45002 .......................................................................................................................... 33

Figure 14. Time-series of simulated water temperature versus observed at 45007 for 1994-1995 .................. 34

Figure 15. Time-series of simulated surface water temperature versus observed at 45002 and

45010 for 1994-1995 ....................................................................................................................... 34

Figure 16. Time series of simulated water temperature versus observed at Muskegon municipal

water intake for May-June 1994 ...................................................................................................... 35

Figure 17. Temporal evolution of simulated versus observed temperature profiles, station 18M .................... 36

Figure 18. Simulated versus observed temperature, transects offshore of Ludington, MI

(upper panel) and Holland, MI (lower panel) .................................................................................. 36

Figure 19. Simulated versus observed short-wave radiation for May-June 1994 at ISWS

station at Beaver Island .................................................................................................................... 36

Figure 20. Stream-function for summer and winter circulation for 1994-1995. ............................................... 37

Figure 21. Level II model grid .......................................................................................................................... 39
5

Tables

PAGE

Table 1. Vertical distribution of σ-levels in the hydrodynamic model ..............................................................13

Table 2 Meteorological station locations ..........................................................................................................17

Table 3. Mooring locations, 1982-83 .................................................................................................................17

Table 4. Mooring locations, 1994-95 .................................................................................................................17

Table 5 Municipal water intake locations ........................................................................................................17

Table 6. Water level gauge locations ................................................................................................................18

Table 7. 1982-83 Hydrodynamic model evaluations for surface temperature at NDBC

buoys (45002 and 45007) and subsurface temperature at GLERL current meter

moorings (28 instruments) ...................................................................................................................23

Table 8. Statistical comparison of 1982-83 observed and computed short term water

level fluctuations (cm) .........................................................................................................................29

Table 9. Statistical comparison of 1982-83 observed and computed wave heights (m) ....................................31

Table 10. Statistical comparison of 1982-83 observed and computed wave periods (s) .....................................31

Table 11. 1994-95 Hydrodynamic model evaluations for surface temperature at NDBC

buoys (45002, 45007, and 45010) and subsurface temperature at GLERL

current meter moorings (10 instruments) .............................................................................................34

Table 12. Statistical comparison of 1994-95 observed and computed short-term water level

fluctuations (cm) ..................................................................................................................................35

Table 13. Statistical comparison of 1994-95 observed and computed wave heights (m) ....................................37

Table 14. Statistical comparison of 1994-95 observed and computed wave periods (s) .....................................37
6

LAKE MICHIGAN MASS BALANCE STUDY:

HYDRODYNAMIC MODELING PROJECT

David J. Schwab

NOAA, Great Lakes Environmental Research Laboratory, Ann Arbor, MI

Dmitry Beletsky

Dept. of Naval Architecture and Marine Engineering, Univ. of Michigan, Ann Arbor, MI

Cooperative Institute for Limnology and Ecosystems Research/NOAA, Great Lakes Environmental

Research Laboratory and the Univ. of Michigan, Ann Arbor, MI

ABSTRACT. This report describes the hydrodynamic modeling framework for the U.S.

EPA Lake Michigan Mass Balance Study. It consists of a three-dimensional lake circulation model, surface flux model for atmospheric input, and a wind wave model. These

models provide a description of the physical environment for sediment resuspension and

transport models, as well as eutrophication and toxic contaminant models. The models

are validated using an extensive array of long-term measurements of temperature, currents, water levels, and wind waves in Lake Michigan during the 1982-83 calibration

period and the 1994-95 Mass Balance Study period.

1. INTRODUCTION

In order to develop improved strategies for management and control of toxic chemicals in the coastal environment, the USEPA Great Lakes National Program Office has undertaken a program to measure and model the

transport, fate, and bioaccumulation of four particular chemicals in Lake Michigan: PCB’s (industrial compounds

once widely used in a variety of products, banned since 1982), trans-nonachlor (a chlorinated hydrocarbon

originally registered as a pesticide in 1948, banned by EPA in 1988), atrazine (the most widely used herbicide in

U.S. corn and sorghum production), and mercury (a toxic element which occurs both naturally and anthropogenically). The project consists of a large monitoring program, including tributary inputs, atmospheric inputs,

sedimentary processes, and biological processes, and also a modeling program.

1.1 LMMBS Modeling Framework

The modeling framework for LMMBS (Figure 1) consists of a series of linked submodels which can be divided

into three groups: computational transport models, mass balance models, and bioaccumulation models. The

computational transport models include a hydrodynamic model to estimate the three dimensional velocity and

temperature fields, a wave model, and a sediment and particulate transport model. A surface flux model is used to

convert raw meteorological data into gridded estimates of heat and momentum flux suitable for use as forcing

functions in the hydrodynamic, wave, and sediment transport models.

The eutrophication, sorbent dynamics, and transport and fate models will be used in conjunction with the results

of the computational transport models and the measured constituent loadings to simulate the seasonal cycle of

primary production in the lake as well as the transport, intermedia exchange, phase distribution, and biogeochemical transformation of the target chemicals through the water column and the sediments. In the original project

design, these models were intended to be applied on a collapsed (Level II) computational grid with considerably

less resolution than the computational transport models (Level III). The drawbacks of this approach are well

known, particularly the loss of spatial resolution as well as the difficulty of providing realistic estimates of
7

dispersion coefficients for chemical constituents in the collapsed grid calculations. In order to minimize these

problems, the original plan was modified to include a provision for applying computational transport and mass

balance models with the same computational grid resolution (5 or 10 km horizontal resolution with 20 vertical

levels). With this approach, results from the transport models may not be as realistic as results with higher grid

resolution, but the mass balance computations should be considerably more realistic than results using the collapsed grid approach.

The data collection phase of LMMBS provided information on atmospheric and tributary loadings, as well as in

situ chemical concentrations for model validation for 1994-1995. The hydrodynamic models were applied for this

entire period, and the results will be used to provide three dimensional transport fields for the mass balance

models. Mass balance models predict chemical concentrations in water and sediment, as well as the bioavailability of toxic chemicals. The food web bioaccumulation model will then be used to estimate the accumulation of

chemical constituents in various elements of the food web, ranging from benthos and zooplankton to forage fish

and the top predators (lake trout and coho salmon).

1.2 Hydrodynamic Modeling Program Overview

The main goal of hydrodynamic modeling program is to calculate three-dimensional fields of currents and temperature, and wind wave characteristics for the period of study. The wind wave model used in this task is the

GLERL/Donelan parametric wind wave model developed by Schwab et al. (1984a, 1984b). The wave model has

been successfully applied to Lake Erie (Schwab et al., 1984a) and Lake Michigan (Liu et al., 1984), as well as the

Baltic Sea and several other lakes and embayments around the world. Schwab et al. (1984a) compared wave

model results to wave measurements from an instrumented tower in Lake Erie and found root mean square

differences on the order of 0.2 m for wave height and 1 sec for wave period. Liu et al. (1984) showed a high

correlation between model results and lake-wide synoptic wave height measurements from an airborne laser

altimeter in Lake Michigan. The GLERL/Donelan model is also used operationally by the NWS (Johnson et al.,

1992) and has proven to be highly accurate when wind forecasts are accurate.

While there has been significant progress in hydrodynamic modeling simulation of wind waves in the Great

Lakes, long-term circulation modeling efforts were rare in the past (Schwab, 1992). For example, since the

pioneering works of Simons (1974, 1980) created the basis for numerical studies of circulation and thermal

structure in the Great Lakes, Lake Michigan experienced only two long-term modeling exercises: Allender and

Advective/

Dispersive

Transport

Food Web

Bioaccumulation

Model

resuspension

flux

Sediment

Transport

Model

(SEDZL)

G.L. Wave

Model

Hydrodynamic

Model (POM/3D)

Eutrophication

Model

Transport &

Fate Model

meteorology

data

constituent

loadings Solids/Sorbent

Dynamics

Collapsed Grid

Transport Model

phytoplankton

benthos

chemical

exposure & availability Atrazine

MercuryPCBs/TNC

trophic

dynamics

Computational

Transport Mass Balance Bioaccumulation

Lake Michigan Modeling Framework

Figure 1. LMMBS modeling framework.
8

Saylor (1979) simulated three-dimensional circulation and thermal structure for an 8-month period, and Schwab

(1983) studied circulation with a two-dimensional barotropic model also for an 8-month period.

Nowadays, with increases in computer power, seasonal variations in thermal structure and circulation can be more

easily studied using three-dimensional hydrodynamic models. The numerical circulation model used in this task is

a three-dimensional ocean circulation model developed at NOAA’s Geophysical Fluid Dynamics Lab at Princeton

University for coastal ocean applications by Blumberg and Mellor (1987) and subsequently adapted for Great

Lakes use at GLERL (Schwab and Bedford, 1994, O’Connor and Schwab, 1994). The Princeton Ocean Model

has been used extensively for coastal and estuarine applications, including the Middle Atlantic Bight, the

South Atlantic Bight, The California Shelf, the Santa Barbara Channel, and New York Harbor. The Great Lakes

version is used operationally in the Great Lakes Forecasting System (Bedford and Schwab, 1990, Schwab and

Bedford, 1994) for Lake Erie. Extensive validation tests with observed currents, water level fluctuations, and

surface temperature distributions have been carried out in the development of the Great Lakes Forecasting System

(Kuan et al., 1994). The physical parameters predicted by the model are the three-dimensional velocity distributions, the temperature field, turbulent kinetic energy, turbulent mixing length, and the free surface water level.

In this study, the Princeton model was applied to Lake Michigan for two 2-year periods: 1982-1983, and 1994-

1995. The first 2-year period was chosen for the model calibration because of an extensive set of observational

data (Figure 2) including surface temperature observations at two NDBC weather buoys, and current and temperature observations during June 1982 - July 1983 at several depths from 15 subsurface moorings (Gottlieb et al.,

1989).

There is no ice modeling component in the present version of the model, which can be a problem for the annual

cycle modeling in general, because ice cover can cause significant changes in winter circulation patterns in a large

lake (Campbell et al., 1987). The Great Lakes are usually at least partially covered with ice from December to

April. Maximum ice extent is normally observed in late February, when ice typically covers 45% of Lake Michigan (Assel et al., 1983). Fortunately enough, it was not the case for the chosen periods of study: the 1982-83 and

1994-95 winters were among the warmest winters of the century and therefore practically ice-free (Assel et al.,

1985).

2. MODEL DESCRIPTION

2.1 Wave Model

The wave model is a numerical finite-difference solution to the two-dimensional wave momentum conservation

equation. The wave energy spectrum is parameterized at each point on a rectilinear computational grid in terms of

total wave energy, peak energy period, and predominant wave direction. Momentum is transferred from the wind

to the waves using Donelan’s (1979) formulation, which depends on the difference between the phase velocity of

the waves and the local wind velocity. The momentum balance equation is

∂

∂

M v M

t

+ •∇ = g wτ (eq. 1)

where ∇ is the horizontal gradient operator, and M and vg are the total momentum vector and the corresponding

group velocity vector, and τw is that part of the momentum input from the wind that produces net wave momentum growth. Assuming equipartition of kinetic and potential wave energy in the wave field, the momentum vector

can be expressed as

M = 



 



 = ( )

( )





 

 ∫∫ 

M

M g

F f

c f d df x

y

ρw

θ θ

θ

θ , cos

sin (eq. 2)
9

Lake Michigan Hydrodynamic Model

 5 km Computational Grid

Figure 2. Model grid.
10

where ρw is the water density, c(f) the phase speed and F(f,θ) the two-dimensional frequency spectrum of wave

energy as a function of frequency f and direction θ. Assuming further that the wave energy has a cosine-squared

angular dependence about the mean angle independent of frequency and that there is no energy for

θθ π − > 0 / 2, we use

Ff Ef , cos θ

π ( ) = ( ) ( ) θ θ − 2 2

0 (eq. 3)

where E(f) is the one-dimensional spectral energy density. Since we are concerned with total momentum, we

define the mean-square surface elevation:

σ θθ

ρ

2 = ( ) = ∫∫ F f dfd

c

g

p

w

, M

(eq. 4)

where MMM = + ( ) x y

2 2 1 2/

 and cp is the deep-water phase speed of the peak of the spectrum, cp =g/(2πf

p).

Applying deep-water linear theory, Schwab et al. (1984a) have shown that

v M

vM vM

v M

c x x

x yyx

y y

p ,

cos

cos sin

sin





















=

+

+





























M

4

1

2

1

2

2

0

0 0

2

0

θ

θ θ

θ

(eq. 5)

We calculate τw as

τw a = −− 0 028 0 83 0 83 Df p p . .( .) ρ U cU c (eq. 6)

where cp is a vector with magnitude cp in the direction of the wave momentum vector, i.e. c MM p p = c / ,

ρa

 is the air density, U the 10 m wind vector, Df the form drag coefficient given by

Df = [ . / ln( / )] 0 4 50 2 σ (eq. 7)

with σ in meters. The factor 0.028 is the empirical fraction of the wind stress that is retained by the waves.

To solve momentum balance equation we still need a relationship between wave momentum and wave height. We

use the following empirical relation derived from JONSWAP relations (Hasselman et al., 1973) and linking σ 2

with peak energy frequency f

p and wind speed U:

σ2 6

10 3 4

2 = × 6 23 10 



 



 −

−

.

/

f U

g

U

g

p

(eq. 8)

The model is thus semi-empirical and parametric. A simple numerical integration scheme can then be applied to

momentum balance equation. Forward time differences are used to calculate the momentum components at the

center of the grid squares, and a combination of upwind and centered differences are used to evaluate the momentum advection terms at the edges of the grid squares. Model output at each grid point consists of significant wave

height (defined by H1/3 = 4σ ), peak-energy wave period and average wave direction.
11

The output from the 5 km wave model (wave height, wave period, and wave direction) will be used in the sediment transport model to estimate bottom shear stress in each computational grid square at each time step.

2.2 Hydrodynamic Model

The Princeton model is a nonlinear, fully three-dimensional, primitive equation, finite difference model that

solves the equations of fluid dynamics. The model is hydrostatic and Boussinesq so that density variations are

neglected except where they are multiplied by gravity in the buoyancy force. The model uses time-dependent

wind stress and heat flux forcing at the surface, zero heat flux at the bottom, free-slip lateral boundary conditions,

and quadratic bottom friction. The drag coefficient in the bottom friction formulation is spatially variable. It is

calculated based on the assumption of logarithmic bottom boundary layer using constant bottom roughness of 1

cm.

In order to simplify the discussion of model physics, we present the dynamical equations in Cartesian coordinates.

The velocity components (u,v,w) are in the (x,y,z) directions. The mass continuity equation is

∇• + = V ∂

∂

w

z

0 (eq. 9)

where V= (u,v) is the horizontal velocity. The horizontal momentum equations are

∂

∂

∂

∂ ρ

∂

∂

∂

∂

∂

∂

∂

∂

∂

∂

∂

∂

∂

∂

u

t

u

w

z

fv p

x x

A u

x y

A u

y z

K u

z

+ •∇ + − =− + MMM





 



 +





 



 + 



 



 V 1

0

(eq. 10)

∂

∂

∂

∂ ρ

∂

∂

∂

∂

∂

∂

∂

∂

∂

∂

∂

∂

∂

∂

v

t

v w

v

z

fu p

y x

A v

x y

A v

y z

K v

z

+ •∇ + + =− + MMM





 



 +





 



 + 



 



 V 1

0

(eq. 11)

where ρ is density, p is pressure, f is the Coriolis parameter, and AM and KM are the horizontal and vertical momentum eddy viscosities, respectively. The Eulerian derivatives at a point are the result of the horizontal and

vertical velocity advections, horizontal pressure gradient force, Coriolis force, and horizontal and vertical momentum diffusion. The model is hydrostatic and Boussinesq, so that density variations are neglected except where

they are multiplied by gravity in the buoyancy force. The internal energy conservation equation with no sources or

sinks of heat, gives an expression for the temperature T :

∂

∂

∂

∂

∂

∂

∂

∂

∂

∂

∂

∂

∂

∂

∂

∂

T

t

T w T

z x

A T

x y

A T

y z

K T

z

+ •∇ + = HHH





 



 +





 



 + 



 



 V (eq. 12)

where AH and KH are the horizontal and vertical thermal diffusivities, respectively.

The finite difference form of the Princeton model is described by Blumberg and Mellor (1987). Horizontal

diffusion is calculated with a Smagorinsky eddy parameterization (with a multiplier C=0.1) to give a greater

mixing coefficient near strong horizontal gradients.

A Cxy u x v x u y v y m = ∆ ∆ [ ] ( ) ∂∂ ∂∂ ∂∂ ∂∂ / / // / + + ( ) + ( ) / 2 22 1 22 (eq. 13)

Horizontal momentum diffusion will be assumed to be equal to horizontal thermal diffusion, as is common

practice in hydrodynamic models where the primary horizontal mixing process is eddy diffusion (Blumberg,
12

1986). The equation of state (Mellor, 1991) calculates the density as a function of temperature, salinity, and

pressure. For applications to the Great Lakes, the salinity is set to a constant value of 0.2 ppt.

The terrain following vertical coordinate system (sigma-coordinate) replaces the vertical coordinate, z, with a

normalized vertical coordinate, σ ηη =− + ( ) /( ) z d , where d is the local depth and η is the surface elevation.

The advantage of this system is that in the transformed coordinate system, the bottom corresponds to a uniform

value of the vertical coordinate (σ = -1), thus simplifying the governing transport and continuity equations. The

disadvantage is that an extra term is introduced in horizontal gradient terms that can lead to artificial vertical

diffusion of heat and momentum, particularly in areas of large topographic gradients. The Princeton model

calculates the bottom friction coefficient, Cz based on the water depth d, and bottom roughness length z0. Using

law of the wall, near-bottom velocity at the height zb can be described by:

u uk zz b b = ( ) * / ln / ( ) 0 (eq. 14)

where k=0.40 is the von Karman constant, and u* is the bed-shear velocity, which meets the following relation:

u Cu * z b

2 2 = (eq. 15)

 The friction coefficient can be found from (14)-(15):

 Cz

 = [k/ln(zb/z0)]2 (eq. 16)

 To provide sufficient bottom friction in deep water, in the POM

Cz = max{0.0025, 0.16/[ln(zb/z0)]2} (eq. 17)

 In shallow water zb can be less than z0 as it was found in Lake Erie simulations by O’Connor and Schwab (1994).

Therefore, the model uses zb = 2z0 in such cases. In earlier versions of POM, bottom roughness length was depthdependent: z0=0.01(1+100/d). In the latest version of POM, z0 is assigned a constant value of 0.01 m. For the

vertically integrated mode of POM, the drag coefficient is taken as a constant 0.0025.

Although the current version of the model can incorporate a curvilinear, coastline-following coordinate system,

this feature is not used in the Great Lakes version. We felt that the additional complications of a curvilinear

coordinate system in (1) the interpolation and analysis of model results, and (2) the effect of variable grid size on

numerical truncation errors were not justified by the potential for increased accuracy in the hydrodynamic model.

The equations are written in flux form, and the finite differencing is done on an Arakawa-C grid using a control

volume formalism. The finite differencing scheme is second order and centered in space and time (leapfrog).

A stratified body of water such as a lake has two types of motions, the barotropic (density independent) mode and

the baroclinic (density dependent) mode. The Princeton model uses a mode splitting technique that solves the

barotropic mode for the free surface and vertically averaged horizontal currents, and the baroclinic mode for the

fully three-dimensional temperature, turbulence, and current structure. This necessitates specifying both a barotropic and baroclinic mode time step in accordance with the Courant-Friedrich-Lewy computational stability criterion.

The model includes the Mellor and Yamada (1982) level 2.5 turbulence closure parameterization. The vertical

mixing coefficients for momentum KM and heat KH are calculated from the variables describing the flow regime.

The turbulence field is described by prognostic equations for the turbulent kinetic energy q2/2 and turbulent length

scale l,
13

∂

∂

∂

∂

∂

∂

∂

∂

∂

∂

∂

∂ ρ

∂ρ

∂

q

t

q wq

z z

K q

z

K u

z

v

z

g K

z

q

B l qMH

2

2

2 2 2 2

0

3

1

2 2 2 + •∇ + = 



 



 + 



 



 + 



 



 



 



V  + − (eq. 18)

∂

∂

∂

∂

∂

∂σ

∂

∂

∂

∂

∂

∂ ρ

∂ρ

∂

q l

t

qlw q l

z

K q l

z

E lK u

z

v

z

E l g K

z

q

B qMH W

2

2

2 2

1

2 2

1

0

3

1

( ) + •∇( ) + ( ) =  ( )















+ 



 



 + 



 



 



 



V  + − ƒ

(eq. 19)

The first term on the right in each equation arises from the vertical diffusion of turbulent kinetic energy, and Kq is

the diffusivity for this variable. The next two terms arise from the production of turbulent kinetic energy from

shear and from buoyancy, respectively. The last terms with q3 represent the dissipation of turbulent energy. The B1

and E1 are empirical constants, g is the acceleration of gravity, andWƒ is the empirical wall proximity function

which approaches unity away from the surface. The problem is closed by expressing the vertical mixing coefficients for momentum, heat, and turbulent kinetic energy in the form

K lq S M m = (eq. 20 )

K lq S H = h (eq. 21)

Kq=lqSq (eq. 22)

where Sm, Sh, and Sq are analytically derived algebraic stability functions. Details of this procedure are given by

Blumberg and Mellor (1987).

The hydrodynamic model of Lake Michigan has 20 vertical levels and a uniform horizontal grid size of 5 km.

Vertical levels were spaced more closely in the upper 30 m of water and near the bottom to better resolve both the

seasonal thermocline and bottom boundary layer (Table 1).

The external mode time step is 30 s, the internal mode time step is 5 min. In the Princeton Ocean model, there is

an option to allow a portion of the short-wave radiation to penetrate the upper part of the water column according

to one of the Jerlov’s (1976) five optical categories (I, IA, IB, II, III). We have found that allowing for short-wave

radiation penetration improved the simulated vertical thermal structure in the multi-annual hydrodynamic model

runs over runs with no short-wave radiation penetration. However, more work needs to be done to determine the

optimal parameterization of short-wave radiation penetration for the Great Lakes. We use the default optical

category (IPSWR=2, NTP=2) in the model. We also set a minimum water temperature in the model (0oC) to

prevent negative water temperatures in winter due to the lack of an ice model. An alternative procedure for

handling water temperatures below 0oC in the model would be to

allow water temperature to become negative, and then eliminate

surface momentum flux in grid squares where surface water

temperature was less than 0oC (A. Blumberg, personal communication). We felt that because ice cover in the 1982-83 and 1994-95

winters was so low (maximum ice cover reached only 17% in

1982-83 and 20% in 1994-95 versus 45% during a normal year)

the first procedure would be adequate and would have negligible

effect on conservation of heat in the model.

Table 1. Vertical

σ-levels used in

hydrodynamic

model.

Level σ Level σ

1 .000 11 -.2270

2 -.0227 12 -.2724

3 -.0454 13 -.3405

4 -.0681 14 -.4313

5 -.0908 15 -.5448

6 -.1135 16 -.6810

7 -.1362 17 -.7945

8 -.1589 18 -.8853

9 -.1816 19 -.9534

10 -.2043 20 -1
14

There is no open boundary in the model, which means that we neglect the influence of tributaries and outflow

through the Straits of Mackinac on large-scale lake circulation. In Lake Michigan the hydraulic flow is at least an

order of magnitude smaller than typical wind-driven and density-driven currents. However, hydraulic flow may be

more important for longer time scale simulations.

The output from the lake circulation model can be used to provide estimates of horizontal advection and bottom

shear stress for the sediment resuspension and transport model, as well as the water quality models. In addition,

the turbulence closure scheme in the circulation model can provide estimates of physical dispersion coefficients

for water quality and toxics models.

2.3 Bathymetric Grids

Three bathymetric grids were developed for LMMBS: 2 km, 5 km, and 10 km. The bathymetry was derived from

the 2 km gridded bathymetric data compiled by Schwab and Sellers (1980). The 5 km gridded depths are slightly

smoothed by adjusting the depths to ensure that the relative depth change between adjacent grid squares was less

than 0.5 while still preserving the volume of the original grid. The 2 km grid was used only in the idealized

upwelling -Kelvin wave studies (Beletsky et al.,1997). The 5 km grid (Figure 2) was used through the course of

our study as the main grid, for which all model evaluations were done. The 10 km grid was developed for prototype testing of the coupled hydrodynamic- biogeochemical model.

3. FORCING FUNCTIONS

In order to calculate heat and momentum flux fields over the water surface for the lake circulation, sediment

transport, and wave models, it is necessary to estimate wind, temperature, dew point, and cloud cover fields at

model grid points. Meteorological data were obtained from the NWS weather stations and buoys as well as

additional marine observations from U.S. Coast Guard stations and ships of opportunity in Lake Michigan (Table

2). These data are routinely collected and quality-controlled at the Cleveland Weather Service Forecast Office. In

addition, data from several meteorological stations in the LMMBS air sampling network around Lake Michigan

were used. The marine observation network is shown in Figures 3 and 4 (in 1982-83 simulations only NWS

weather stations and buoys data were available). These observations form the basis for generating gridded

overwater wind, temperature, dew point, and cloud cover fields. Figures 3 and 4 also show observation networks

of currents and temperature at moorings (Tables 3 and 4), temperature at municipal water intakes (Table 5), and

water levels (Table 6) that were used for model evaluation.

Three main steps are required to develop overwater fields from the marine observation data base: (1) height

adjustments, (2) overland/overlake adjustment, and (3) interpolation. First, measurements must be adjusted to a

common anemometer height. Ship observations are usually obtained at considerably higher distances above the

water surface than buoy measurements. Measurements are adjusted to a common 10 m height above the water

surface using profile methods developed by Schwab (1978) and described more thoroughly by Liu and Schwab

(1987). The wind and temperature profiles are represented as

u z u k zz w m ( ) ( / )[ln( / ) ] = − * 0 Ψ (eq. 23)

Tz T T k zz a h ( ) ( / )[ln( / ) ] =+ − 0 0 * Ψ (eq. 24)

where u is wind speed; z is the vertical coordinate; u* is friction velocity, k=0.4 is the von Karman constant; Ta is

air temperature; T0 = Ta(0) is the surface temperature, T* is the scaling temperature; and ψm and ψh are functions

of dimensionless stability height given by Long and Shaffer (1975). In conjunction with the Charnock (1955)

relation for overwater surface roughness, z0=0.045 u*

2 /g, these equations can be solved iteratively to obtain z0 and

u*, yielding profiles for uw(z) and Ta(z).
15

Figure 3. Observation network for 1982-1983.
16

Figure 4. Observation network for 1994-1995.
17

Station Location Lat Lon

45002 NDBC Data Buoy 02 45.30 86.40

45007 NDBC Data Buoy 07 42.70 87.00

45010 NDBC Data Buoy 10 43.00 87.80

SGNW3 Sheboygan Brkwtr WI 43.75 87.69

0Y2 Sturgeon Bay CG WI 44.80 87.31

17C Ludington CG MI 43.95 86.40

20C St. Joseph CG MI 42.12 86.48

Y09 Charlevoix CG MI 45.32 85.27

GSRW3 Gills Rock WI 45.30 86.98

BEH Benton Harbor MI 42.13 86.43

CIU Sault Ste. Marie MI 46.25 84.48

ESC Escanaba MI 45.75 87.10

GRB Green Bay WI 44.48 88.13

ISWBVI Beaver Island MI 45.73 85.54

ISWIDN Indiana Dunes IN 41.63 87.09

ISWIIT Illinois Inst Tch IL 41.83 87.62

ISWKEN Kenosha WI 42.51 87.81

ISWMSK Muskegon MI 43.23 86.34

ISWMTW Manitowoc WI 44.08 87.68

ISWSBR Sleeping Br Dunes MI 44.76 86.06

ISWSHV South Haven MI 42.46 86.17

MKE Milwaukee WI 42.95 87.90

MKG Muskegon MI 43.17 86.23

ORD OHare Field Chi IL 41.98 87.90

SBN South Bend IN 41.70 86.32

TVC Traverse City MI 44.73 85.58

Y62 Sault Ste. Marie MI 46.50 84.30

Table 3. Mooring locations, 1982-83.

Table 2. Meteorological station locations.

Table 4. Mooring locations, 1994-95.

Table 5. Municipal water intake locations.

Mooring Lat Lon

1 42.69 86.42

2 42.05 87.04

3 42.69 87.52

4 43.32 87.77

5 43.89 87.41

6 43.75 87.35

7 43.61 87.29

8 43.49 87.21

9 43.37 87.17

10 43.33 87.14

11 43.46 86.73

12 43.86 86.86

13 44.51 86.83

23 42.71 87.06

33 45.30 86.30

Station Location Lat Lon

101 Benton Harbor, MI 42.13 86.49

108 Bridgman, MI 41.94 86.60

121 Grand Haven, MI 43.06 86.27

109 Holland, MI 42.80 86.23

102 Ludington, MI 43.96 86.47

103 Muskegon, MI 43.21 86.35

110 Muskegon Heights, MI 43.18 86.32

111 St. Joseph, MI 42.10 86.50

112 Cudahy, WI 42.96 87.82

113 Green Bay, WI 44.49 87.47

114 Kenosha, WI 42.59 87.79

115 Manitowoc, WI 44.08 87.62

116 Marinette, WI 45.10 87.60

104 Sheboygan, WI 43.67 87.67

117 Oak Creek, WI 42.87 87.83

105 East Chicago, IN 41.66 87.40

106 Borman Park, IN 41.64 87.34

118 Ogden Dunes, IN 41.63 87.20

119 Michigan City, IN 41.74 86.90

107 Whiting, IN 41.46 87.49

120 Chicago, IL 41.90 87.59

Mooring Lat Lon

 1 43.01 87.80

 2 43.00 87.86

 3 43.08 87.84

 4 43.01 87.57

 5 43.04 87.03

 6 43.05 86.65

 7 42.72 87.06

 8 42.29 86.64

10 42.00 86.69

19 43.05 86.66

23 43.23 86.42

24 43.19 86.39

27 43.19 86.43

33 42.93 86.29

34 42.93 86.26

35 42.93 86.33
18

The second problem in dealing with the combination of

overland and overwater measurements is that overland

wind speeds generally underestimate overwater values

because of the marked transition from higher aerodynamic roughness over land to much lower aerodynamic

roughness over water. This transition can be very abrupt

so that wind speeds reported at coastal stations are often

not representative of conditions only a few kilometers

offshore. Schwab and Morton (1984) found that wind

speeds from overland stations could be adjusted by

empirical methods to obtain fair agreement with overlake

wind speeds measured from an array of meteorological

buoys in Lake Erie. For meteorological stations that are

more representative of overland than overwater conditions, namely airports and other “surface stations” in the

marine observation network (Figures 3-4), we apply the empirical overland-overlake wind speed adjustment from

Resio and Vincent (1977). For wind speed, these take the form

uw = ul

F1(ul)F2(∆T) (eq. 25)

where F1(ul

) =1.2 + 1.85/ ul (m/s), and FT TT T u 2 w

1 3 () / / , 1 1 1920 /

∆ ∆∆ ∆ = − ( ){ } − ( ) is the over water wind

speed; ul is over land wind speed; and ∆T = Ta-Tw (Degrees C).

For heat flux calculations, we also need to know overwater humidity to calculate latent heat flux. Dew point

observations are only available from land stations. Phillips and Irbe (1978) used a simple empirical formula to

estimate overwater dew point temperature from overland values:

Tdw = Tdl-c1(Tdl-Tw) (eq. 26)

where Tdw is the dewpoint temperature over water; Tdl is the dewpoint temperature over land; and c1 is on the

order of 0.35 for neutral stability. Air temperature reports from overland stations are adjusted with similar empirical formula:

Ta = 0.4Tal+0.6Tw (eq. 27)

where Ta is the air temperature over water; Tal is the air temperature over land; and Tw is lake-averaged surface

water temperature.

Finally, in order to interpolate meteorological data observed at irregular points in time and space to a regular grid

so that it can be used for input into numerical wave, sediment transport, and circulation models, some type of

objective analysis technique must be used. The complexity of the analysis technique should be compatible with

the complexity of the observed data, i.e., if observations from only a few stations are available, a best-fit linear

variation of wind components in space might be an appropriate method. If more observations can be incorporated

into the analysis, spatial weighting techniques can be used. For LMMBS we used the nearest-neighbor technique,

with the addition of a spatial smoothing step (with a specified smoothing radius). The nearest neighbor technique

assigns the value of the nearest measurement station to each point in the regular grid, similar to the Thiessen

polygon weighting scheme (Thiessen, 1911). The spatial smoothing step replaces each value on the regular grid

with the average of all grid points within the specified smoothing radius. In the nearest neighbor technique, we

also consider observations from up to three hours before the interpolation time to three hours after the interpolation time. In the nearest-neighbor distance calculations, the distance from a grid point to these observation points

is increased by the product of the time difference multiplied by a scaling speed. The interpolation scaling speed is

taken as 10 km/hr. Interpolation smoothing distance is 30 km. We found that the nearest-neighbor technique

provided results comparable to results from the inverse power law or negative exponential weighing functions.

Table 6. Water level gauge locations.

Station Location Lat Lon

7023 Ludington 43.95 86.45

7031 Holland 42.47 86.20

7044 Calumet Harbor 41.73 87.53

7057 Milwaukee 43.00 87.88

7068 Kewanee 44.47 87.50

7072 Sturgeon Bay 44.80 87.32

7079 Green Bay 44.53 88.00

7096 Port Inland 45.97 85.87
19

Figure 5 shows the average of the first 30 days of hourly wind fields for 1982 generated by this procedure.

After we have produced hourly gridded overwater fields for wind, dew point, air temperature, and cloud cover, the

momentum flux and heat flux can be calculated at each horizontal grid square in the three-dimensional lake

circulation model at each model time step. To calculate momentum flux, the profile theory described above for

anemometer height adjustment is used at each grid square at each time step to estimate surface stress, using the

surface water temperature from the circulation model. This procedure provides estimates of bulk aerodynamic

transfer coefficients for momentum and heat. Surface heat flux, H, is calculated as

H = Hsr + Hs + Hl + Hlr (eq. 28)

where Hsr is short-wave radiation from the sun, Hs

 is sensible heat transfer, Hl is latent heat transfer, and Hlr is

long wave radiation. The heat flux procedure follows the methods described by McCormick and Meadows (1988)

for mixed-layer modeling in the Great Lakes. Hsr is calculated based on latitude and longitude of the grid square,

time of day, day of year, and cloud cover (CL).

Hsr = Hcs F3(CL) (eq. 29)

where Hcs is a clear sky value and F3 is a cubic function of cloud cover that ranges from 1.0 for clear sky to 0.36

for total cloud cover. Hs

 and Hl are calculated using the bulk aerodynamic transfer formulas:

Hs

 = ChCpρauw ∆T (eq. 30)

Hl = Cdqlρauw (ha-hw) (eq. 31)

where Ch is bulk heat coefficient, Cp is specific heat of air at constant pressure, ∆T is water-air temperature

difference, Cd is drag coefficient, ql is latent heat of vaporization, ha is specific humidity of air, and hw is specific

humidity at water surface. Hlr is calculated as a function of Ta, T, and cloud cover according to Wyrtki (1965).

Figure 6 shows the average of the first 30 days of hourly net surface heat flux fields for 1982.

McCormick and Meadows (1988) showed that this procedure works quite well for modeling mixed layer depth in

the Great Lakes. The gross thermal structure generated in the three dimensional model using these heat flux fields

is similar to the profile that would be obtained from a one dimensional model. However, there is considerable

horizontal variability in the three dimensional temperature field due mainly to wind forcing.

4. MODEL RESULTS AND COMPARISON WITH OBSERVATIONS

4.1 1982-83 Period

To initialize the model, we used surface temperature observations at two buoys (45007 and 45002) located in the

southern and northern parts of the lake respectively. The model run starts on March 31, 1982. Vertical temperature

gradients are very small because of convection during that time of the year when the water temperature is less

than the temperature of maximum density (4oC). Therefore, we set vertical temperature gradients to zero, but

retained horizontal gradients. The horizontal temperature distribution at the beginning of the run depends on depth

(d) and latitude (converted to j- coordinate) according to the formula:

T = max( 0.0161 d - 0.0167 j, 0) (where 1 < j < 102, 5 km grid) (eq. 32)

This formula was chosen so that initial temperature distribution matched observed surface water temperature at

NDBC buoys 45002 and 45007 and also exhibited decreasing temperatures toward shallow water

consistent with existing observations for that time of year.
20

Figure 5. 30-day average wind field, April 1982.
21

Figure 6. 30-day average total surface heat flux, April 1982.
22

The typical range of water temperature in Lake Michigan for that time of year is from 0 to 3oC. Extensive experiments with different initial conditions for the 1994-95 period revealed that statistically, model results were

insensitive to specified initial temperature conditions (necessarily simplified due to insufficient observations) if

the mean lake temperature difference between different model runs was within approximately 0.5-1oC. Based on

that fact, it would probably be fair to say that for the 1982-83 period, any chosen realistic three-dimensional

temperature field yielding mean lake temperature about 1.5oC will provide results close to the ones that used

initialization scheme in (32). For both the 1982-83 and 1994-95 Lake Michigan simulations, the initial velocity

field is set to zero. Because of the strongly wind-driven character of the circulation in the lake, the effect of the

initial condition for velocity on calculated currents disappears within the first few weeks of the simulation or even

first few days if there is a strong wind event. The model run ends on November 20, 1983, after 600 days of

integration.

4.1.1 Temperature

The most distinctive feature of the physical limnology of the Great Lakes is a pronounced annual thermal cycle

(Boyce et al., 1989). By the end of fall, the lakes usually become vertically well-mixed from top to bottom at

temperatures near or below the temperature of maximum density for freshwater, about 4oC. Further cooling

during winter can lead to inverse stratification and ice cover. Springtime warming tends to heat and stratify

shallower areas first, leaving a pool of cold water (less than 4oC and vertically well-mixed because of convection)

in the deeper parts of the lake. In spring, stratified and homogeneous areas of the lake are separated by a sharp

thermal front, commonly known as the thermal bar. (See surface temperature plots on the Final Report CD.)

Depending on meteorological conditions and depth of the lake, the thermal bar may last for a period of from 1 to 3

months. Stratification eventually covers the entire lake, and a well-developed thermocline generally persists

throughout the summer. In the fall, decreased heating and stronger vertical mixing tend to deepen the thermocline

until the water column is again mixed from top to bottom. When the nearshore surface temperature falls below the

temperature of maximum density, the fall thermal bar starts its propagation from the shoreline toward the deeper

parts of the lake. Thermal gradients are much smaller during this period than during the springtime thermal bar.

The model was able to reproduce all of the basic features of thermal structure of Lake Michigan during the 600

day period of study: spring thermal bar, full stratification, deepening of the thermocline during the fall cooling,

and finally an overturn in the late fall (Figure 7). Observed temperatures from surface buoys and subsurface

moorings were compared to model output (Figure 8). The comparison is quite good for the horizontal distribution

and time evolution of the surface and bottom temperature, but it is worse in the thermocline area where internal

waves are also much less pronounced than in observations. We think that because the model tends to generate

excessive vertical diffusion, the modeled thermocline is too diffuse and hence temperature fluctuations are

decreased. On the other hand, the simulation of the surface temperature is much more accurate, which shows

correct calculation of heat fluxes near the surface. We should also note that the accuracy of surface temperature

predictions is similar in the first and second year of simulations (we do not have subsurface observations for the

second year summer) which we attribute to the rapid adjustment of the surface temperature field to the boundary

conditions. Statistics of temperature field validation are presented in Table 7. To insure comparability with the

1994-95 period, only temperature and current observations longer than 300 days were used. In the table, RMSD is

the root mean square difference (error) between observed and computed temperatures, Max Error is the maximum

temperature difference, Avg. is the arithmetic mean, and CC is the correlation coefficient. In Table 7, the correlation coefficient gives the usual statistical indication of the strength of the linear relationship between computed

and observed variables. We have found it to be a good measure of the agreement of the timing of events in computed and observed time series. Values greater than 0.5 usually indicate significant correlation in timing.

4.1.2 Currents

Wind-driven transport is a dominant feature of circulation in the lakes. As shown by Bennett (1974), Csanady

(1982), and others, the response of an enclosed basin with a sloping bottom to a uniform wind stress consists of

longshore, downwind currents in shallow water, and a net upwind return flow in deeper water. The streamlines of
23

Table 7. 1982-83 Hydrodynamic model evaluations for surface temperature at NDBC buoys (45002 and

45007) and subsurface temperature at GLERL current meter moorings (28 instruments).

RMSD Max Error Avg. (obs) Avg. (comp) CC

Surface 1.2 6.6 12.1 12.1 0.99

Subsurface

epilimnion

hypolimnion

2.5

0.7

10.6

3.3

7.1

4.2

6.4

4.3

0.87

0.78

Figure 7. Simulated mean temperature profile for 1982-1983.

Simulated Average Temperature of L. Michigan

100 200 300 400 500 600

JULIAN DAYS, 1982-83

100

80

60

40

20

0

Depth (meters)

100 200 300 400 500 600

JULIAN DAYS, 1982-83

100

80

60

40

20

0

Depth (meters)

Simulated Average Temperature of L. Michigan 2 2 4

4

4 4

4 4

6

6

6

8

6

8 8

10

8

10

10

12

1416

1214 16

1820
24

Surface Temperature at 45007

100 200 300 400 500 600

0

5

10

15

20

25

30

Temperature (C)

Surface Temperature at 45002

100 200 300 400 500 600

0

5

10

15

20

25

30

Temperature (C)

 Temperature at 15 m

100 200 300 400 500 600

0

5

10

15

20

25

30

Temperature (C)

 Temperature at 15 m

100 200 300 400 500 600

0

5

10

15

20

25

30

Temperature (C)

 Temperature at 50 m

100 200 300 400 500 600

0

5

10

15

20

25

30

Temperature (C)

 Temperature at 50 m

100 200 300 400 500 600

0

5

10

15

20

25

30

Temperature (C)

 Temperature at 148 m

100 200 300 400 500 600

JULIAN DAYS, 1982-83

0

5

10

15

20

25

30

Temperature (C)

 Temperature at 154 m

100 200 300 400 500 600

JULIAN DAYS, 1982-83

0

5

10

15

20

25

30

Temperature (C)

Figure 8. Time-series of simulated water temperature versus observed at 45007 and 45002 for

1982-1983. Red line is observation, blue is model simulation.
25

the flow field form two counter-rotating closed gyres, a cyclonic gyre to the right of the wind and an anticyclonic

gyre to the left (in the northern hemisphere). As the wind relaxes, the two-cell streamline pattern rotates cyclonically within the basin, with a characteristic period of about 4 days corresponding to the lowest mode topographic

wave of the basin (Saylor et al., 1980). Numerical models approximating actual lake geometry have proven to be

effective in explaining observed short-term circulation patterns in lakes (Simons, 1980, Schwab, 1992). The

results of these modeling exercises show that the actual bathymetry of each of the Great Lakes tends to act as a

combination of bowl-shaped sub-basins, each of which tends to support its own two-gyre circulation pattern.

Besides bathymetry and geometry, two other important factors tend to modify the simple two-gyre lake circulation

model described above, namely nonuniform wind forcing and stratification. Thus, during the stratified period,

longshore currents frequently form a single cyclonic gyre circulation pattern driven by onshore-offshore density

gradients. The effect of horizontal variability in the wind field enters through the curl of the wind stress field (Rao

and Murthy, 1970). Any vorticity in the forcing field is manifest as a tendency of the resulting circulation pattern

toward a single gyre streamline pattern, with the sense of rotation corresponding to the sense of rotation of the

wind stress curl. Because of the size of the lakes, and their considerable heat capacity, it is not uncommon to see

lake-induced mesoscale circulation systems superimposed on the regional meteorological flow, a meso-high in the

summer (Lyons, 1971) and a meso-low in the winter (Petterssen and Calabrese, 1959). There are also indications

that nonlinear interactions of topographic waves can contribute to the mean single gyre cyclonic circulation

(Simons, 1985).

Recent long-term current observations in Lake Michigan suggested a cyclonic large-scale circulation pattern, with

cyclonic circulation within each subbasin, and anticyclonic circulations in ridge areas (Gottlieb et al., 1989). Our

model results coincide with their conclusions. To study seasonal changes in circulation patterns, we averaged

model results over two 6-month periods: from May to October (summer period), and from November to April

(winter period). The selected averaging periods roughly correspond to the stratified and non-stratified periods.

Circulation is more organized and more cyclonic in winter than in summer, which is in agreement with Gottlieb et

al. (1989) and earlier findings of Saylor et al. (1980). One can also notice the presence of the computational noise

in the summer circulation by the small scale wiggles in the streamfunction contour lines. Because winter circulation is stronger than summer circulation, annual circulation looks very similar to winter circulation (Figure 9).

We also compared progressive vector diagrams of simulated and observed currents (Figure 10). The largest

currents occur in the fall and winter, when temperature gradients are lowest, but winds are strongest. Nearshore

currents appear to be much stronger than offshore currents, in agreement with existing conceptual models and

observations (Csanady, 1982). Still, even driven by the nearshore currents similar to ones that are depicted in

Figure 10, it takes a passive particle about a year to complete a round trip of Lake Michigan, which would be

about 1000 km. The point to point comparison of currents was most successful in the southern basin, which is

characterized by a relatively smooth bathymetry. It was more successful in fall-winter months than in summer,

most probably because the horizontal resolution of the model is not adequate for proper simulation of baroclinic

processes with horizontal length scales comparable to the Rossby deformation radius (which is around 5 km for

summer months). In addition, model resolution was too coarse to describe precisely the dynamics in the areas of

strong depth gradients, even in the fall and winter when lake dynamics are essentially barotropic.

The statistics of simulated currents evaluation is presented in the form of the Fourier norms (rms difference). The

Fourier norm of time-series of observed current vectors vo and computed vc is defined as

vv v v oc o c ,

/

= − 









 =

∑1 2 1 2

M t t

M t

∆

∆

(eq. 33)

We use a normalized Fourier norm: Fn oc o = vv vo , /, . In the case of perfect prediction Fn=0. In the case

 0 < Fn < 1, model predictions are better than no prediction at all (zero currents). Using Fn not only allows us to
26

Figure 9. Stream-function for summer and winter circulation for 1982-1983.

-0.6

-0.6

-0.6

-0.4

-0.4

-0.2

-0.2

0.0

0.0

0.0

0.0

May-October 1982

-0.6

-0.6

-0.4

-0.4

-0.2

-0.2

0.0

0.0

0.0

0.0

0.0

0.0

November 1982-April 1983

-0.6

-0.6

-0.6

-0.4 -0.4

-0.4

-0.2

-0.2

0.0

0.0

0.0

May-October 1983

Normalized Stream Function

Blue - counterclockwise rotation

Yellow - clockwise rotation

Figure 10. Progressive vector

diagrams of simulated

currents versus observed at

stations 1,2,3, and 23 for

1982-1983. Black line is

observations, green is model

simulation.
27

use one number for characterization of model skills in predicting vector quantities, but also to compare our model

results more objectively with previous model results. For example, in one of the earlier modeling exercises,

Allender (1977) obtained 1.00 < Fn < 1.11. Later, Schwab calculated 0.79 < Fn < 1.01. In the 1982-83 simulations 0.70 < Fn < 0.98. The normalized Fourier norm can also be thought of as the relative percentage of variance

in the observed currents that is unexplained by the calculated currents. So, for the 1982-83 simulations, the hourly

computed currents account for 2-30% of the total variance observed in the hourly current meter records.

4.1.3 Water levels

Lake water levels are routinely measured on an hourly basis by NOAA National Ocean Service at 8 sites around

Lake Michigan as shown in Figure 3. Because of its considerable depth, Lake Michigan exhibits considerably

smaller short term water level fluctuations than Lake Erie, where the wind-induced changes in water level often

exceed 1 m (Schwab, 1978). In Lake Michigan, typical wind-induced water level fluctuations are only in the

range of 10-20 cm, which makes it more difficult to use them for model validation purposes.

The free surface water level fluctuation is one of the dynamic variables simulated by the hydrodynamic model.

The hydrodynamic model simulation does not include seasonal changes in lakewide mean water level due to

precipitation and evaporation. Therefore, in order to compare measured water levels to output from the hydrodynamic model, it is necessary to remove the seasonal fluctuation of mean lake level from the individual water level

records. Hourly water levels were assembled and an hourly average was computed for all stations except Green

Bay and Sturgeon Bay. Green Bay was excluded because it is not on the main part of the lake and is highly

sensitive to local influences, and Sturgeon Bay was excluded because it is on a canal connecting Green Bay with

the lake proper and again may not be representative of conditions in the main part of the lake. The average time

series was smoothed with a 12 hour running mean filter and the subtracted from the observed hourly values at all

stations. The resulting water level fluctuations represent the combined effects of wind, atmospheric pressure

gradients, and local bathymetry. Since the hydrodynamic model forcing does not include atmospheric pressure

gradients, and the hydrodynamic model resolution of 5 km is too coarse to resolve many local effects (harbor

resonance, edge waves, etc.), the comparison of modeled water levels to observations is expected to be less

satisfactory than for Lake Erie, where the wind-induced part of the water level fluctuation is much larger than the

part due to atmospheric pressure gradients, or local effects.

A sample of the observed and computed water level fluctuation for Calumet Harbor, Milwaukee, and Green Bay

for the period JD 265-325, 1983 is shown in Figures 11a-c. Calumet Harbor and Green Bay show the highest

fluctuations (up to 50 cm), while Milwaukee does not exhibit as high an amplitude. At all stations there is a

significant amount of high frequency (> 0.5/hr) fluctuation with amplitude on the order of 5 cm apparent in the

observed water level record which is not seen in the computed water levels. We believe this difference is due

mainly to local water level fluctuations in the harbors and coastal areas where the water level gauges are located

which cannot be represented on the 5 km hydrodynamic grid.

The statistical comparison between observed and computed water level fluctuations is presented in Table 8. In the

table, Npts is the number of water level observations available at that station (out of a possible 14400), Avg is the

arithmetic mean, RMSV is the root mean square value, Diff. is the difference in the means, RMSD is the root

mean square difference, and CC is the correlation coefficient. It can be seen from Table 8 that the RMSV of the

observations is significantly greater than the RMSV of the computed water level fluctuations at all stations except

Green Bay. This is probably due to the high frequency fluctuations which are seen in the observations, but not as

much in the model. The RMSD between model results and observations is also reflective of this difference.

Correlation coefficients are highest at Calumet Harbor and Port Inland which are located at the south and north

ends of the main lake respectively, and which exhibit the largest amplitude wind-induced water level fluctuations.

4.1.4 Waves
28

(a)

(b)

(c)

Figure 11. Time series of simulated water levels versus observed for October-November 1983.

Green Bay, WI (a); Milwaukee, WI (b); Calumet Harbor, IL (c).
29

Table 8. Statistical comparison of 1982-83 observed and computed short-term water level fluctuations (cm).

NOAA routinely operates two weather buoys in Lake Michigan during the ice-free season which provide hourly

observations of surface meteorology, as well as wave height and wave period. Wave data from the buoys (45002

and 45007 in Figure 13) were assembled and compared to calculated wave height and wave period. As shown in

Table 9, calculated wave heights compared well to observed wave heights, with root mean square differences of

0.34 and 0.30 m and correlation coefficients of 0.89 and 0.92 at buoys 45002 and 45007 respectively. The model

has a tendency to underestimate the wave period as shown in Table 10. Observed mean wave periods were 0.71

and 0.64 s greater than calculated periods at 45002 and 45007. Monthly averaged simulated wave height for April,

1982 is presented in Figure 12. Figure 13 shows an example of the model output wave height, period, and direction, and the buoy observations at 45002 for the period April to May 1982.

4.2 1994-95 period

The first three months of 1994 were more challenging for modeling because of the extensive ice cover during the

winter of 1993-94. Maximum ice cover reached 78% (Assel et al., 1996). In the absence of an ice model, we

chose to keep the water temperature steady at 2oC until April when the ice is gone. Therefore, the model was run

in a diagnostic mode (constant temperature) until April 1. We experimented with different initial conditions based

on the AVHRR-derived lake surface temperature and vertical temperature profile at one mooring in southern Lake

Michigan. Those experiments suggested that qualitatively and quantitatively neither method of initialization

provided better results in 1994-95. Therefore, we switched the model to a prognostic mode on April 1 keeping the

same uniform temperature distribution.

4.2.1 Temperature

As in the 1982-83 case, the second year winter was very mild and practically ice-free. Therefore, water temperature was close to the temperature of maximum density (4oC). The result of that lake preconditioning was a much

shorter thermal bar period in the second spring. Observed temperatures from surface buoys and subsurface

Gauge Npts Avg

(Obs.)

RMSV

(Obs.)

Avg

(Comp.)

RMSV

(Comp.)

 Diff. RMSD CC

 7023 14371 -0.28 3.45 0.59 0.95 -0.87 3.32 0.39

 7031 13730 0.33 2.73 0.67 1.64 -0.35 2.39 0.51

 7044 14253 0.63 7.03 0.58 4.67 0.05 5.50 0.62

 7057 14233 -0.27 3.46 0.12 1.64 -0.39 3.22 0.39

 7068 13123 -0.07 3.42 -0.25 0.97 0.18 3.39 0.17

 7072 14366 -0.13 2.97 -0.17 1.01 0.04 2.92 0.22

 7079 14368 0.27 9.72 0.46 12.20 -0.18 12.30 0.39

 7096 14283 -0.32 5.70 0.33 2.73 -0.65 4.61 0.61
30

Figure 12. 30-day average wind wave height, April 1982.
31

Table 10. Statistical comparison of 1982-83 observed and computed wave periods (s).

Table 9. Statistical comparison of 1982-83 observed and computed wave heights (m).

Buoy Npts Avg

(Obs.)

Avg

(Comp.)

Diff. RMSD CC

45002 10356 0.85 0.69 0.17 0.34 0.89

45007 10217 0.82 0.70 0.12 0.30 0.92

Buoy Npts Avg

(Obs.)

Avg

(Comp.)

Diff. RMS CC

45002 10356 4.07 3.36 0.71 1.31 0.67

45007 10217 4.00 3.37 0.64 1.39 0.67

Figure 13. Time series of simulated wave height versus observed for April-May

1982, NDBC buoy 45002.
32

moorings were compared to model output. Again, the comparison is better for the surface and bottom temperature,

but it is worse in the thermocline area where internal waves are also much less pronounced than in observations

(Figure 14). In 1994-95 an additional set of surface temperature data was available from the NDBC buoy 45010

located in the coastal zone off Milwaukee, WI (Figure 15). Statistics of temperature field validation are presented

in Table 11. As with the 1982-83 period, only time-series longer than 300 days were used. The accuracy of surface

temperature predictions was similar in the first and the second year of simulations, which is in agreement with the

1982-83 findings. On the other hand, the accuracy of subsurface temperature predictions does improve in the

second year, especially in the thermocline area, as was shown from data collected from the only mooring with

1994-95 summer observations (Figure 14). Another piece of information on nearshore thermal structure was

obtained from 23 municipal water intakes around Lake Michigan. These data were used mostly for qualitative

assessment (Figure 16) because unlike typical in-situ measurements, water temperature is measured at the municipal water plants, not in the lake itself. There were also seven large-scale GLNPO temperature surveys of Lake

Michigan during the 1994-95 period that allowed us to compare observed and simulated evolution of temperature

profiles at some 20 locations (one example is shown in Figure 17) and several nearshore transects done by USGS

(two examples are shown in Figure 18). Note that model bathymetry is somewhat smoothed as discussed in 2.3.

Finally, solar radiation measurements at several ISWS stations were available for comparison with calculated

solar radiation. That comparison showed good agreement between measured and calculated short-wave radiation

(one example is shown in Figure 19).

4.2.2 Currents

Currents measurements in 1994-95 were not as comprehensive as in 1982-83. Statistics of simulated currents

evaluation showed that 0.80 < Fn < 1.08, which is slightly worse than in the 1982-83 case. Calculated circulation

patterns were somewhat different from the 1982-83 case. Winter circulation was predominantly cyclonic (Figure

20). Summer circulations exhibit cyclonic circulations around the deepest portions of northern and southern basin,

but also an anticyclonic circulation in the mid-lake ridge area. Analysis of wind field vorticity showed that it has

its maximum cyclonic vorticity during the winter, and minimum during summer. It is likely that as wind field

cyclonic vorticity approaches zero, more and more anticyclonic circulations emerge. More research is needed to

support this hypothesis.

4.2.3 Water Levels

In 1994-95, the same eight water level stations were available for comparison with model output. The results of

the comparisons are shown in Table 12. The same conclusions apply as for the 1982-83 water levels, although the

correlation coefficients between observed and computed water level fluctuations are somewhat lower in 1994-95.

4.2.4 Waves

In 1994-95, a third NDBC buoy was deployed near Milwaukee in Lake Michigan (Figure 4) as part of a limitedterm nearshore hydrodynamics experiment funded by NOAA. Wave height and period from this buoy as well as

the other two NDBC buoys are compared to calculated wave height and period in Tables 13 and 14. The root

mean square differences in wave height were somewhat lower in 1994-95 than in 1982-83, but correlation coefficients were not any higher. There is still a tendency for the model to underestimate wave period, although not as

much at 45010.

5. AGGREGATION OF HYDRODYNAMIC MODEL OUTPUT FOR LEVEL II

WATER QUALITY MODELING

In order to use computed hydrodynamic transport in LMMBS Level II water quality models, we developed

aggregation protocols for linking the Level II and Level III model grids. The Level II (aggregated) grid would be

created as an overlay on the Level III (hydrodynamic) grid. The Level II grid can be specified either by indicating

the boundaries between segments or by describing the area of the segments in terms of Level III grid. We have
33

Figure 14. Time-series of simulated water temperature versus

observed at 45007 for 1994-1995.

Red line is observation, blue line is

model simulation.

Surface Temperature at 45007

100 200 300 400 500 600

0

5

10

15

20

25

30

Temperature (C)

 Temperature at 12 m

100 200 300 400 500 600

0

5

10

15

20

25

30

Temperature (C)

 Temperature at 42 m

100 200 300 400 500 600

0

5

10

15

20

25

30

Temperature (C)

 Temperature at 102 m

100 200 300 400 500 600

JULIAN DAYS, 1994-95

0

5

10

15

20

25

30

Temperature (C)
34

Table 11. 1994-95 Hydrodynamic model evaluations for surface temperature at NDBC buoys (45002,

45007, and 45010) and subsurface temperature at GLERL current meter moorings (10 instruments).

Figure 16. Time series of simulated water temperature versus observed at

Muskegon municipal water intake for May-June 1994.

Figure 15. Time-series of simulated surface water temperature versus observed at 45002 and

45010 for 1994-1995. Red line is observation, blue is model simulation.

Surface Temperature at 45002

100 200 300 400 500 600

JULIAN DAYS, 1994-95

0

5

10

15

20

25

30

Temperature (C)

Surface Temperature at 45010

100 200 300 400 500 600

JULIAN DAYS, 1994-95

0

5

10

15

20

25

30

Temperature (C)

RMSD Max Error Avg (obs) Avg (comp) CC

Surface 1.5 6.1 13.1 13.3 0.96

Subsurface

epilimnion

hypolimnion

2.4

1.3

9.2

5.2

7.3

4.5

7.7

5.3

0.93

0.87
35

Table 12. Statistical comparison of 1994-95 observed and computed short term water level fluctuations (cm.)

Figure 17. Temporal evolution of simulated

versus observed temperature profiles, station

18M.

Gauge Npts Avg

(Obs.)

RMSV

(Obs.)

Avg

(Comp.)

RMSV

(Comp.)

 Diff. RMSD CC

 7023 17249 2.06 3.80 0.14 1.21 1.92 4.08 0.32

 7031 16571 0.69 2.83 0.65 1.95 0.04 2.55 0.48

 7044 16705 -0.29 6.95 0.81 4.23 -1.10 6.06 0.52

 7057 17243 0.34 3.69 0.07 2.01 0.27 3.45 0.39

 7068 15746 -0.72 3.53 -0.31 1.19 -0.41 3.42 0.28

 7072 17272 -0.99 3.12 -0.24 1.34 -0.74 3.04 0.34

 7079 16827 -0.38 10.11 -0.09 12.82 -0.29 13.58 0.32

 7096 16147 -2.27 5.74 0.06 3.09 -2.33 5.52 0.49
36

Figure 19. Simulated versus observed short-wave radiation for May-June 1994 at ISWS station at Beaver Island.

Figure 18. Simulated

versus observed

temperature, transects

offshore of Ludington,

MI (upper panel) and

Holland, MI (lower

panel).
37

Table 14. Statistical comparison of 1994-95 observed and computed wave periods (s).

Table 13. Statistical comparison of 1994-95 observed and computed wave heights (m).

-0.6

-0.6

-0.4

-0.4

-0.2

-0.2

-0.2

0.0

0.0

0.0

0.0

May-October 1994

-0.8

-0.6

-0.4

-0.4

-0.2

-0.2

-0.2

0.0

0.0

0.0

0.0

0.0

0.0

0.0

0.0

November 1994-April 1995 -0.8

-0.6

-0.6

-0.4

-0.4

-0.2

-0.2

0.0

0.0

0.0

0.0

0.0

May-October 1995

Normalized Stream Function

Blue - counterclockwise rotation

Yellow - clockwise rotation

Figure 20. Stream-function for summer and winter circulation for 1994-1995.

Buoy Npts Avg

(Obs.)

Avg

(Comp.)

Diff. RMS CC

45002 6993 0.77 0.85 -0.08 0.31 0.87

45007 7848 0.67 0.71 -0.04 0.23 0.91

45010 4529 0.57 0.55 0.02 0.23 0.88

Buoy Npts Avg

(Obs.)

Avg

(Comp.)

Diff. RMS CC

45002 6993 3.89 3.77 0.12 1.47 0.61

45007 7848 4.19 3.48 0.71 1.24 0.65

45010 4529 3.30 3.02 0.28 1.70 0.59
38

also developed programs to convert from one description to the other.

A ten segment Level II grid was developed (Fig. 21). Aggregated average surface heat flux (on an hourly time

scale) and average vertical temperature profiles (on a six-hourly time scale) were computed for each of the 10

segments for both the 1982-83 period and the 1994-95 period. In addition, inter-segment transports averaged over

1 day and 6 day intervals were computed.

We have also developed programs for aggregating the 3-d hydrodynamic model results for 3-d Level II grids. We

provided LLRS with daily averages and 6-day averages of horizontal and vertical inter-segment transport for the

10 segment Level II grid with 5 vertical layers: layer one 0-10m, layer two 10-20, layer three 20-30, layer four

30-50, and layer five 50-bottom.

All data sets and computer programs associated with Level II aggregation are on the Final Report CD.

6. INTERPRETATION/DISCUSSION

Our general conclusion is that overall, the models simulated the large scale thermal structure, circulation, and

waves quite realistically on the 5 km grid. We want to mention some problems, though. These problems reflect

principal limitations of hydrodynamical modeling approach: limitations of model physics, coarse spatial

discretization, errors in specification of forcing functions and initial conditions. First, the lack of an ice model

will be a serious problem if the model is applied during a year with normal or severe ice conditions. It will cause

both significant violations of lake’s heat balance and errors in calculating transfer of momentum from air to water

because of the difference in surface roughness of ice and water and momentum absorption by the ice.

Second, the exclusion of tributaries and outflow through the Straits of Mackinac (while having negligible effect

on lake circulation) can potentially influence the mass balance of constituents over a long period of time. The

flow in the vicinity of the Straits of Mackinac is probably the biggest challenge due to the lack of detailed

information on currents and temperature, which is needed on the open boundary. As Saylor and Sloss (1976),

Schelske et al. (1976), and Beeton and Saylor (1995) showed, there is a significant intraseasonal and interseasonal

variability of currents and temperature in that area. In summer, there is a persistent flow of deep water from Lake

Huron into Lake Michigan caused by temperature gradients in the vicinity of the Straits. It is not clear how far

Lake Huron water penetrates into Lake Michigan, although it is probably fair to say that it should effect the Lake

Michigan’s northern basin water quality since the amount of water in this inflow is sufficient to decrease the

residence time of Lake Michigan from 137 to 69 years (Quinn, 1977). Therefore, accurate time-dependent

simulation of water exchange through the Straits of Mackinac will most probably require a simulation of the

entire Lake Michigan - Lake Huron system since neither 1982-83 nor 1994-95 current and temperature observations are available. Other alternatives for including the Straits of Mackinaw in the model are 1) to specify the

measured water level at that grid cell and allow the model to compute inflow and outflow or 2) to specify the

time-dependent vertical profile of current and temperature at the Straits (this would require that appropriate

measurements be available).

Third, the model did not perform as well in the thermocline area as it performed near the surface. The simulated

thermocline was too diffuse. To study the effect of vertical resolution on the vertical temperature gradients, we

carried out a model run with 39 sigma levels, e.g. double the vertical resolution. In this run we noticed only

moderate improvement in the thermocline area. We also ran the model with a zero horizontal diffusion to test for

artificial diffusion along sigma surfaces. Again, we did not notice a significant improvement in model results. On

the other hand, experiments with the 1-dimensional version of the model showed that the Mellor-Yamada scheme

can provide a sharp thermocline (which is also sensitive to the choice of extinction coefficient which possess

significant spatial and temporal variability in large lakes but was kept constant in our calculations). The extinction coefficient determines how much incoming heat is absorbed in each vertical layer in each grid square in the

hydrodynamic model. Extinction coefficient vary with concentration of chlorophyll, suspended solids/turbidity,
39

Figure 21. Level II model grid.
40

and dissolved organic carbon. At the time of these model runs, there was not sufficient information about the

temporal and spatial variability of these parameters to try to include a spatially and temporally variable extinction

coefficient. We think that the problem of the thermocline being too diffuse will most likely be solved by increasing the horizontal resolution. We feel that a horizontal resolution of 2 km should improve the results for Lake

Michigan, but even higher resolution may be required for accurate simulation of vertical thermal structure.

Unfortunately, 2 km and finer horizontal resolution is still beyond the limits of available computer resources,

especially for water quality models. Another option would be to use a z-level type of model instead of

sigma-level model. There are only a few direct comparisons between different types of models. For a problem of

wind-induced upwelling and wave propagation on the existing thermocline, the results from z-level and σ-level

models were very similar (Beletsky et al., 1997). Although we cannot prove it, we believe that a sigma-level

model is suited better for a simulation of the developing (decaying) thermocline on the sloping bottom because of

the superior vertical resolution in the coastal areas of the Great Lakes.

The horizontal resolution of the model was only sufficient for a description of large-scale circulation patterns.

Point to point comparison of observed and simulated currents (as indicated by smallest Fourier norms, Fn~0.7-0.8)

was successful only in the fall-winter months in the southern basin, which is characterized by a smooth bathymetry. Comparison of observed and predicted progressive vector diagrams (Figure 10) showed that the modeled

currents agreed qualitatively with the observed long-term flow in the southern basin, with modeled current speeds

approximately 10-30% higher than observed speeds. Obviously, 5 km horizontal grid resolution is too coarse to

adequately resolve baroclinic motions in summer. The point to point comparison was the worst (Fn~1.0) in the

areas of strong depth gradients where model resolution was inadequate. The circulation was stronger in winter

than in summer, and also more organized and cyclonic, which is in agreement with observations. Since the lake is

essentially homogeneous in winter, there are two most probable explanations: existence of stronger cyclonic wind

vorticity in winter, or existence of residual mean cyclonic circulation driven by nonlinear interactions of topographic waves. More research is needed to clarify the role of each factor.

Finally, there is potential for improvement in the specification of forcing functions and initial conditions. Currently, meteorological observations are 30-50 km apart, which means that we may miss details of mesoscale

meteorology. The 1994-95 observation network had more stations than 1982-83, but model predictions did not

improve noticeably. Perhaps the reason is that almost all new added stations were land stations, while the most

significant changes in the air characteristics should occur over water. Initialization of 3-dimensional temperature

field also presents a significant challenge because only surface temperature from two NDBC buoys is usually

available for that purpose supplemented by episodic satellite observations.

7. CONCLUSIONS

We have described the hydrodynamic modeling components of the LMMBS including the surface flux model, the

wave model, and the three dimensional circulation model. We expect that these models will provide estimates of

the physical environment in the lake with a spatial resolution of tens of kilometers and a temporal resolution of

hours. These space and time scales are significant for many of the important chemical and biological processes

that are being modeled as part of the LMMBS, particularly sediment resuspension and transport and lower food

web modeling. Hydrodynamic models generally have very few system-specific adjustable parameters, so they are

much easier to apply to new conditions, outside the conditions encountered during the project years. It is our

conclusion that hydrodynamic and water quality models are essential for assessing regional impacts of the types

of pollutants being studied in LMMBS. As a result of this study, we have identified several areas for future

research including coupling of near-field and far-field models, closer linkages between hydrodynamic and water

quality models, and more coordinated sampling programs for physical, biological, and chemical parameters. It is

our hope that the LMMBS modeling framework, which is based on the hydrodynamic models, will provide a

sound basis for the development of management tools for Lake Michigan.

8. REFERENCES
41

Allender, J.H. Comparison of model and observed currents in Lake Michigan. Journal of Physical Oceanography

7:711-718 (1977).

Allender, J.H., and J.H. Saylor. Model and observed circulation throughout the annual temperature cycle of Lake

Michigan. Journal of Physical Oceanography 9:573-579 (1979).

Assel, R.A., F.H. Quinn, G.A. Leshkevich, and S.J. Bolsenga. Great Lakes Ice Atlas. National Oceanic and

Atmospheric Administration, Great Lakes Environmental Research Laboratory, Ann Arbor, MI (1983).

Assel, R. A., C.R. Snider, and R. Lawrence. Comparison of 1983 Great Lakes winter weather and ice conditions

with previous years. Monthly Weather Review 113:291-303 (1985).

Assel, R. A., J. E. Janowiak, S. Young, and D. Boyce. Winter 1994 weather and ice conditions for the Laurentian

Great Lakes. Bulletin of the American Meteorological Society 77:71-88 (1996).

Bedford, K.W., and D.J. Schwab, 1990. Preparation of real-time Great Lakes Forecasts. Cray Channels Summer

1990:14-17 (1990).

Beeton, A.M., and J.H. Saylor. Limnology of Lake Huron. In: The Lake Huron Ecosystem: Ecology, Fisheries and

Management, M. Munawar, T. Edsall and J. Leach (eds.), SPB Academic Publishing, Amsterdam, p. 1-37 (1995).

Beletsky, D., W.P. O’Connor, D.J. Schwab, and D.E. Dietrich. Numerical simulation of internal Kelvin waves and

coastal upwelling fronts. Journal of Physical Oceanography 27(7):1197-1215 (1997).

Bennett, J.R. On the dynamics of wind-driven lake currents. Journal of Physical Oceanography 4(3):400-414

(1974).

Blumberg, A.F. Turbulent mixing processes in lakes, reservoirs and im poundments. In Physics-based modeling

of lakes, reservoirs, and impoundments, W.G. Gray, (Ed.). ASCE, New York, NY, pp. 79-104 (1986).

Blumberg, A.F., and G.L. Mellor. A description of a three-dimensional coastal ocean circulation model. In:

Three dimensional Coastal Ocean Models, Coastal and Estuarine Sciences, 5, N.S. Heaps (ed.). American

Geophysical Union, Washington, DC, pp. 1-16 (1987).

Boyce, F.M., M.A. Donelan, P.F. Hamblin, C.R. Murthy, and T.J. Simons. Thermal structure and circulation in the

Great Lakes. Atmosphere-Ocean 27(4):607-642 (1989).

Campbell, J. E., A.H. Clites, and G.M. Green. Measurements of ice motion in Lake Erie using satellite-tracked

drifter buoys. NOAA Data Report, ERL GLERL-30, 22pp. (1987).

Charnock, H. Wind stress on a water surface. Quarterly Journal of the Royal Meteorological Society 81:639

(1955).

Csanady, G.T. Circulation in the Coastal Ocean. D. Reidel Co., Dordrecht, Holland, 279 pp. (1982).

Donelan, M.A. On the fraction of wind momentum retained by waves. In: Marine Forecasting, Predictability

and Modeling in Ocean Hydrodynamics, J.C.Nihoul (ed.) Elsevier, Amsterdam (1979).

Gottlieb, E.S., J.H. Saylor, and G.S. Miller. Currents and temperatures observed in Lake Michigan from June

1982 to July 1983. NOAA Technical Memorandum ERL GLERL-71, 45 pp. (1989).

Hasselmann, K., T.P. Barnett, E. Bouws, H. Carlson, D.E. Cartwright, K. Enke, J.A. Ewing, H. Gienapp, D.E.
42

Hasselmann, P. Kruseman, A. Merrburg, P. Muller, D.J. Olbers, K. Richter, W. Sell, and H. Walden. Measurements of wind-wave growth and swell decay during the Joint North Sea Wave Project (JONSWAP). Dtsch.

Hydrogr. Z., A12, 95 pp. (1973).

Jerlov, N.G. Marine Optics. Elsevier, Amsterdam, 231 pp. (1976).

Johnson, F.R., D.E. Boyce, J.A. Bunn, and J.L. Partain. In search of the perfect wave - a new method to forecast

waves on the Great Lakes. NWS Eastern Region Technical Attachment No. 92-9A, NOAA National Weather

Service, Silver Spring, MD. 12 pp. (1992).

Kuan, C., K.W. Bedford, and D.J. Schwab. A preliminary credibility analysis of the Lake Erie portion of the

Great Lakes Forecasting System for springtime heating conditions. In: Quantitative Skill Assessment for Coastal

Ocean Models number 47 in Coastal and Estuarine Studies, D.R. Lynch and A.M. Davies (eds.). American

Geophysical Union, pp. 397-423 (1994).

Liu, P.C., and D.J. Schwab. A comparison of methods for estimating u* from given uz and air-sea temperature

differences. Journal of Geophysical Research 92(C6):6488-6494 (1987).

Liu, P.C., D.J. Schwab, and J.R. Bennett. Comparison of a two-dimensional wave prediction model with synoptic

measurements in Lake Michigan. Journal of Physical Oceanography 14(9):1514-1518 (1984).

Long, P.E., and W.A. Shaffer. Some physical and numerical aspects of boundary layer modeling. NOAA Technical

Memorandum NWS TDL-56 (COM75-10980) (1975).

Lyons, W.A. Low level divergence and subsidence over the Great Lakes in summer. Proceedings of the 14th

Conference on Great Lakes Research, International Association for Great Lakes Research, pp. 467-487 (1971).

McCormick, M.J., and G.A. Meadows. An intercomparison of four mixed layer models in a shallow inland sea.

Journal of Geophysical Research 93(C6):6774-6788 (1988).

Mellor, G.L. An equation of state for numerical models of oceans and estauries. Journal of Atmospheric and

Oceanic Technology 8:609-611 (1991).

Mellor, G.L., and T. Yamada. Development of a turbulence closure model for geophysical fluid problems. Rev. of

Geophysical Space Physics 20(4):851-875 (1982).

O’Connor, W.P., and D.J. Schwab. Sensitivity of Great Lakes Forecasting System Nowcasts to Meteorological

Fields and Model Parameters. Proceedings, 3rd International Conference on Estuarine and Coastal Modeling.

ASCE Waterway, Port, Coastal and Ocean Division, 149-157 (1994).

Pettersen, S., and P.A. Calabrese. On some weather influences due to warming of the air by the Great Lakes in

winter. Journal of Meteorology 16:646-652 (1959).

Phillips, D.W., and J.G. Irbe. Lake to land comparison of wind, temperature, and humidity on Lake Ontario

during the International Field Year for the Great Lakes (IFYGL). Atmospheric Environment Service Report CLI2-77, Downsview, Ontario (1978).

Quinn, F.H. Annual and seasonal flow variations through the Straits of Mackinac. Water Resources Research

13:137-144 (1977).

Rao, D.B., and T.S. Murty. Calculation of the steady-state wind-driven circulation in Lake Ontario. Arch. Meteor.

Geophys. Bioklim., A19:195-210 (1970).
43

Resio, D.T., and C.L. Vincent. Estimation of winds over the Great Lakes. Journal Waterway Port Coastal Ocean

Division, ASCE, 102:265-283 (1977).

Saylor, J.H., and P.W. Sloss. Water volume transport and oscillatory current flow through the Straits of Mackinac. Journal of Physical Oceanography 6(2):229-237 (1976).

Saylor, J.H., J.C.K. Huang, and R.O. Reid. Vortex modes in Lake Michigan. Journal of Physical Oceanography

10(11):1814-1823 (1980).

Schelske, C.L., E.U. Stoermer, J.E. Gannon, and M.S. Simmons. Biological, chemical, and physical relationships

in the Straits of Mackinac. EPA-600/3-76-095, Environmental Research Laboratory-Duluth, Duluth, MN (1976).

Schwab, D.J. Simulation and forecasting of Lake Erie storm surges. Monthly Weather Review 106:1476-1487

(1978).

Schwab, D.J. Numerical simulation of low-frequency current fluctuations in Lake Michigan. Journal of Physical

Oceanography 13(12):2213-2224 (1983).

Schwab, D.J. A review of hydrodynamic modeling in the Great Lakes from 1950-1990 and prospects for the

1990’s. In: A.P.C. Gobas and A. McCorquodale (eds.), Chemical Dynamics in Fresh Water Ecosystems. Lewis

Publishers, Ann Arbor, MI, 41-62 (1992).

Schwab, D.J., and K.W. Bedford. Initial implementation of the Great Lakes Forecasting System: a real-time

system for predicting lake circulation and thermal structure. Water Pollution Research Journal of Canada 29(2/

3):203-220 (1994).

Schwab, D.J., J.R. Bennett, P.C. Liu, and M.A. Donelan. Application of a simple numerical wave model to Lake

Erie. Journal of Geophysical Research 89(C3):3586-3592 (1984a).

Schwab, D.J., J.R. Bennett, and E.W. Lynn. A two-dimensional lake wave prediction system. NOAA Technical

Memorandum ERL GLERL-51, Great Lakes Environmental Research Laboratory, Ann Arbor, MI, 70 pp. (1984b).

Schwab, D.J., and J. A. Morton. Estimation of overlake wind speed from overland wind speed: a comparison of

three methods. Journal of Great Lakes Research 10(1):68-72 (1984).

Schwab, D.J., and D.L. Sellers. Computerized bathymetry and shorelines of the Great Lakes. NOAA Data Report,

ERL GLERL-16, 13 pp. (1980).

Simons, T.J. Verification of numerical models of Lake Ontario: I. Circulation in spring and early summer. Journal

of Physical Oceanography 4: 507-523 (1974).

Simons, T.J. Circulation models of lakes and inland seas. Canadian Bulletin of Fisheries and Aquatic Science

203, 146 pp. (1980).

Simons, T.J. Reliability of circulation models. Journal of Physical Oceanography 15:1191-1204 (1985).

Thiessen, A.H. Precipitation average for large areas. Monthly Weather Review 39:1082- 1084 (1911).

Wyrtki, K. The average annual heat balance of the North Pacific Ocean and its relation to ocean circulation.

Journal of Geophysical Research 70:7547-4599 (1965).
44

9. PRODUCTS

9.1 Publications

Beletsky, D., J.H. Saylor, and D.J. Schwab. Mean circulation in the Great Lakes. Journal of Great Lakes Research (in review) (1998).

Beletsky, D., K.K. Lee, and D.J. Schwab. Large-scale circulation. In: D. Lam (ed.) Climatic Effects on Lake

Hydrodynamics and Water Quality (In review) (1998).

Schwab, D.J. and D. Beletsky. Propagation of Kelvin waves along irregular coastlines in finite-difference

models. Advances in Water Resources ( in press) (1998).

Beletsky, D., and D.J. Schwab. Modeling thermal structure and circulation in Lake Michigan. Estuarine and

Coastal Modeling, Proceedings, 5th International Conference, October 22-24, 1997, Alexandria, VA, p. 511-522

(1998).

Schwab, D.J., and D. Beletsky. Hydrodynamic modeling program in the Lake Michigan Mass Balance Project.

Workshop on research needs for coastal pollution in urban areas, October 16-17, Milwaukee, WI, p. 14-22 (1997).

Beletsky, D., W.P. O’Connor, D.J. Schwab, and D.E. Dietrich. Numerical simulation of internal Kelvin waves and

coastal upwelling fronts. Journal of Physical Oceanography 27:1197-1215 (1997) .

Beletsky, D., K.K. Lee, and D.J. Schwab. Recent advances in hydrodynamic modeling of the Great Lakes. Proceedings of the XXVII IAHR Congress, August15-18, San Francisco, CA, p. 925-930 (1997).

Beletsky, D., W.P. O’Connor, and D.J. Schwab. Hydrodynamic modeling for the Lake Michigan Mass Balance

Project. In G. Delic and M.F. Wheeler (eds.), Next Generation Environmental Models Computational Methods,

Proceedings of a U.S.EPA sponsored workshop at the National Environmental Supercomputing Center, August

7-9, 1995, Bay City, MI. SIAM, Philadelphia, PA, p.125-128 (1997).

Schwab, D.J., D. Beletsky, W.P. O’Connor, and D.E. Dietrich. Numerical simulation of internal Kelvin waves

with z-level and sigma level models. In: M.L. Spaulding and R.T. Cheng (eds.), Estaurine and Coastal Modeling,

Proceedings, 4th International Conference, October 26-28, 1995, San Diego, CA. American Society of Civil

Engineers, New York, NY, p. 298-312 (1996).

9.2 Presentations

Beletsky, D. Climate and large lakes dynamics (Overview of CILER Task II). CILER Formal Review, July 15-16,

1998. Ann Arbor, MI.

Schwab, D. J., and D. Beletsky. Hydrodynamic modeling in the Lake Michigan Mass Balance Project. LMMBP

Modeling Program External Review, June 23-25, 1998. Southgate, MI.

Beletsky, D., and D.J. Schwab. Simulation of the interannual variability of circulation and thermal structure in

Lake Michigan. The 41st Conference of the International Association for Great Lakes Research, May 18-22,

1998, Hamilton, Ontario, Canada.

Schwab, D. J., and D. Beletsky. Modeling transport and mixing in the Lake Michigan Mass Balance Project. The

1998 Ocean Science Meeting, February 9-13, 1998, San Diego, CA.
45

Schwab, D. J. and D. Beletsky. Hydrodynamic modeling program in the Lake Michigan Mass Balance Project.

Workshop on research needs for coastal pollution in urban areas, October 16-17, 1997, Milwaukee, WI.

Beletsky, D., and D.J. Schwab. Modeling thermal structure and circulation in Lake Michigan. The 5th International Conference for Estuarine and Coastal Modeling, October 22-24, 1997, Alexandria, VA.

Schwab, D.J., and D. Beletsky. Internal Kelvin Wave Propagation in Finite Difference Models. The 5th International Conference for Estuarine and Coastal Modeling, October 22-24, 1997, Alexandria, VA.

Beletsky, D., and D.J. Schwab. Toward coupled ice-circulation model of Lake Michigan. Remote Sensing and

Modeling Great Lakes Ice Workshop, October 8-9, 1997, Alexandria, VA.

Beletsky, D., K.K. Lee, and D.J. Schwab. Recent advances in hydrodynamic modeling of the Great Lakes. Proceedings of the XXVII IAHR Congress, August 15-18, 1997. San Francisco, CA.

Schwab, D.J., and D. Beletsky. Hydrodynamic modeling of Lake Michigan. The 40th Conference of the International Association for Great Lakes Research, June 1-5, 1997. Buffalo, NY.

Beletsky, D., J.H. Saylor, and D.J. Schwab. On the mean circulation in the Great Lakes. The 40th Conference of

the International Association for Great Lakes Research, June 1-5, 1997. Buffalo, NY.

Beletsky, D., and D.J. Schwab. Visualization techniques for lake hydrodynamics studies. CGLAS/CILER Mini

Symposium, February 4, 1997, Ann Arbor, MI.

Beletsky, D., and D.J. Schwab. Numerical modeling of circulation in Lake Michigan. CoastWatch users meeting,

September 18, 1996, Ludington, MI - Manitowoc, WI.

Schwab, D.J., and D. Beletsky. Application of POM in the Great Lakes. Princeton Ocean Model Users Meeting,

June 9-12, 1996. Princeton, NJ.

Beletsky, D., and D.J. Schwab. Modeling of the annual cycle of thermal structure and circulation in Lake Michigan. The 1996 Ocean Science Meeting, February 12-16, 1996. San Diego, CA.

Beletsky, D., and D.J. Schwab. Developments in the Great Lakes circulation modeling. CGLAS/CILER Mini

Symposium, February 2, 1996. Ann Arbor, MI.

Schwab, D.J, D. Beletsky, W.P. O’Connor, and D.E.Dietrich. Numerical simulation of internal Kelvin waves with

z-level, and sigma level models. The 4th International Conference on Estuarine and Coastal Modeling, October

26-28, 1995. San Diego, CA.

Beletsky, D., W.P. O’Connor, and D.J. Schwab. Hydrodynamic modeling for the Lake Michigan mass balance

project. U.S. EPA Workshop on Next Generation Environmental Models Computational Methods, August 7-9,

1995, Bay City, MI.

Beletsky, D., W.P. O’Connor, and D.J. Schwab, 1995. Numerical simulation of internal Kelvin waves in lakes.

The 38th Conference of the International Association for Great Lakes Research, May 28-June 1, 1995, East

Lansing, MI.

O’Connor, W.P., D. Beletsky, and D.J. Schwab. Internal Kelvin waves in lakes. CGLAS/CILER Mini Symposium,

January 2, 1995, Ann Arbor, MI.
46

9.3 Final Report CD: Graphics, Programs,Observations, and Selected Model Results

 (including averages and aggregation)

The Final Report CD supplements this report with extensive collection of plots of model results and observations

in the form of HTML and GIF files. It also contains all IDL, FORTRAN, HTML, and C-shell programs, observed

data, selected model results (mainly various types of averages), and this report in the pdf format. The Final Report

CD contains the following directories:

data observation files

graphics GIF files

html HTML files

model_results 1D files, averages, aggregates

programs aggregate, graphics, interpolation, pomgl, wmgl and XDR programs

report this report

The full CD directory structure is as follows:

 |---data various observation data

 | |

 | |-ctd temperature observations (GLNPO CTD surveys,1994-95)

 | |-intakes temperature and turbidity at water intakes (1994-95)

 | |-marobs marobs data files (1982-83; 1994-95)

 | |-moorings temperature and currents at moorings (1982-83; 1994-95)

 | |-ndbc NDBC data files (1982-83; 1994-95)

 | |-transects temperature observations (USGS CTD surveys, 1994-95)

 | \-wlevels water level data (1982-83; 1994-95)

 |

 |---graphics plots generated for this study

 | |

 | |---averages 30 day average plots

 | | |

 | | |-meteo plots of meteorological data

 | | \-model plots of model data

 | |

 | |-ctd GLNPO survey temperature plots

 | |-profiles USGS survey temperature transect plots

 | |

 | \--time_series 60 day time series plots

 | |

 | |-currents mooring current and temperature plots

 | |-meteo meteorological plots

 | |-temperature water intake and NDBC temperature plots

 | |-waves wave height model results vs observed plots

 | \-wlevels water level model results vs observed plots

 |

 |---html html files

 | |

 | \---images images used in the html pages

 | |

 | \-stations locations of meteorological stations

 |

 |-index.html the starting point for your web browser

 |

 |--model_results selected model results

 | |
47

 | |-1d XDR files with time-series of simulated variables

 | |-aggregate XDR files with aggregated transports

 | \-averages XDR files used in 30 day averages plots

 |

 |---programs programs used in the lmmbs study

 | |

 | |---idl idl programs

 | | |

 | | |-averages programs used to plot 30 day averages

 | | |-buoys programs to plot temperature at NDBC buoys

 | | |-ctd programs to plot temperature obtained in GLNPO CTD surveys

 | | |-currents program used to generate 60 day temp & current meter plots

 | | |-intakes programs to plot temperature at municipal water intakes

 | | |-maps programs to plot observation networks

 | | |-meteo program used to generate meteorological plots

 | | |-miscellan observations coordinates and miscellaneous subroutines

 | | |-solar programs to plot short-wave solar radiation

 | | |-transects program used to generate temp. profiles

 | | |-waves programs used to generate 60 day wave plots

 | | \-wlevels programs used to generate 60 day water level plots

 | |

 | |---interp programs for interpolation of meteorological observations

 | |---pomgl hydrodynamic model code and shell scripts

 | |---wmgl wave model code and shell scripts

 | \---xdr XDR I/O programs

 |

 \---report final report in Adobe Acrobat (pdf) format

9.3.1 Model Codes

Lake Michigan bathymetric grids are in files mich2.dat, mich5.dat, and mich10.dat in the programs/pomgl

subdirectory. Hydrodynamic model code includes files pomgl-main.f, pomgl-subs.f, pom-subs.f, and pomgl.h

(common block file). Hydrodynamic model initialization code is in the file pomglrin.f (common block file

pomglrin.h). Interpolation model code is in the file interpgl.f. Idealized boundary conditions (wind stress and heat

flux) can be created with the program pomglflx.f. Wave model code is in the file wmgl.f in the programs/wmgl

subdirectory. A description of a sample single period model run is given below. Section 9.3.2 describes how a

multi-year simulation (consisting of several periods) was implemented.

As an example of how to run the hydrodynamic model, we will describe the steps required for one 60-day period.

The example covers the 60 day period from May 1, 1994 0000Z to June 30, 1994 0000Z. Meteorological forcing

functions are based on station observations from airports, NDBC buoys, and LMMBS atmospheric deposition

monitoring stations. The initial condition is read in from the file m94061.rst.

There are two steps involved in the hydrodynamic model simulation. First, the station meteorological data is

interpolated to the hydrodynamic grid. Second, the hydrodynamic model is run.

The wave model can also be run for the same time period using the interpolated meteorological data created in

step 1. The wave model run is independent of the hydrodynamic model run.
48

Approximate execution times on HP C160 workstation:

interpce: 10 mins.

pomgl: 6 hrs.

wmgl: 10 mins.

------------------------------------------------------------------------------

Step 0. Compile C interface to XDR I/O routines on your system

-------------------------------------------------------------------------------

0.1 Copy source code to working directory: hpxdrc.c (or sunxdrc.c) from programs/xdr subdirectory

0.2 Rename to xdrc.c

0.3 Compile interface routines: cc -c xdrc.c

-----------------------------------------------------------------------------------

Step 1. Interpolate station meteorological data to hydrodynamic grid

-----------------------------------------------------------------------------------

1.1 Copy source code and input files to working directory from programs/interp subdirectory:

interpce.f

interpsubs.f

xdrsubs.f

interp.in

mich5.dat

stations.dat

m94121.mmo

1.2 Compile interpolation program ‘interpce’.

f77 -o interpce -O interpce.f interpsubs.f xdrsubs.f xdrc.o

1.3 Run program interpce.

interpce < interp.in > interp.out

1.4 The following files are created:

interp.adj *

interp.e *

interp.s *

m94121.at

m94121.cl

m94121.dp

m94121.wu

m94121.wv

m94121.wt *

* = temporary file which is not required for pomgl model run
49

----------------------------------------

Step 2. Run hydrodynamic model

----------------------------------------

2.1 Copy source code and input files for hydrodynamic model to working directory from programs/pomgl

subdirectory:

pomgl-main.f

pomgl-subs.f

pom-subs.f

pomgl.h

xdrc.o

pomgl.in

mich5.dat

2.2 Compile hydrodynamic model

f77 -o pomgl -O pomgl-main.f pomgl-subs.f pom-subs.f xdrc.o

2.3 Run hydrodynamic model

pomgl < pomgl.in > pomgl.out

2.4 The following xdr files are created:

m94121.1d

m94121.utm

m94121.vtm

m94121.rst

m94121.tmz

-------------------------------------------------------------

Step 3. Run wave model (Does not require step 2)

-------------------------------------------------------------

3.1 Copy source code and input files for wave model to working directory from programs/wmgl subdirectory:

wmgl.f

wmgl.in

xdrc.o

mich5.dat

3.2 Compile wave model

f77 -o wmgl -O wmgl xdrc.o

3.3 Run wave model

wmgl < wmgl.in > wmgl.out

3.4 The following xdr files are created:

m94121.wvh

m94121.wvp

m94121.wvd

wmgl.rst

9.3.2 Shell Scripts

UNIX C-shell scripts were used during sorting/interpolation of meteorological data, and model runs. Interpolation script is in the file interp.csh in the programs/interp subdirectory.
50

Model run script (in the file s_m82h_crunch1 or s_m94z11_crunch1) in the programs/pomgl subdirectory makes

changes in the dimensions of arrays, archives current model code, and then compiles the program. A typical

model run consists of several separate 60-day periods, therefore, several control files (pomgl.in) are needed. The

first control file is created in the beginning of the script, and is used as a template for creation of subsequent

control files by means of stream editor sed. After that, the script starts an initialization program to specify initial

temperature and velocity fields in the lake, and begins execution after that. In the 1994-95 case, we ran the model

in a diagnostic mode for the first three months. Therefore, the second 60-day period consists of a 30-day barotropic run followed by a separate 30-day baroclinic run. The script combines two 30-day periods into a regular size

period and the proceeds to normal 60-day period runs. In the end of the script we combine all files with extension

‘1d’ and ‘tmz’.

9.3.3 Programs for input/output analysis (IDL)

 Programs that were used to create graphic output on Final Report CD in programs/idl subdirectory:

a82atp.pro - plots 30-day average air temperature (1982-83)

a82clp.pro - plots 30-day average cloud cover (1982-83)

a82cp.pro - plots 30-day average currents (1982-83)

a82dpp.pro - plots 30-day average dew point (1982-83)

a82hfp.pro - plots 30-day average total surface heat flux (1982-83)

a82mwp.pro - plots 30-day maximum wave height (1982-83)

a82sfp.pro - plots 30-day average transport stream function (1982-83)

a82tp.pro - plots 30-day average surface water temperature (1982-83)

a82uvp.pro - plots 30-day average wind field (1982-83)

a82whp.pro - plots 30-day average wave height (1982-83)

a82wsp.pro - plots 30-day average wind stress (1982-83)

a94atp.pro - plots 30-day average air temperature (1995-95)

a94clp.pro - plots 30-day average cloud cover (1995-95)

a94cp.pro - plots 30-day average currents (1995-95)

a94dpp.pro - plots 30-day average dew point (1995-95)

a94hfp.pro - plots 30-day average total surface heat flux (1995-95)

a94mwp.pro - plots 30-day maximum wave height (1995-95)

a94sfp.pro - plots 30-day average transport stream function (1995-95)

a94tp.pro - plots 30-day average surface water temperature (1995-95)

a94uvp.pro - plots 30-day average wind field (1995-95)

a94whp.pro - plots 30-day average wave height (1995-95)

a94wsp.pro - plots 30-day average wind stress (1995-95)

plot82buoy.pro - plots 60-day time-series of surface water temperature at NDBC buoys 45002 and 45007

(1982-83)

plot94buoy.pro - plots 60-day time-series of surface water temperature at NDBC buoys 45002 and 45010

(1994-95)

plot94buoy2.pro - plots 60-day time-series of surface water temperature at NDBC buoy 45007(1994-95)

ctd.pro - plots temperature profiles at CTD survey locations (1994-95)

tcurrent.pro plot_intake.pro- plots 60-day time-series of observed versus simulated temperatures at water intakes

locations (1994-95)

plot_marobs.pro - plots 60-day time-series of observed air temperature, water temperature, dew point, cloud

cover, wind speed and direction

plot_solar.pro - plots 60-day time-series of observed versus simulated short-wave solar radiation (1994-95)

plot_transect.pro - plots observed versus simulated water temperature at several transects (1994-95)

plot_82wave.pro - plots 60-day time-series of observed versus simulated wind wave height, period and direction

(1982-83)

 plot_94wave.pro - plots 60-day time-series of observed versus simulated wind wave height, period and direction

(1994-95)
51

plot_wl.pro - plots 60-day time-series of observed versus simulated water level elevations

table_01.txt - location of moorings, water level gauges and NDBC buoys (1982-83)

table_03.txt - location of moorings, water level gauges and NDBC buoys (1994-95)

There are several files in directory programs/idl/miscellaneous with coordinates of other observations:

ctdlog.dat - CTD surveys

intakes.dat - water intakes

stations1.dat - meteorological observations

9.4 Model input and output on CD - binary (XDR) and ASCII

There are 10 CDs with 1982-83 and 12 CDs with 1994-95 binary (XDR) model input and output, and also ASCII

input and output model files. Each CD corresponds to a single 60 day period and contains three directories:

circulation, meteorology, and waves. Fields are stored at hourly (meteorology and waves), 3-hourly (currents),

and 6-hourly (water temperature) intervals on the 5 km Lake Michigan hydrodynamic grid. For the 1982-83

simulation, the period covered is from 31 March 1982 through 20 November 1983 (600 days). There are 10 files

for each field, corresponding to the 10 consecutive 60 day periods during this interval. For the 1994-95 simulation, the period covered is from 1 Jan. 1994 to 21 Dec. 1995 (720 days). There are 12 files for each field,

corresponding to the 12 consecutive 60 day periods during this interval. The file names are mYYJJJ.ppp where

YY is the year for the first data field in the file (94 or 95), JJJ is the day of the year (1-365), and ppp indicates

which parameter is stored in the file as shown above. The filenames for the 1982-83 simulations are: m82090,

m82150, m82210, m82270m, m82330, m83025, m83085, m83145, m83205, and m83265. The filenames for the

1994-95 simulations are: m94001, m94061,m94121, m94181, m94181, m94241, m94301, m94361, m95056,

m95116, m95176, m95236, and m95296. We illustrate the contents of the CD set using first 60 day period of

1994-95 simulations as a template.

The directory 'meteorology' contains four ASCII files:

m94001.mmo - meteorological observations file

m94001.in - input file for interpolation program interpce.f

m94001i.out - output of the interpolation program

m94001.wt - lake-averaged climatological surface water temperature.

The rest of the directory contains XDR files with hourly meteorological data at all grid points on the 5 km Lake

Michigan grid:

m94001.at - air temperature at 10m height (C)

m94001.cl - cloud cover (percent)

m94001.dp - dew point at 10m height (C)

m94001.wu - eastward component of wind velocity at 10m height (m/sec)

m94001.wv - northward component of wind velocity at 10m height (m/sec)

The directory 'waves' containes 2 ASCII files:

m94001.win - wave model input file

m94001w.out - wave model output file.

There are also three XDR files with hourly wave data:

m94001.wvh - significant wave height from wave model (m)

m94001.wvp - peak energy wave period from wave model (sec)

m94001.wvd- wave direction from wave model (degrees, 0 = toward east)
52

The directory 'circulation' contains four ASCII files:

m94001.in - hydrodynamic model input file

m94001.out - hydrodynamic model output file

m94001.err - hydrodynamic model error message file

m94001.outrin -output of the initialization program

The last file occurs only in the beginning of the model run as well as the XDR initialization file m94001.ini which

is an equivalent of the XDR restart file m94001.rst generated in the end of each model run.

There are 11 hourly XDR files (we omit the model run suffix here for brevity):

m94001.t.000.000.101 - surface water temperature (C)

m94001.tx - X-component of momentum flux (m**2/s**2)

m94001.ty - Y-component of momentum flux (m**2/s**2)

m94001.hf - Net surface heat flux out of lake (W/m**2)

m94001.swr - Short wave radiation (W/m**2)

m94001.lwr - Long wave radiation (W/m**2)

m94001.shf - Sensible heat flux (W/m**2)

m94001.lhf - Latent heat flux (W/m**2)

m94001.eg - Surface water level (m)

m94001.1d - Miscellaneous model output at specified locations

m94001.tmz - Lake-averaged vertical temperature profile (C)

Currents are saved 3-hourly, and water temperature is saved 6-hourly in the following XDR files:

m94001.t - water temperature (C)

m94001.um - X component of 3D current field

m94001.vm - Y component of 3D current field

m94001.utm - X component of 2D current field (depth-averaged)

m94001.vtm - Y component of 2D current field (depth-averaged)

For the 1994-95 simulation, the water temperature in the hydrodynamic model of Lake Michigan is held constant

at 2 degrees for the first 90 days of the run, so surface water temperature fields will indicate a constant 2 degrees

at all grid points during this period. During the first 90 days the model does not calculate surface heat flux, so

these values are 0.

The wu, wv, at, dp and cl files all start at hour 00Z on the indicated day and include 1441 hourly fields. All the

other fields are based on model output and start at 01Z on the indicated day and include 1440 hourly fields.

All fields except tx,um,utm and ty,vm,vtm represent values at the centers of the non-land grid cells in the Lake

Michigan 5 km hydrodynamic grid. There are 2318 non-land cells. The tx, um and utm fields are calculated at upoints in the computational lattice (the center of the west side of each grid cell). The ty, vm, and vtm

fields are calculated at v-points in the computational lattice (the south side of each grid cell). There are 2162

interior u-points and 2223 interior v-points respectively. Note also that the sign convention for tx and ty is the

opposite of wu and wv, that is a positive tx is a westward stress and a positive ty is a southward stress.

The data are stored in XDR format, which is a machine independent method for representing data. Each type of

UNIX machine has its own subroutines for converting XDR data into native format. You can usually find out

more about XDR on your UNIX machine by typing 'man xdr'. The UNIX system routines to read and write XDR

files are usually callable only from a C program, so if you want to use a Fortran program to read the files it is

necessary to provide a set of C routines that are callable from Fortran which will link to the system's native XDR

routines. We have provided a set of these routines for several different UNIX systems in the files hpxdrc.c (HP),

sunxdrc.c (SUN), and decxdrc.c (DEC) in the programs/xdr area on the graphics CD. These routines need to be

compiled on your system with the C compiler and the object file saved for inclusion with your Fortran program.
53

The command for compiling the XDR linkages is ususally 'cc -c hpxdrc.c', or the name of the appropriate routines

for your system. This procedure should create an object file named 'hpxdrc.o' (or the appropriate name for your

system) which can be included in the compilation of your Fortran program to provide a link to the system XDR

routines.

We are providing an example of a FORTRAN program (readxdr.f) in the programs/xdr area on the Final Report

CD to extract the data from the files on the model output CDROMs and create either time series of data from a

single point, or the entire gridded fields for each time step. You will have to include the XDR linkage routines

when you compile this program. On our system, the compilation command is 'f77 -o readxdr readxdr.f hpxdrc.o'.

The data are stored in the XDR files only for the lake (non-land) grid squares in the 5 km hydrodynamic grid.

There are 2318 such squares in the 53 x 102 grid. The example program includes the subroutines needed to map

the data back into the 2-dimensional hydrodynamic grid. We are also providing a program (gridinfo.f) to list the

geographic coordinates (latitude and longitude) of each grid square, the 5 km grid itself (depth at each grid

square), and a vector representation of the shoreline (lat-lon) that you can use in your plotting programs. The grid

and shoreline are in the file mich5.dat.